!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Structured Markov Chains Solver      [  SMCSolver  ]       !
!  Dario Bini, Beatrice Meini, Sergio Steffe'                 !
!  bini@dm.unipi.it, meini@dm.unipi.it, steffe@dm.unipi.it    !
!  Dipartimento di Matematica "Leonida Tonelli"               !
!  Largo Pontecorvo 5                                         !
!  56127 Pisa                                                 !
!  Italy                                                      !
!  Version 2.1 - June  2009                                   !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   fcalls.f90   main calls from C to  Fortran computation    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! erroflag =0 no errors, 1000 allocation errors, 2000 numeric
!    (lapack) errors  > 3000 various (5000= non stocastic )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !        globals in ponte_f_f.mod
  ! real(KIND=8),dimension(:,:,:),allocatable:: A,B
  ! real(KIND=8),dimension(:,:),allocatable:: G,R,U,Pi
  ! real(KIND=8),dimension(:,:),allocatable:: Pi0
  ! logical debug  !if true debugging output to std out 
  ! logical verb   !if true more detailed output to NotePad
  ! character (len=256) wout  !line of output
  !
  ! integer,parameter :: dp=kind(0.d0)
  ! integer :: info=0 info /=0 will catch allocation errors
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem QBD Algorithm Cyclic Reduction
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine qbd_cr_solve(algflag, goalflag,errorflag,errmax,itermax)
use ponte_f_f
use smc_int
implicit none
integer algflag   ! input:   0=basic,  1=shift acceleration,  2=tau shift 3=duoble shift  4=diagonal adjustment
integer goalflag  ! input:   1= G only, 2= G and R, 3= G, R and U,  4= G, R, U and Pi
integer errorflag ! output:  0= no error else some error did force premature termination
integer itermax   ! input:   max number of allowed iterations, default=50
integer errmax    ! input:   eps = epsilon(1.0d0)* 10 ** errmax
logical ::   dg  ! algorithm flags
integer ::  ds      ! algorithm flags
real(dp)::drift 
integer m         ! for allocating matrices
!
if (debug) write (*,*) "qbd_cr_solve: entering with algflag=",algflag, &
          " goalflag=",goalflag," errorflag=",errorflag,"eps=",errmax,"niter=",itermax
!
errorflag=0
finish=.false.
info =0; fdrift=0
!
!
! deallocate G,R,U,Pi,Pi0 if allocated
!
if (allocated(G)) deallocate(G)
if (allocated(R)) deallocate(R)
if (allocated(U)) deallocate(U)
if (allocated(Pi)) deallocate(Pi)
if (allocated(Pi0)) deallocate(Pi0)
!
!set the algorithm flags
!
select case (algflag)
case (0) !basic no shift
   ds=algflag
   dg=.false.
case (1) !shift
   ds=algflag
   dg=.false. 
case (2) !tau shift 
   ds=algflag
   dg=.false.
case (3) !double shift 
   ds=algflag
   dg=.false.
case (4) !diagonal adjustment
   ds=0  !no shift
   dg=.true.
case default
   ds=0
   dg=.false.
end select

m=size(A,1) 

select case(goalflag)
case (1) ! G only
     allocate(G(m,m),STAT=info)
     if (debug) write (*,*) "qbd_cr_solve:  Allocating G, stat=",info
     if (info /= 0 ) then
     errorflag=1000
     return
     endif
     call crqbd(am1=A(:,:,1),a0=A(:,:,2),a1=A(:,:,3),G=G,doshift=ds,&
dogth=dg,drift=drift, nerror=errmax ,maxit=itermax)
     write(wout,*)"drift=",drift
     call print_it
     if (info /= 0 ) then
     if (debug) write (*,*) "qbd_cr_solve:  allocation error executing crqbd, G Goal"
     errorflag=3000
     return
     end if
     if (debug) write (*,*) "qbd_cr_solve:  Cyclic Reduction  did find G"
     errorflag=info
     
     finish=.true.
case (2) ! G and R
     allocate(G(m,m),STAT=info)
     if (debug) write (*,*) "qbd_cr_solve:  Allocating G, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     allocate(R(m,m),STAT=info)
     if (debug) write (*,*) "qbd_cr_solve:  Allocating R, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     call crqbd(am1=A(:,:,1),a0=A(:,:,2),a1=A(:,:,3),G=G,R=R,doshift=ds,&
dogth=dg,drift=drift, nerror=errmax ,maxit=itermax)
     write(wout,*)"drift=",drift
     call print_it
     if (info /= 0 ) then
     if (debug) write (*,*) "qbd_cr_solve:  allocation error executing crqbd, G and R Goal"
     errorflag=info
     return
     endif
     if (debug) write (*,*) "qbd_cr_solve:  Cyclic Reduction  did find G and R"
     errorflag=info
     finish=.true.
case (3) ! G R and U
     allocate(G(m,m),STAT=info)
     if (debug) write (*,*) "qbd_cr_solve:  Allocating G, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     allocate(R(m,m),STAT=info)
     if (debug) write (*,*) "qbd_cr_solve:  Allocating R, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     allocate(U(m,m),STAT=info)
     if (debug) write (*,*) "qbd_cr_solve:  Allocating U, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     call crqbd(am1=A(:,:,1),a0=A(:,:,2),a1=A(:,:,3),G=G,R=R,U=U,doshift=ds,&
dogth=dg,drift=drift, nerror=errmax ,maxit=itermax)
     write(wout,*)"drift=",drift
     call print_it
     if (info /= 0 ) then
     if (debug) write (*,*) "qbd_cr_solve:  allocation error executing crqbd, G, R and U Goal"
     errorflag=info
     return
     endif
     if (debug) write (*,*) "qbd_cr_solve:  Cyclic Reduction  did find G, R and U"
     errorflag=info
     finish=.true. 
case (4) ! G R U and Pi
     allocate(G(m,m),STAT=info)
     if (debug) write (*,*) "qbd_cr_solve:  Allocating G, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     allocate(R(m,m),STAT=info)
     if (debug) write (*,*) "qbd_cr_solve:  Allocating R, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     allocate(U(m,m),STAT=info)
     if (debug) write (*,*) "qbd_cr_solve:  Allocating U, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     call crqbd(am1=A(:,:,1),a0=A(:,:,2),a1=A(:,:,3),G=G,R=R,U=U,doshift=ds,&
dogth=dg,drift=drift, nerror=errmax ,maxit=itermax)
     write(wout,*)"drift=",drift
     call print_it
      if (info /= 0 ) then
        if (debug) write (*,*) "qbd_cr_solve:  allocation error executing crqbd, G, R ,U and PiGoal"
        errorflag=info
        return
     endif
     if (debug) write (*,*) "qbd_cr_solve:  Cyclic Reduction  did find G, R and U"
     errorflag=info
     finish=.true.
     fdrift=drift
case default
if(debug) write(*,*) "qbd_cr_solve: strange goalflag=",goalflag
end select
if (debug) write (*,*) "qbd_cr_solve: exiting with algflag=",algflag,&
" goalflag=",goalflag,"errorflag=",errorflag," finish =",finish
end subroutine qbd_cr_solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem QBD Algorithm Logarithmic Reduction
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine qbd_lr_solve(algflag, goalflag,errorflag,errmax,itermax)
use ponte_f_f
use smc_int
implicit none
integer algflag   ! input:   0=basic,  1=shift acceleration, 2=tau shift, 3=double shift,  4= diagonal adjustment
integer goalflag  ! input:   1= G only, 2= G and R, 3= G, R and U,  4= G, R, U and Pi
integer errorflag ! output:  0= no error else some error did force premature termination
integer itermax   ! input:   max number of allowed iterations, default=50
integer errmax    ! input:   eps = epsilon(1.0d0)* 10 ** errmax
logical ::  dg  ! algorithm flags
integer ::  ds ! algorithm flags
real(dp)::drift    
integer m         
!
!
if (debug) write (*,*) "qbd_lr_solve: entering with algflag=",algflag,&
        " goalflag=",goalflag," errorflag=",errorflag,"eps=",errmax,"niter=",itermax
!
errorflag=0
finish=.false.
info =0; fdrift=0
!
! deallocate G,R,U,Pi,Pi0 if allocated
!
if (allocated(G)) deallocate(G)
if (allocated(R)) deallocate(R)
if (allocated(U)) deallocate(U)
if (allocated(Pi)) deallocate(Pi)
if (allocated(Pi0)) deallocate(Pi0)
!
!set the algorithm flags
select case (algflag)
case (0) !basic no shift
   ds=algflag
   dg=.false.
case (1) !shift
   ds=algflag
   dg=.false. 
case (2) !tau shift 
   ds=algflag
   dg=.false.
case (3) !double shift 
   ds=algflag
   dg=.false.
case (4) !diagonal adjustment
   ds=0  !no shift
   dg=.true.
case default
   ds=0
   dg=.false.
end select


m=size(A,1)

select case(goalflag)
case (1) ! G only
     allocate(G(m,m),STAT=info)
     if (debug) write (*,*) "qbd_lr_solve:  Allocating G, stat=",info
     if (info /= 0 ) then 
     errorflag=info
     return
     endif
     call lrqbd(am1=A(:,:,1),a0=A(:,:,2),a1=A(:,:,3),G=G,doshift=ds,&
dogth=dg,drift=drift, nerror=errmax ,maxit=itermax)
     write(wout,*)"drift=",drift
     call print_it
      if (debug) write (*,*) "qbd_lr_solve:  Logarithmic Reduction  did find G"
     errorflag=info
     finish=.true.
case (2) ! G and R
     allocate(G(m,m),STAT=info)
     if (debug) write (*,*) "qbd_lr_solve:  Allocating G, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     allocate(R(m,m),STAT=info)
     if (debug) write (*,*) "qbd_lr_solve:  Allocating R, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     call lrqbd(am1=A(:,:,1),a0=A(:,:,2),a1=A(:,:,3),G=G,R=R,doshift=ds,&
dogth=dg,drift=drift, nerror=errmax ,maxit=itermax)
     write(wout,*)"drift=",drift
     call print_it
      if (debug) write (*,*) "qbd_lr_solve:  Logarithmic Reduction  did find G and R"
     errorflag=info
     finish=.true.
case (3) ! G R and U
     allocate(G(m,m),STAT=info)
     if (debug) write (*,*) "qbd_lr_solve:  Allocating G, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     allocate(R(m,m),STAT=info)
     if (debug) write (*,*) "qbd_lr_solve:  Allocating R, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     allocate(U(m,m),STAT=info)
     if (debug) write (*,*) "qbd_lr_solve:  Allocating U, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     call lrqbd(am1=A(:,:,1),a0=A(:,:,2),a1=A(:,:,3),G=G,R=R,U=U,doshift=ds,&
dogth=dg,drift=drift, nerror=errmax ,maxit=itermax)
     write(wout,*)"drift=",drift
     call print_it
      if (debug) write (*,*) "qbd_lr_solve:  Logarithmic Reduction  did find G, R and U"
     errorflag=info
     finish=.true.
case (4) ! G R U and Pi
     allocate(G(m,m),STAT=info)
     if (debug) write (*,*) "qbd_lr_solve:  Allocating G, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     allocate(R(m,m),STAT=info)
     if (debug) write (*,*) "qbd_lr_solve:  Allocating R, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     allocate(U(m,m),STAT=info)
     if (debug) write (*,*) "qbd_lr_solve:  Allocating U, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     call lrqbd(am1=A(:,:,1),a0=A(:,:,2),a1=A(:,:,3),G=G,R=R,U=U,doshift=ds,&
dogth=dg,drift=drift, nerror=errmax ,maxit=itermax)
     write(wout,*)"drift=",drift
     call print_it
      if (debug) write (*,*) "qbd_lr_solve:  Logarithmic Reduction  did find G, R and U"
     errorflag=info
     finish=.true.
     if(debug) write(*,*) "qbd_lr_solve: computation of G R U and Pi not yet complete ..."
case default
if(debug) write(*,*) "qbd_lr_solve: strange goalflag=",goalflag
end select

fdrift=drift
if (debug) write (*,*) "qbd_lr_solve: exiting with algflag=",algflag,&
" goalflag=",goalflag,"errorflag=",errorflag," finish =",finish

end subroutine qbd_lr_solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem QBD Algorithm Functional Iteration
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine qbd_fi_solve(algflag, goalflag,errorflag,errmax,itermax)
use ponte_f_f
use fi_int
use smc_int
implicit none
integer algflag   ! input:   100=natural 200=traditional 300=U-based  
                  !             +1=Shift +2=TauShift +3=Double Shift
                  !             +30=x(0)=0 +20=x(0)=I +10=x(0)=user X
integer goalflag  ! input:   1= G only, 2= G and R, 3= G, R and U,  4= G, R, U and Pi
integer errorflag ! output:  0= no error else some error did force premature termination
integer itermax   ! input:   max number of allowed iterations, default=50
integer errmax    ! input:   eps = epsilon(1.0d0)* 10 ** errmax
real(dp)::drift,norm1,myeps     !
integer m
real(dp),dimension(:),allocatable :: v
! logical ::  ds    ! do shift
integer :: method ! 1 natural, 2 traditional, 3 U-based (default)
integer :: i
!
!
!
if (debug) write (*,*) "qbd_fi_solve: entering with algflag=",algflag,&
        " goalflag=",goalflag," errorflag=",errorflag," eps=",errmax," niter=",itermax
!
errorflag=0
finish=.false.
info =0; fdrift=0
myeps=epsilon(1.d0) * 10.0d0 ** errmax
!
!
! deallocate G,R,U,Pi,Pi0 if allocated
!
if (allocated(G)) deallocate(G)
if (allocated(R)) deallocate(R)
if (allocated(U)) deallocate(U)
if (allocated(Pi)) deallocate(Pi)
if (allocated(Pi0)) deallocate(Pi0)
if (allocated(v)) deallocate(v)

m=size(A,1)
!
!
allocate(G(m,m),v(m),STAT=info)
if (debug) write (*,*) "qbd_fi_solve:  Allocating G, v, stat=",info
if (info /= 0 ) then
     errorflag=info
     return
endif

call drft(am1=a(:,:,1),a0=a(:,:,2),a1=a(:,:,3),v=v,drift=drift)
fdrift=drift

method=3   
!
!setup flags
!
if(algflag >= 300) then !U-based (default)
   method=3 
   algflag=algflag-300
endif
if(algflag >= 200) then !Traditional
   method=2 
   algflag=algflag-200
endif
if(algflag >= 100) then !Natural
   method=1
   algflag=algflag-100
endif

if(algflag .GE. 30) then !  x(0)=0 (default)
      if (allocated(X0)) deallocate(X0)
      allocate(X0(m,m),STAT=info)
      if (debug) write (*,*) "qbd_fi_solve:  Allocating X0, stat=",info
      if (info /= 0 ) then
          errorflag=info
          return
      endif
      X0=0.0d0
      algflag=algflag-30
endif
if(algflag .GE. 20) then !  x(0)=I 
      if (allocated(X0)) deallocate(X0)
      allocate(X0(m,m),STAT=info)
      if (debug) write (*,*) "qbd_fi_solve:  Allocating X0, stat=",info
      if (info /= 0 ) then
          errorflag=info
          return
      endif
      X0=0.0d0
      do i=1,m
      X0(i,i)=1.0d0
      enddo
      algflag=algflag-20
endif
if (algflag .GE. 10) then !user defined x(0)
   algflag=algflag-10
   if (debug) write (*,*) "qbd_fi_solve: used defined x(0)"
endif
!
! setup goals
!
call fi(A=A, ds=algflag, method=method, eps=myeps, maxit=itermax, x0=X0, G=G, drift=drift, err=norm1)
if ( .not. verb) then
   write(wout,*)" "
   call print_it
endif
write(wout,*)"G-residual=",norm1
call print_it
select case(goalflag)
case(2) ! G and R
allocate(u(m,m),r(m,m),stat=info)
if(info/=0)then
write(wout,*)"qbd_fi_solve: info=",info
call print_it
end if
call gtou(A0=A(:,:,2),A1=A(:,:,3),G=G,U=U)
call utor(A1=A(:,:,3),u=u,r=r)
call qbdrres(am1=a(:,:,1),a0=a(:,:,2),a1=a(:,:,3),r=r,res=norm1)
write(wout,*)"R-residual=",norm1
call print_it
case(3) ! G  R U
allocate(u(m,m),r(m,m),stat=info)
if(info/=0)then
write(wout,*)"qbd_fi_solve: info=",info
call print_it
end if
call gtou(A0=A(:,:,2),A1=A(:,:,3),G=g,U=u)
call utor(A1=A(:,:,3),u=u,r=r)
call qbdrres(am1=a(:,:,1),a0=a(:,:,2),a1=a(:,:,3),r=r,res=norm1)
write(wout,*)"R-residual=",norm1
call print_it
call qbdures(am1=a(:,:,1),a0=a(:,:,2),a1=a(:,:,3),u=u,res=norm1)
write(wout,*)"U-residual=",norm1
call print_it
case(4) ! pi G  R U
allocate(u(m,m),r(m,m),stat=info)
if(info/=0)then
write(wout,*)"qbd_fi_solve: info=",info
call print_it
end if
call gtou(A0=A(:,:,2),A1=A(:,:,3),G=g,U=u)
call utor(A1=A(:,:,3),u=u,r=r)
call qbdrres(am1=a(:,:,1),a0=a(:,:,2),a1=a(:,:,3),r=r,res=norm1)
write(wout,*)"R-residual=",norm1
call print_it
call qbdures(am1=a(:,:,1),a0=a(:,:,2),a1=a(:,:,3),u=u,res=norm1)
write(wout,*)"U-residual=",norm1
call print_it
end select
write(wout,*)"drift=",fdrift
call print_it

!!
if (debug) write (*,*) "qbd_fi_solve: exiting with algflag=",algflag,&
" goalflag=",goalflag,"errorflag=",errorflag," finish =",finish
end subroutine qbd_fi_solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem QBD Algorithm Invariant Subspace
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine qbd_is_solve(algflag, goalflag,errorflag,errmax,itermax)
use ponte_f_f
use is_int
use smc_int
implicit none
integer algflag   ! input:   0=basic, 1=Balzer 2=Shur
integer goalflag  ! input:   1= G only, 2= G and R, 3= G, R and U,  4= G, R, U and Pi
integer errorflag ! output:  0= no error else some error did force premature termination
integer itermax   ! input:   max number of allowed iterations, default=50
integer errmax    ! input:   eps = epsilon(1.0d0)* 10 ** errmax
real(dp)::drift,norm1,myeps     !
real(dp),dimension(:),allocatable :: v
integer m  
!
!
!
if (debug) write (*,*) "qbd_is_solve: entering with algflag=",algflag,&
         " goalflag=",goalflag," errorflag=",errorflag," eps=",errmax," niter=",itermax
!
errorflag=0
finish=.false.
info =0; fdrift=0
myeps=epsilon(1.d0) * 10.0d0 ** errmax
!
!
! deallocate G,R,U,Pi,Pi0 if allocated
!
if (allocated(G)) deallocate(G)
if (allocated(R)) deallocate(R)
if (allocated(U)) deallocate(U)
if (allocated(Pi)) deallocate(Pi)
if (allocated(Pi0)) deallocate(Pi0)
if (allocated(v)) deallocate(v)

m=size(A,1)
     allocate(G(m,m),v(m),STAT=info)
     if (debug) write (*,*) "qbd_is_solve:  Allocating G, v, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif

call drft(am1=a(:,:,1),a0=a(:,:,2),a1=a(:,:,3),v=v,drift=drift)
fdrift=drift

call is(A, method=algflag, eps=myeps, maxit=itermax, G=G, drift=drift, err=norm1);

if (debug) write (*,*) "qbd_is_solve: exit with info=",info

if ( .not. verb) then
   write(wout,*)" "
   call print_it
endif
write(wout,*)"G-residual=",norm1
call print_it

if (info/=0) then
  write(wout,*)"qbd_is_solve: ERROR, INFO=",info," exiting now !"
  call print_it
  errorflag=info
  return
endif

select case(goalflag)
case(2) ! G and R
allocate(u(m,m),r(m,m),stat=info)
if(info/=0)then
write(wout,*)"qbd_is_solve: info=",info
call print_it
end if
call gtou(A0=A(:,:,2),A1=A(:,:,3),G=G,U=U)
call utor(A1=A(:,:,3),u=u,r=r)
call qbdrres(am1=a(:,:,1),a0=a(:,:,2),a1=a(:,:,3),r=r,res=norm1)
write(wout,*)"R-residual=",norm1
call print_it
case(3) ! G  R U
allocate(u(m,m),r(m,m),stat=info)
if(info/=0)then
write(wout,*)"qbd_is_solve: info=",info
call print_it
end if
call gtou(A0=A(:,:,2),A1=A(:,:,3),G=g,U=u)
call utor(A1=A(:,:,3),u=u,r=r)
call qbdrres(am1=a(:,:,1),a0=a(:,:,2),a1=a(:,:,3),r=r,res=norm1)
write(wout,*)"R-residual=",norm1
call print_it
call qbdures(am1=a(:,:,1),a0=a(:,:,2),a1=a(:,:,3),u=u,res=norm1)
write(wout,*)"U-residual=",norm1
call print_it
case(4) ! pi G  R U
allocate(u(m,m),r(m,m),stat=info)
if(info/=0)then
write(wout,*)"qbd_is_solve: info=",info
call print_it
end if
call gtou(A0=A(:,:,2),A1=A(:,:,3),G=g,U=u)
call utor(A1=A(:,:,3),u=u,r=r)
call qbdrres(am1=a(:,:,1),a0=a(:,:,2),a1=a(:,:,3),r=r,res=norm1)
write(wout,*)"R-residual=",norm1
call print_it
call qbdures(am1=a(:,:,1),a0=a(:,:,2),a1=a(:,:,3),u=u,res=norm1)
write(wout,*)"U-residual=",norm1
call print_it
end select
 write(wout,*)"drift=",fdrift
 call print_it


if (debug) write (*,*) "qbd_is_solve: exiting with algflag=",algflag,&
" goalflag=",goalflag,"errorflag=",errorflag," finish =",finish
end subroutine qbd_is_solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem M/G/1 Algorithm Cyclic Reduction
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine mg1_cr_solve(algflag, goalflag,errorflag,errmax,itermax,interpmax,dual)
  use ponte_f_f
  use pwcr_interface
  use smc_int

  implicit none
  integer algflag   ! input:   0=basic,  1=shift acceleration, 2= tau shift 3= double shift  4= diagonal adjustment
  integer goalflag  ! input:   1= G only, 2= G and Pi
  integer errorflag ! output:  0= no error else some error did force premature termination
  integer itermax   ! input:   max number of allowed iterations, default=50
  integer errmax    ! input:   eps = epsilon(1.0d0)* 10 ** errmax
  integer interpmax ! input:   man number of allowed interpolations, default=1  #interpolations= 256*2**interpmax
  integer dual      ! input:   0=automatic, 1=ramaswami, 2=brigth 3= both 4=no
  real(dp)::norm1,myeps,tau     
  integer:: m, nba,ldual,i
  real(kind(0.d0)),dimension(:,:,:),allocatable:: aa
  real(kind(0.d0)),dimension(:),allocatable:: v,vi
  !
  !
  !

  if (debug) write (*,*) "mg1_cr_solve: entering with algflag=",algflag,&
       " goalflag=",goalflag," errorflag=",errorflag," eps=",errmax," niter=",&
       itermax," interpmax=",interpmax," dual=",dual
  !
  errorflag=0
  finish=.false.
  info =0; fdrift=0
  myeps=epsilon(1.d0) * 10.0d0 ** errmax
  !
  !
  ! deallocate G,R,U,Pi,Pi0 if allocated
  !
  if (allocated(G)) deallocate(G)
  if (allocated(R)) deallocate(R)
  if (allocated(U)) deallocate(U)
  if (allocated(Pi)) deallocate(Pi)
  if (allocated(Pi0)) deallocate(Pi0)
  if (allocated(AA)) deallocate(AA)
  if (allocated(v)) deallocate(v)
  if (allocated(vi)) deallocate(vi)

  m=size(A,1)
  nba=size(A,3)
  allocate(AA(m,m,nba),STAT=info)
  if (debug) write (*,*) "mg1_cr_solve:  Allocating AA, stat=",info
  AA=A

  ldual=dual ! local dual

  do i=1,nba
     aa(:,:,i)=transpose(aa(:,:,i))
  end do
  fdrift=drftmg1(aa)
  do i=1,nba
     aa(:,:,i)=transpose(aa(:,:,i))
  end do

  if(ldual==0) then
     if (fdrift.ge.0) then ! M/G/1 transient
        ldual=3 !both duals
     else
        ldual=4 !no dual
     end if
  end if

  if(ldual.eq.3)then 
     allocate(v(m),vi(m),STAT=info)
     if (debug) write (*,*) "mg1_cr_solve:  Allocating w, stat=",info
     do i=1,nba
        aa(:,:,i)=transpose(a(:,:,i))
     end do
     call mg1tomg1(AA, AA,tau, v)
     do i=1,nba
        aa(:,:,i)=transpose(aa(:,:,i))
     end do
  end if

  select case(goalflag)
  case (1) ! G only
     allocate(G(m,m),STAT=info)
     if (debug) write (*,*) "mg1_cr_solve:  Allocating G, stat=",info
     if (info /= 0 ) then
        errorflag=info
        return
     endif
     if (algflag==0) call  pwcr(A=AA,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==1) call  spwcr(A=AA,ds=1,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax) 
     if (algflag==2) call  spwcr(A=AA,ds=2,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==3) call  spwcr(A=AA,ds=3,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==4) call  pwcr(A=AA,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax) 
     if (info /= 0 ) then
        errorflag=info
        if (info == 5000) then
           write(wout,*) "===> A FAILS Stochaticity check: Algorithm max error bound exceeded"
           call print_it
        endif
        return
     endif
     if ( .not. verb) then
        write(wout,*)" "
        call print_it
     endif
!!!!!!!! ldual=3 !!!!!!!!!!!!!!
     if (ldual.eq.3) then
        vi=1.d0/v
        do i=1,m
           G(:,i)=G(:,i)*vi(i)
           G(i,:)=G(i,:)*v(i)
        end do
        G=G*tau
        g=transpose(g)
        call residual(a,g,norm1)
        g=transpose(g)
     end if
!!!!!!!!!!!!!!!!!!!!
     if(dual.eq.0)then
        if (ldual.eq.3) then
           write(wout,*)"Automatic dual: both duals"
           call print_it
        end if
        if (ldual.eq.4) then
           write(wout,*)"Automatic dual: no dual"
           call print_it
        end if
     end if

     write(wout,*)"G-residual=",norm1
     call print_it

     write(wout,*)"drift=",fdrift
     call print_it


     if (debug) write (*,*) "mg1_cr_solve:  Cyclic Reduction  did find G"
     errorflag=info
     finish=.true.
  case (2) ! G and Pi
     allocate(G(m,m),STAT=info)
     if (debug) write (*,*) "mg1_cr_solve:  Allocating G, stat=",info
     if (info /= 0 ) then
        errorflag=info
        return
     endif
     if (algflag==0) call  pwcr(A=AA,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax) 
     if (algflag==1) call  spwcr(A=AA,ds=1,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==2) call  spwcr(A=AA,ds=2,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==3) call  spwcr(A=AA,ds=3,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==4) call  pwcr(A=AA,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (info /= 0 ) then
        errorflag=info
        if (info == 5000) then
           write(wout,*) "===> A FAILS Stochaticity check: Algorithm max error bound exceeded"
           call print_it
        endif
        return
     endif
     if (.not. verb) then 
        write(wout,*)" "
        call print_it
     endif
!!!!!!!!! ldual=3 !!!!!!!!!!!!!!!!!!!!!!!!!!
     if (ldual.eq.3) then
        vi=1.d0/v
        do i=1,m
           G(:,i)=G(:,i)*vi(i)
           G(i,:)=G(i,:)*v(i)
        end do
        G=G*tau
        g=transpose(g)
        call residual(a,g,norm1)
        g=transpose(g)
     end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     if(dual.eq.0)then
        if (ldual.eq.3) then
           write(wout,*)"Automatic dual: both duals"
           call print_it
        end if
        if (ldual.eq.4) then
           write(wout,*)"Automatic dual: no dual"
           call print_it
        end if
     end if

     write(wout,*)"G-residual=",norm1
     call print_it
     write(wout,*)"drift=",fdrift
     call print_it

     if (debug) write (*,*) "mg1_cr_solve:  Cyclic Reduction  did find G"
     errorflag=info
     finish=.true.
  case default
     if(debug) write(*,*) "mg1_cr_solve: strange goalflag=",goalflag
  end select

  if (debug) write (*,*) "mg1_cr_solve: exiting with algflag=",algflag,&
       " goalflag=",goalflag,"errorflag=",errorflag," finish =",finish," dual=",dual
end subroutine mg1_cr_solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem M/G/1 Algorithm Functional Iteration
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine mg1_fi_solve(algflag, goalflag,errorflag,errmax,itermax,dual)
use ponte_f_f
use fi_int
use smc_int
use pwcr_interface
implicit none
integer algflag   ! input:   100=natural 200=traditional 300=U-based  
                  !             +1=Shift +2=TauShift +3=Double Shift
                  !             +30=x(0)=0 +20=x(0)=I +10=x(0)=user X
integer goalflag  ! input:   1= G only, 2= G  and Pi
integer errorflag ! output:  0= no error else some error did force premature termination
integer itermax   ! input:   max number of allowed iterations, default=50
integer errmax    ! input:   eps = epsilon(1.0d0)* 10 ** errmax
integer dual      ! input:   0=automatic, 1=ramaswami, 2=brigth 3= both 4=no
real(dp)::drift,norm1,myeps,tau     !
integer m
!logical ::  ds    ! do shift
integer :: method ! 1 natural, 2 traditional, 3 U-based (default)
integer :: i, ldual,nba
  real(kind(0.d0)),dimension(:,:,:),allocatable:: aa
  real(kind(0.d0)),dimension(:),allocatable:: v,vi
!
!
!
if (debug) write (*,*) "mg1_fi_solve: entering with algflag=",algflag,&
           " goalflag=",goalflag," errorflag=",errorflag," eps=",errmax,&
           " niter=",itermax," dual=",dual
!
errorflag=0
finish=.false.
info =0; fdrift=0
myeps=epsilon(1.d0) * 10.0d0 ** errmax
!
!
! deallocate G,R,U,Pi,Pi0 if allocated
!
if (allocated(G)) deallocate(G)
if (allocated(R)) deallocate(R)
if (allocated(U)) deallocate(U)
if (allocated(Pi)) deallocate(Pi)
if (allocated(Pi0)) deallocate(Pi0)
  if (allocated(AA)) deallocate(AA)
  if (allocated(v)) deallocate(v)
  if (allocated(vi)) deallocate(vi)
!
m=size(A,1)
  nba=size(A,3)
  allocate(AA(m,m,nba),STAT=info)
  if (debug) write (*,*) "fi_cr_solve:  Allocating AA, stat=",info
  AA=A
!
!
allocate(G(m,m),STAT=info)
if (debug) write (*,*) "mg1_fi_solve:  Allocating G, stat=",info
if (info /= 0 ) then
     errorflag=info
     return
endif


  ldual=dual ! local dual
  fdrift=drftmg1(a)
  if(ldual==0) then
     if (fdrift.ge.0) then ! M/G/1 transient
        ldual=3 !both duals
     else
        ldual=4 !no dual
     end if
  end if

  if(ldual.eq.3)then 
     allocate(v(m),vi(m),STAT=info)
     if (debug) write (*,*) "mg1_fi_solve:  Allocating v, stat=",info
     call mg1tomg1(A, AA,tau, v)
  end if

method=3   
!
!setup flags
!
if(algflag >= 300) then !U-based (default)
   method=3 
   algflag=algflag-300
endif
if(algflag >= 200) then !Traditional
   method=2 
   algflag=algflag-200
endif
if(algflag >= 100) then !Natural
   method=1
   algflag=algflag-100
endif

if(algflag .GE. 30) then !  x(0)=0 (default)
      if (allocated(X0)) deallocate(X0)
      allocate(X0(m,m),STAT=info)
      if (debug) write (*,*) "mg1_fi_solve:  Allocating X0, stat=",info
      if (info /= 0 ) then
          errorflag=info
          return
      endif
      X0=0.0d0
      algflag=algflag-30
endif
if(algflag .GE. 20) then !  x(0)=I 
      if (allocated(X0)) deallocate(X0)
      allocate(X0(m,m),STAT=info)
      if (debug) write (*,*) "mg1_fi_solve:  Allocating X0, stat=",info
      if (info /= 0 ) then
          errorflag=info
          return
      endif
      X0=0.0d0
      do i=1,m
      X0(i,i)=1.0d0
      enddo
      algflag=algflag-20
endif
if (algflag .GE. 10) then !user defined x(0)
   algflag=algflag-10
   if (debug) write (*,*) "mg1_fi_solve: used defined x(0)"
endif
!
! setup goals
! only G for now
!
call fi(A=AA, ds=algflag, method=method, eps=myeps, maxit=itermax, x0=X0, G=G, drift=drift, err=norm1)
if ( .not. verb) then
       write(wout,*)" "
       call print_it
     endif
!!!!!!!! ldual=3 !!!!!!!!!!!!!!
     if (ldual.eq.3) then
        vi=1.d0/v
        do i=1,m
           G(:,i)=G(:,i)*vi(i)
           G(i,:)=G(i,:)*v(i)
        end do
        G=G*tau
        do i=1,nba
           aa(:,:,i)=transpose(a(:,:,i))
        end do
        call gresidual(a,g,norm1)
     end if
!!!!!!!!!!!!!!!!!!!!
     if(dual.eq.0)then
        if (ldual.eq.3) then
           write(wout,*)"Automatic dual: both duals"
           call print_it
        end if
        if (ldual.eq.4) then
           write(wout,*)"Automatic dual: no dual"
           call print_it
        end if
     end if

write(wout,*)"G-residual=",norm1
call print_it
write(wout,*)"drift=",fdrift
 call print_it

!!
if (debug) write (*,*) "mg1_fi_solve: exiting with algflag=",algflag,&
" goalflag=",goalflag,"errorflag=",errorflag," finish =",finish," dual=",dual
end subroutine mg1_fi_solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem M/G/1 Algorithm Invariant Subspace
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine mg1_is_solve(algflag, goalflag,errorflag,errmax,itermax,dual)
use ponte_f_f
use is_int
use fi_int
use smc_int
implicit none
integer algflag   ! input:   0=basic, 1=Balzer 2=Shur
integer goalflag  ! input:   1= G only, 2= G and Pi
integer errorflag ! output:  0= no error else some error did force premature termination
integer itermax   ! input:   max number of allowed iterations, default=50
integer errmax    ! input:   eps = epsilon(1.0d0)* 10 ** errmax
integer dual      ! input:   0=automatic, 1=ramaswami, 2=brigth 3= both 4=no
real(dp)::drift,norm1,myeps,tau    !
integer i,m, ldual,nba
 real(kind(0.d0)),dimension(:,:,:),allocatable:: aa
  real(kind(0.d0)),dimension(:),allocatable:: v,vi
!
!
!
if (debug) write (*,*) "mg1_is_solve: entering with algflag=",algflag,&
        " goalflag=",goalflag," errorflag=",errorflag," eps=",errmax," niter=",itermax," dual=",dual
!
errorflag=0
finish=.false.
info =0; fdrift=0
myeps=epsilon(1.d0) * 10.0d0 ** errmax
!
!
! deallocate G,R,U,Pi,Pi0 if allocated
!
if (allocated(G)) deallocate(G)
if (allocated(R)) deallocate(R)
if (allocated(U)) deallocate(U)
if (allocated(Pi)) deallocate(Pi)
if (allocated(Pi0)) deallocate(Pi0)
  if (allocated(AA)) deallocate(AA)
  if (allocated(v)) deallocate(v)
  if (allocated(vi)) deallocate(vi)

m=size(A,1)
  nba=size(A,3)
  allocate(AA(m,m,nba),STAT=info)
  if (debug) write (*,*) "is_cr_solve:  Allocating AA, stat=",info
  AA=A

     allocate(G(m,m),STAT=info)
     if (debug) write (*,*) "mg1_is_solve:  Allocating G, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif

  ldual=dual ! local dual
  fdrift=drftmg1(a)
  if(ldual==0) then
     if (fdrift.ge.0) then ! M/G/1 transient
        ldual=3 !both duals
     else
        ldual=4 !no dual
     end if
  end if

  if(ldual.eq.3)then 
     allocate(v(m),vi(m),STAT=info)
     if (debug) write (*,*) "mg1_fi_solve:  Allocating v, stat=",info
      call mg1tomg1(A, AA,tau, v)
  end if

call is(A=AA, method=algflag, eps=myeps, maxit=itermax, G=G, drift=drift, err=norm1)

!!!!!!!! ldual=3 !!!!!!!!!!!!!!
     if (ldual.eq.3) then
        vi=1.d0/v
        do i=1,m
           G(:,i)=G(:,i)*vi(i)
           G(i,:)=G(i,:)*v(i)
        end do
        G=G*tau
        do i=1,nba
           aa(:,:,i)=transpose(a(:,:,i))
        end do
        call gresidual(a,g,norm1)
     end if
!!!!!!!!!!!!!!!!!!!!
     if(dual.eq.0)then
        if (ldual.eq.3) then
           write(wout,*)"Automatic dual: both duals"
           call print_it
        end if
        if (ldual.eq.4) then
           write(wout,*)"Automatic dual: no dual"
           call print_it
        end if
     end if

write(wout,*)"G-residual=",norm1
call print_it
write(wout,*)"drift=",fdrift
 call print_it

if (debug) write (*,*) "mg1_is_solve: exiting with algflag=",algflag," goalflag=",&
goalflag,"errorflag=",errorflag," finish =",finish," dual=",dual
end subroutine mg1_is_solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem GI/M/1 Algorithm Cyclic Reduction
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine gim1_cr_solve(algflag, goalflag,errorflag,errmax,itermax,interpmax,dual)
use ponte_f_f
use pwcr_interface
use smc_int 
implicit none
integer algflag   ! input:   0=basic,  1=shift acceleration, 2=tau shift 3= double shift  4= diagonal adjustment
integer goalflag  ! input:   1= R only, 2= R and Pi
integer errorflag ! output:  0= no error else some error did force premature termination
integer itermax   ! input:   max number of allowed iterations, default=50
integer interpmax ! input:   man number of allowed interpolations, default=1  #interpolations= 256*2**interpmax
integer errmax    ! input:   eps = epsilon(1.0d0)* 10 ** errmax
integer dual      ! input:   0=automatic, 1=ramaswami, 2=brigth 3= both
real(dp)::norm1,myeps, tau
integer m
integer:: nba,  i, ldual
 real(dp),dimension(:),allocatable :: vi,v
 real(dp),dimension(:,:,:),allocatable :: AA
if (debug) write (*,*) "gim1_cr_solve: entering with algflag=",algflag,&
        " goalflag=",goalflag," errorflag=",errorflag," eps=",errmax,&
        " niter=",itermax," interpmax=",interpmax," dual=",dual
errorflag=0
finish=.false.
info =0; fdrift=0
myeps=epsilon(1.d0) * 10.0d0 ** errmax
!
!
! deallocate G,R,U,Pi,Pi0 if allocated
!
if (allocated(G)) deallocate(G)
if (allocated(R)) deallocate(R)
if (allocated(AA)) deallocate(AA)
if (allocated(Pi)) deallocate(Pi)
if (allocated(Pi0)) deallocate(Pi0)
if (allocated(v)) deallocate(v)

m=size(A,1)
nba=size(A,3)
allocate(R(m,m),AA(m,m,nba),v(m),stat=info)
if(info/=0)then
   return
end if

ldual=dual ! local dual
fdrift=-drftmg1(a)

if(ldual==0) then
    if (fdrift.ge.0) then ! G/M/1 transient
      ldual=1 !ram dual
   else ! G/M/1 recurrent
      ldual=2 !bright dual
   end if
end if

if (ldual==1) call gm1tomg1(A,AA,v)
if (ldual==2) call gm1tomg1b(A, AA,tau, v)

do i=1,nba
   aa(:,:,i)=transpose(aa(:,:,i))
end do

select case(goalflag)
case (1) ! G only
     allocate(G(m,m),STAT=info)
     if (debug) write (*,*) "gm1_cr_solve:  Allocating G, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     if (algflag==0) call  pwcr(A=AA,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==1) call  spwcr(A=AA,ds=1,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==2) call  spwcr(A=AA,ds=2,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==3) call  spwcr(A=AA,ds=3,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==4) call  pwcr(A=AA,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax) 
     if (info /= 0 ) then
       errorflag=info
        if (info == 5000) then
          write(wout,*) "===> A FAILS Stochaticity check: Algorithm max error bound exceeded"
          call print_it
        endif
        return
     endif
     if ( .not. verb) then
       write(wout,*)" "
       call print_it
     endif
       r=transpose(g)
     if (allocated(vi)) deallocate(vi)
     allocate(vi(m),stat=info)
     if (info/=0)then
        return
     end if
     vi=1.d0/v
     do i=1,m
        R(:,i)=R(:,i)*v(i)
        R(i,:)=R(i,:)*vi(i)
     end do
     if (ldual==2) then
        r=r*tau
     end if
     if(dual.eq.0)then
        if (ldual.eq.1) then
           write(wout,*)"Automatic dual: Ramaswami dual"
           call print_it
        end if
        if (ldual.eq.2) then
           write(wout,*)"Automatic dual: Bright dual"
           call print_it
        end if
     end if

     call rresidual(A,R,norm1)
     write(wout,*)"R-residual=",norm1
     call print_it
     write(wout,*)"drift=",fdrift
     call print_it

     if (debug) write (*,*) "gm1_cr_solve:  Cyclic Reduction  did find G"
     errorflag=info
     finish=.true.
case (2) ! G and Pi
     allocate(G(m,m),STAT=info)
     if (debug) write (*,*) "gm1_cr_solve:  Allocating G, stat=",info
     if (info /= 0 ) then
     errorflag=info
     return
     endif
     if (algflag==0) call  pwcr(A=AA,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax) 
     if (algflag==1) call  spwcr(A=AA,ds=1,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==2) call  spwcr(A=AA,ds=2,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==3) call  spwcr(A=AA,ds=3,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==4) call  pwcr(A=AA,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (info /= 0 ) then
       errorflag=info
        if (info == 5000) then
          write(wout,*) "===> A FAILS Stochaticity check: Algorithm max error bound exceeded"
          call print_it
        endif
        return
     endif
     if (.not. verb) then 
       write(wout,*)" "
       call print_it
     endif
     r=transpose(g)
     if (allocated(vi)) deallocate(vi)
     allocate(vi(m),stat=info)
     if (info/=0)then
        return
     end if
     vi=1.d0/v
     do i=1,m
        R(:,i)=R(:,i)*v(i)
        R(i,:)=R(i,:)*vi(i)
     end do
     if (ldual==2) then
        r=r*tau
     end if
     if(dual.eq.0)then
        if (ldual.eq.1) then
           write(wout,*)"Automatic dual: Ramaswami dual"
           call print_it
        end if
        if (ldual.eq.2) then
           write(wout,*)"Automatic dual: Bright dual"
           call print_it
        end if
     end if

     call rresidual(A,r,norm1)
     write(wout,*)"R-residual=",norm1
     call print_it
     write(wout,*)"drift=",fdrift
     call print_it

     if (debug) write (*,*) "gm1_cr_solve:  Cyclic Reduction  did find G"
     errorflag=info
     finish=.true.
case default
if(debug) write(*,*) "gm1_cr_solve: strange goalflag=",goalflag
end select






if (debug) write (*,*) "gim1_cr_solve: exiting with algflag=",algflag,&
       " goalflag=",goalflag,"errorflag=",errorflag," finish =",finish," dual=",dual
end subroutine gim1_cr_solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem GI/M/1 Algorithm Logarithmic Reduction
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine gim1_lr_solve(algflag, goalflag,errorflag,errmax,itermax,dual)
use ponte_f_f
implicit none
integer algflag   ! input:   0=basic,  1=shift acceleration, 2=tau shift 3=double shift  4= diagonal adjustment
integer goalflag  ! input:   1= R only, 2=  R and Pi
integer errorflag ! output:  0= no error else some error did force premature termination
integer itermax   ! input:   max number of allowed iterations, default=50
integer errmax    ! input:   eps = epsilon(1.0d0)* 10 ** errmax
integer dual      ! input:   0=automatic, 1=ramaswami, 2=brigth 3= both
!
!
!
if (debug) write (*,*) "gim1_lr_solve: entering with algflag=",algflag,&
       " goalflag=",goalflag," errorflag=",errorflag," eps=",errmax," niter=",itermax," dual=",dual
!
errorflag=0
finish=.false.
info =0; fdrift=0
!
!
! deallocate G,R,U,Pi,Pi0 if allocated
!
if (allocated(G)) deallocate(G)
if (allocated(R)) deallocate(R)
if (allocated(U)) deallocate(U)
if (allocated(Pi)) deallocate(Pi)
if (allocated(Pi0)) deallocate(Pi0)

write(wout,*)"********** algorithm not ready ********** work in progress **********"
call print_it
!!
if (debug) write (*,*) "gim1_lr_solve: exiting with algflag=",algflag,&
          " goalflag=",goalflag,"errorflag=",errorflag," finish =",finish," dual=",dual

end subroutine gim1_lr_solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem GI/M/1 Algorithm Functional Iteration
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine gim1_fi_solve(algflag, goalflag,errorflag,errmax,itermax,dual)
use ponte_f_f
use smc_int
use fi_int
implicit none
integer algflag   ! input:   100=natural 200=traditional 300=U-based  
                  !             +1=Shift +2=TauShift +3=Double Shift
                  !             +30=x(0)=0 +20=x(0)=I +10=x(0)=user X
integer goalflag  ! input:   1= R only, 2= R and Pi
integer errorflag ! output:  0= no error else some error did force premature termination
integer itermax   ! input:   max number of allowed iterations, default=50
integer errmax    ! input:   eps = epsilon(1.0d0)* 10 ** errmax
integer dual      ! input:   0=automatic, 1=ramaswami, 2=brigth 3= both
real(dp)::drift,norm1,myeps,tau     !
integer m
!logical ::  ds    ! do shift
integer :: method ! 1 natural, 2 traditional, 3 U-based (default)
integer :: i,nba,ldual
real(dp),dimension(:),allocatable :: vi,v
real(dp),dimension(:,:,:),allocatable :: AA
!
!
!
nba=size(a,3)
if (debug) write (*,*) "gim1_fi_solve: entering with algflag=",algflag,&
      " goalflag=",goalflag," errorflag=",errorflag," eps=",errmax," niter=",itermax," dual=",dual
!
errorflag=0
finish=.false.
info =0; fdrift=0
myeps=epsilon(1.d0) * 10.0d0 ** errmax
!
!
! deallocate G,R,U,Pi,Pi0 if allocated
!
if (allocated(G)) deallocate(G)
if (allocated(R)) deallocate(R)
if (allocated(Pi)) deallocate(Pi)
if (allocated(Pi0)) deallocate(Pi0)
if (allocated(vi)) deallocate(vi)
if (allocated(v)) deallocate(v)
if (allocated(AA)) deallocate(AA)

m=size(A,1)

allocate(G(m,m),STAT=info)
if (debug) write (*,*) "gm1_fi_solve:  Allocating G, stat=",info
if (info /= 0 ) then
     errorflag=info
     return
endif
allocate(R(m,m),STAT=info)
if (debug) write (*,*) "gm1_fi_solve:  Allocating R, stat=",info
if (info /= 0 ) then
     errorflag=info
     return
endif


allocate(AA(m,m,nba),STAT=info)
if (debug) write (*,*) "gm1_fi_solve:  Allocating AA, stat=",info
if (info /= 0 ) then
     errorflag=info
     return
endif

allocate(v(m),STAT=info)
if (debug) write (*,*) "gm1_fi_solve:  Allocating v, stat=",info
if (info /= 0 ) then
     errorflag=info
     return
endif

method=3   

!
!setup flags
!
if(algflag >= 300) then !U-based (default)
   method=3 
   algflag=algflag-300
endif
if(algflag >= 200) then !Traditional
   method=2 
   algflag=algflag-200
endif
if(algflag >= 100) then !Natural
   method=1
   algflag=algflag-100
endif

if(algflag .GE. 30) then !  x(0)=0 (default)
      if (allocated(X0)) deallocate(X0)
      allocate(X0(m,m),STAT=info)
      if (debug) write (*,*) "gm1_fi_solve:  Allocating X0, stat=",info
      if (info /= 0 ) then
          errorflag=info
          return
      endif
      X0=0.0d0
      algflag=algflag-30
endif
if(algflag .GE. 20) then !  x(0)=I 
      if (allocated(X0)) deallocate(X0)
      allocate(X0(m,m),STAT=info)
      if (debug) write (*,*) "gm1_fi_solve:  Allocating X0, stat=",info
      if (info /= 0 ) then
          errorflag=info
          return
      endif
      X0=0.0d0
      do i=1,m
      X0(i,i)=1.0d0
      enddo
      algflag=algflag-20
endif
if (algflag .GE. 10) then !user defined x(0)
   algflag=algflag-10
   if (debug) write (*,*) "gm1_fi_solve: used defined x(0)"
endif



!
! setup goals
! only G for now
!
ldual=dual ! local dual
fdrift=-drftmg1(a)


if(ldual==0) then
    if (fdrift.ge.0) then ! G/M/1 transient
      ldual=1 !ram dual
   else ! G/M/1 recurrent
      ldual=2 !bright dual
   end if
end if

if (ldual==1) call gm1tomg1(A,AA,v)
if (ldual==2) call gm1tomg1b(A, AA,tau, v)


call fi(A=AA, ds=algflag, method=method, eps=myeps, maxit=itermax,&
 x0=X0, G=G, drift=drift, err=norm1)
       r=transpose(g)
     if (allocated(vi)) deallocate(vi)
     allocate(vi(m),stat=info)
     if (info/=0)then
        return
     end if
     vi=1.d0/v
     do i=1,m
        R(:,i)=R(:,i)*v(i)
        R(i,:)=R(i,:)*vi(i)
     end do
     if (ldual==2) then
        r=r*tau
     end if
     if(dual.eq.0)then
        if (ldual.eq.1) then
           write(wout,*)"Automatic dual: Ramaswami dual"
           call print_it
        end if
        if (ldual.eq.2) then
           write(wout,*)"Automatic dual: Bright dual"
           call print_it
        end if
     end if

     call rresidual(A,R,norm1)
     write(wout,*)"R-residual=",norm1
     call print_it
     write(wout,*)"drift=",fdrift
     call print_it
     if (debug) write (*,*) "gm1_fi_solve:  FI  did find R"
     errorflag=info
     finish=.true.

if ( .not. verb) then
       write(wout,*)" "
       call print_it
     endif
!!
if (debug) write (*,*) "gim1_fi_solve: exiting with algflag=",algflag,&
            " goalflag=",goalflag,"errorflag=",errorflag," finish =",finish," dual=",dual
end subroutine gim1_fi_solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem GI/M/1 Algorithm Invariant Subspace
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine gim1_is_solve(algflag, goalflag,errorflag,errmax,itermax,dual)
use ponte_f_f
use smc_int
use is_int

implicit none
integer algflag   ! input:   0=basic 1=Balzer 2=Shur
integer goalflag  ! input:   1= R only, 2= R and Pi
integer errorflag ! output:  0= no error else some error did force premature termination
integer itermax   ! input:   max number of allowed iterations, default=50
integer errmax    ! input:   eps = epsilon(1.0d0)* 10 ** errmax
integer dual      ! input:   0=automatic, 1=ramaswami, 2=brigth 3= both
real(dp)::drift,norm1,myeps,tau     !
integer m
integer :: nba, i, ldual
real(dp),dimension(:,:,:),allocatable:: AA
real(dp),dimension(:),allocatable:: v,vi
!
!
!
nba=size(a,3)
if (debug) write (*,*) "gim1_is_solve: entering with algflag=",algflag,&
        " goalflag=",goalflag," errorflag=",errorflag," eps=",errmax," niter=",itermax," dual=",dual
!
errorflag=0
finish=.false.
info =0; fdrift=0
myeps=epsilon(1.d0) * 10.0d0 ** errmax
!
!
! deallocate G,R,U,Pi,Pi0 if allocated
!
if (allocated(G)) deallocate(G)
if (allocated(R)) deallocate(R)
if (allocated(Pi)) deallocate(Pi)
if (allocated(Pi0)) deallocate(Pi0)

m=size(A,1)
allocate(G(m,m),STAT=info)
if (debug) write (*,*) "gm1_is_solve:  Allocating G, stat=",info
if (info /= 0 ) then
   errorflag=info
   return
endif

allocate(r(m,m),STAT=info)
if (debug) write (*,*) "gm1_is_solve:  Allocating R, stat=",info
if (info /= 0 ) then
   errorflag=info
   return
endif


allocate(AA(m,m,nba),STAT=info)
if (debug) write (*,*) "gm1_is_solve:  Allocating AA, stat=",info
if (info /= 0 ) then
     errorflag=info
     return
endif

allocate(v(m),STAT=info)
if (debug) write (*,*) "gm1_is_solve:  Allocating v, stat=",info
if (info /= 0 ) then
     errorflag=info
     return
endif

ldual=dual ! local dual
fdrift=-drftmg1(a)


if(ldual==0) then
    if (fdrift.ge.0) then ! G/M/1 transient
      ldual=1 !ram dual
   else ! G/M/1 recurrent
      ldual=2 !bright dual
   end if
end if

if (ldual==1) call gm1tomg1(A,AA,v)
if (ldual==2) call gm1tomg1b(A, AA,tau, v)

call is(A=AA, method=algflag, eps=myeps, maxit=itermax, G=G, drift=drift, err=norm1)
!!
r=transpose(g)
     if (allocated(vi)) deallocate(vi)
     allocate(vi(m),stat=info)
     if (info/=0)then
        return
     end if
     vi=1.d0/v
     do i=1,m
        R(:,i)=R(:,i)*v(i)
        R(i,:)=R(i,:)*vi(i)
     end do
     if (ldual==2) then
        r=r*tau
     end if
     if(dual.eq.0)then
        if (ldual.eq.1) then
           write(wout,*)"Automatic dual: Ramaswami dual"
           call print_it
        end if
        if (ldual.eq.2) then
           write(wout,*)"Automatic dual: Bright dual"
           call print_it
        end if
     end if

     call rresidual(A,R,norm1)
     write(wout,*)"R-residual=",norm1
     call print_it
     write(wout,*)"drift=",fdrift
     call print_it
     if (debug) write (*,*) "gm1_is_solve:  IS  did find R"
     errorflag=info
     finish=.true.

if ( .not. verb) then
       write(wout,*)" "
       call print_it
     endif

if (debug) write (*,*) "gim1_is_solve: exiting with algflag=",algflag,&
              " goalflag=",goalflag,"errorflag=",errorflag," finish =",finish," dual=",dual
end subroutine gim1_is_solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine  qbd_pi_ponte(errorflag,maxnc,epspi,pinc,boptflag,dual)
use ponte_f_f
use pi_int
!        globals in ponte_f_f.mod
! real(KIND=8),dimension(:,:,:),allocatable:: A,B
! real(KIND=8),dimension(:,:),allocatable:: G,R,U,Pi
! real(KIND=8),dimension(:,:),allocatable:: Pi0
! logical debug  !if true debugging output to std out
! logical verb   !if true more detailed output to NotePad
! character (len=256) wout  !line of output
!
! integer,parameter :: dp=kind(0.d0)
! integer :: info=0 info /=0 will catch allocation errors
!
!        globals in pi_int   
! 
implicit none
integer errorflag ! output:  0= no error else some error did force premature termination
integer :: maxnc, pinc  ! max and computed number of block components of Pi
integer :: boptflag ! 1 if b is present  0 if not present
integer :: epspi    ! input:   myeps = epsilon(1.0d0)* 10 ** epspi
integer dual      ! input:   0=automatic, 1=ramaswami, 2=brigth 3= both
real(dp):: myeps  
integer m,n
!integer :: m0     ! size of the block B0
!
errorflag=0
finish=.false.
info =0
myeps=epsilon(1.d0) * 10.0d0 ** epspi
errorflag=info ! what is this ??
!
!
if (debug) write (*,*) "qbd_pi_ponte: entering with errorflag=",&
errorflag,"maxnc=",maxnc,"epspi=",epspi,"myeps=",myeps,"bn1flag=",boptflag
!
!

if (allocated(Pi)) deallocate(Pi)
if (allocated(Pi0)) deallocate(Pi0)

m=size(a,1)
n=size(b0,1)
allocate(pi(m,maxnc),STAT=info)
if (debug) write (*,*) "qbd_pi_ponte:  Allocating pi, stat=",info
if (info /= 0 ) then
   errorflag=info
   return
endif

allocate(pi0(n),STAT=info)
if (debug) write (*,*) "qbd_pi_ponte:  Allocating pi0, stat=",info
if (info /= 0 ) then
   errorflag=info
   return
endif
if (boptflag==1) then
call qbdpi(am1=A(:,:,1),a0=A(:,:,2),a1=A(:,:,3), bm1=BN1, b0=B0, b1=B(:,:,1),  R=R,&
 maxnc=maxnc, epspi=myeps, pi0=pi0, pi=pi, pinc=pinc)
else
call qbdpi(am1=A(:,:,1),a0=A(:,:,2),a1=A(:,:,3),bm1=BN1,  b0=B0,  R=R,&
 maxnc=maxnc, epspi=myeps, pi0=pi0, pi=pi, pinc=pinc)
end if
if (verb) write (wout,*) " Number of computed blocks of components of pi= ",pinc
call print_it
end subroutine qbd_pi_ponte

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem MG1  Pi computations
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine  mg1_pi_ponte(errorflag,maxnc,epspi,pinc,bn1flag)
use ponte_f_f
use pi_int

implicit none
integer errorflag ! output:  0= no error else some error did force premature termination
integer :: maxnc, pinc  ! max and computed number of block components of Pi
integer :: epspi    ! input:   myeps = epsilon(1.0d0)* 10 ** epspi
integer :: bn1flag ! 1 if bn1 is present  0 if not present 
real(dp):: myeps   
integer m,n
!integer :: m0     ! size of the block B0
!
finish=.false.
info =0
myeps=epsilon(1.d0) * 10.0d0 ** epspi
errorflag=info
!
!
if (debug) write (*,*) "mg1_pi_ponte: entering with errorflag=",errorflag,"maxnc=",&
maxnc,"epspi=",epspi,"myeps=",myeps,"bn1flag",bn1flag
!
!

if (allocated(Pi)) deallocate(Pi)
if (allocated(Pi0)) deallocate(Pi0)


m=size(a,1)
n=size(b0,1)

allocate(pi(m,maxnc),STAT=info)
if (debug) write (*,*) "mg1_pi_ponte:  Allocating pi, stat=",info
if (info /= 0 ) then
   errorflag=info
   return
endif

allocate(pi0(n),STAT=info)
if (debug) write (*,*) "mg1_pi_ponte:  Allocating pi0, stat=",info
if (info /= 0 ) then
   errorflag=info
   return
endif

if (bn1flag==1) then
call mg1pi(A=A, B0=B0, B=B,  Bm1=Bn1, G=G, maxnc=maxnc, epspi=myeps, pi0=pi0, pi=pi, pinc=pinc)
else
call mg1pi(A=A, B0=B0, B=B,  G=G, maxnc=maxnc, epspi=myeps, pi0=pi0, pi=pi, pinc=pinc)
end if
if (verb)  write (wout,*) " Number of computed blocks of components of pi= ",pinc
end subroutine mg1_pi_ponte

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem GIM1  Pi computations
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine  gim1_pi_ponte(errorflag,maxnc,epspi,pinc,bn1flag)
use ponte_f_f
use pi_int

implicit none
integer errorflag ! output:  0= no error else some error did force premature termination
integer :: maxnc, pinc  ! max and computed number of block components of Pi
integer :: epspi    ! input:   myeps = epsilon(1.0d0)* 10 ** epspi
integer :: bn1flag ! 1 if bn1 is present  0 if not present 
real(dp):: myeps  
integer m,n

finish=.false.
info =0
myeps=epsilon(1.d0) * 10.0d0 ** epspi
errorflag=info ! what is this ???
!
!
if (debug) write (*,*) "gim1_pi_ponte: entering with errorflag=",errorflag,"maxnc=",&
maxnc,"epspi=",epspi,"myeps=",myeps,"bn1flag",bn1flag
!
!

if (allocated(Pi)) deallocate(Pi)
if (allocated(Pi0)) deallocate(Pi0)


m=size(a,1)
n=size(b0,1)

allocate(pi(m,maxnc),STAT=info)
if (debug) write (*,*) "gim1_pi_ponte:  Allocating pi, stat=",info
if (info /= 0 ) then
   errorflag=info
   return
endif

allocate(pi0(n),STAT=info)
if (debug) write (*,*) "gim1_pi_ponte:  Allocating pi0, stat=",info
if (info /= 0 ) then
   errorflag=info
   return
endif

if (bn1flag==1) then
call gm1pi(A=A, B=B, B0=B0, Bm1=BN1, R=R, maxnc=maxnc, epspi=myeps, pi0=pi0, pi=pi, pinc=pinc)
else
call gm1pi(A=A, B=B, B0=B0,  R=R, maxnc=maxnc, epspi=myeps, pi0=pi0, pi=pi, pinc=pinc)
end if
if (verb)  write (wout,*) " Number of computed blocks of components of pi= ",pinc
end subroutine gim1_pi_ponte
