!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  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 1.3 - March 2007                                   !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   fcalls.f90   main calls from C to  Fortran computation    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem QBD Algorithm Cyclic Reduction
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine qbd_cr_solve(algflag, goalflag,errorflag,errmax,itermax)
use ponte_f_f
use smc_tools
use smc_int
! 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
!
!        prototypes in smc_int
! subroutine drft(am1, a0, a1, v, drift)
! subroutine gth(a,p)
! subroutine crqbd(Am1, A0, A1, G, R, U, doshift, dogth, verbose,  drift, debugging)
! subroutine lrqbd(Am1, A0, A1,  G, R, U, doshift, dogth, verbose, drift)
! subroutine bgth(v,a,b)
!
!        globals in smc_tools
! integer,parameter :: dp=kind(0.d0)
! integer :: info=0 info /=0 will catch allocation errors
! 
implicit none
integer algflag   ! input:   0=basic,  1=shift acceleration,  2= 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 ::  ds, dg  ! algorithm flags
real(dpr)::drift 
integer m         ! for allocating matrices
! integer :: maxnc  ! max number of block components of Pi
!
!
!
if (debug) write (*,*) "qbd_cr_solve: entering with algflag=",algflag, &
          " goalflag=",goalflag," errorflag=",errorflag,"eps=",errmax,"niter=",itermax
!
finish=.false.
info =0; fdrift=0
!
!
! disallocate 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
   ds=.false.
   dg=.false.
case (1) !shift
   ds=.true.
   dg=.false. 
case (2) !diagonal adjustment
   ds=.false.
   dg=.true.
case default
   ds=.false.
   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=info 
     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=info
     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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!PI PI PI
!     allocate(pi0(m0),STAT=info)
!     if (debug) write (*,*) "qbd_cr_solve:  Allocating pi0, stat=",info
!     if (info /= 0 ) then
!     errorflag=info
!     return
!     endif
!     allocate(pi(m,maxnc),STAT=info)
!     if (debug) write (*,*) "qbd_cr_solve:  Allocating pi, stat=",info
!     if (info /= 0 ) then
!     errorflag=info
!     return
!     endif
!     Bm1=Am1; B0=A0+Am1
!     m0=m
!     call qbd_pi(am1=A(:,:,1),a0=A(:,:,2),a1=A(:,:,3), bm1=A(:,:,1), b0=A(:,:,1)+A(:,:,2), b1=A(:,:,3),  R=R, maxnc=maxnc, epspi=myeps, pi0=pi0, pi=pi, pinc=maxnc)
!     call qbd_pi(Am1, A0, A1, Bm1, B0, B1=B1, R, maxnc, epspi, pi0, pi)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!PI PI PI

!     if(debug) write(*,*) "qbd_cr_solve: computation of G R U and Pi not yet complete ..."
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_tools
use smc_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
!
!        prototypes in smc_int
! subroutine drft(am1, a0, a1, v, drift)
! subroutine gth(a,p)
! subroutine crqbd(Am1, A0, A1, G, R, U, doshift, dogth, verbose,  drift, debugging)
! subroutine lrqbd(Am1, A0, A1,  G, R, U, doshift, dogth, verbose, drift)
! subroutine bgth(v,a,b)
!
!        globals in smc_tools
! integer,parameter :: dp=kind(0.d0)
! integer :: info=0 info /=0 will catch allocation errors
! 
implicit none
integer algflag   ! input:   0=basic,  1=shift acceleration,  2= 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 ::  ds, dg  ! 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
!
finish=.false.
info =0; fdrift=0
!
! disallocate 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
   ds=.false.
   dg=.false.
case (1) !shift
   ds=.true.
   dg=.false.
case (2) !diagonal adjustment
   ds=.false.
   dg=.true.
case default
   ds=.false.
   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
!
! added by sergio
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 smc_tools
use fi_int
use smc_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
!
!        globals in smc_tools
! integer,parameter :: dp=kind(0.d0)
! integer :: info=0 info /=0 will catch allocation errors
!
!        globals inpwcr_interface
! subroutine pwcr(a,eps,g,err)  
! 
implicit none
integer algflag   ! input:   0=basic x(0)=0,  1=basic x(0)=1, 2=shift acceleration
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
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
!
finish=.false.
info =0; fdrift=0
myeps=epsilon(1.d0) * 10.0d0 ** errmax
!
!
! disallocate 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)

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

ds=.false.
method=3   
!
! setup flags
! algflag=0,1,2,3,4,5,6,7,8
!
select case(algflag)
case (0) !Natural, Basic, x(0)=0 
      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
      method = 1
case (1) !Natural, Basic, x(0)=1
      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
      method = 1
      do i=1,m
      X0(i,i)=1.0d0
      enddo
case (2) !Natural, Shift Acceleration
      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
      method = 1
      ds=.true.
case (3) !Traditional, Basic, x(0)=0
      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
      method = 2
case (4) !Traditional, Basic, x(0)=1
      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
      method = 2
      do i=1,m
      X0(i,i)=1.0d0
      enddo
case (5) !Traditional, Shift Acceleration
      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
      method = 2
      ds=.true.
case (6) !U-based, Basic, x(0)=0
      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
case (7) !U-based, Basic, x(0)=1
      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
case (8) !U-based, Shift Acceleration
      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
      ds=.true.
case(9) !Natural, Basic, x(0) user
      method = 1
case(10) !Traditional, Basic, x(0) user
      method = 2
case(11) !U-based, Basic, x(0) user
case(12) !Natural, Shift Acceleration, x(0) user
      method = 1
      ds=.true.
case(13) !Traditional, Shift Acceleration, x(0) user
      method = 2
      ds=.true.
case(14) !U-based, Shift Acceleration, x(0) user
      ds=.true.
case default
if(debug) write(*,*) "qbd_fi_solve: strange algflag=",algflag
end select
!
! setup goals
!
!
! setup goals
!
call fi(A=A, doshift=ds, 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=",drift
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 smc_tools
use is_int
use smc_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
!
!        globals in smc_tools
! integer,parameter :: dp=kind(0.d0)
! integer :: info=0 info /=0 will catch allocation errors
!
!        globals in pwcr_interface
! subroutine pwcr(a,eps,g,err)  
! 
implicit none
integer algflag   ! input:   0=basic, ignore it !
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  
!
!
!
if (debug) write (*,*) "qbd_is_solve: entering with algflag=",algflag,&
         " goalflag=",goalflag," errorflag=",errorflag," eps=",errmax," niter=",itermax
!
finish=.false.
info =0; fdrift=0
myeps=epsilon(1.d0) * 10.0d0 ** errmax
!
!
! disallocate 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)

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

!!
!call is(A, method=algflag, eps=myeps, maxit=itermax, G=G, drift=drift, err=norm1);
!!write(wout,*)"********** algorithm not ready ********** work in progress **********"
!!call print_it
!!

!!
call is(A, method=algflag, eps=myeps, maxit=itermax, G=G, drift=drift, err=norm1);
!!write(wout,*)"********** algorithm not ready ********** work in progress **********"
!!call print_it
!!
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=",drift
 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)
use ponte_f_f
use smc_tools
use pwcr_interface
! 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
!
!        globals in smc_tools
! integer,parameter :: dp=kind(0.d0)
! integer :: info=0 info /=0 will catch allocation errors
!
!        globals in pwcr_interface
! subroutine pwcr(a,eps,g,err)  
! 
implicit none
integer algflag   ! input:   0=basic,  1=shift acceleration,  2= 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
real(dp)::norm1,myeps     
integer m
!
!
!

if (debug) write (*,*) "mg1_cr_solve: entering with algflag=",algflag,&
    " goalflag=",goalflag," errorflag=",errorflag," eps=",errmax," niter=",&
      itermax," interpmax=",interpmax
!
finish=.false.
info =0; fdrift=0
myeps=epsilon(1.d0) * 10.0d0 ** errmax
!
!
! disallocate 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)

m=size(A,1)

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=A,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==1) call  spwcr(A=A,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax) 
     if (algflag==2) call  pwcr(A=A,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax) 
     if (info /= 0 ) then
       errorflag=info
        if (info == 500) 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
     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=A,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax) 
     if (algflag==1) call  spwcr(A=A,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==2) call  pwcr(A=A,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (info /= 0 ) then
       errorflag=info
        if (info == 500) 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
     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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!PI PI PI
!     na=size(a,3)
!     m0=size(a,1)!!!!pasticcio
!     m=size(a,1)
!     maxnc=900
!     allocate(pi0(m0),STAT=info)
!     if (debug) write (*,*) "qbd_cr_solve:  Allocating pi0, stat=",info
!     if (info /= 0 ) then
!     errorflag=info
!     return
!     endif
!     allocate(pi(m,maxnc),STAT=info)
!     if (debug) write (*,*) "qbd_cr_solve:  Allocating pi, stat=",info
!     if (info /= 0 ) then
!     errorflag=info
!     return
!     endif
!     call mg1_pi(a=a, b0=a(:,:,2), B=A(:,:,2:na), G=G, maxnc=maxnc,epspi=1.d-14,pi0=pi0,pi=pi)!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!PI PI PI

     finish=.true.
!     if (debug) write (*,*) "mg1_cr_solve:  G and Pi not yet implemented"
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
end subroutine mg1_cr_solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem M/G/1 Algorithm Functional Iteration
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine mg1_fi_solve(algflag, goalflag,errorflag,errmax,itermax)
use ponte_f_f
use smc_tools
use fi_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
!
!        globals in smc_tools
! integer,parameter :: dp=kind(0.d0)
! integer :: info=0 info /=0 will catch allocation errors
!
!        globals inpwcr_interface
! subroutine pwcr(a,eps,g,err)  
! 
implicit none
integer algflag   ! input:   0=basic x(0)=0, 1=basic x(0)=1,  1=shift acceleration
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
real(dp)::drift,norm1,myeps     !
integer m
logical ::  ds    ! do shift
integer :: method ! 1 natural, 2 traditional, 3 U-based (default)
integer :: i
!
!
!
if (debug) write (*,*) "mg1_fi_solve: entering with algflag=",algflag,&
           " goalflag=",goalflag," errorflag=",errorflag," eps=",errmax,&
           " niter=",itermax
!
finish=.false.
info =0; fdrift=0
myeps=epsilon(1.d0) * 10.0d0 ** errmax
!
!
! disallocate 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)
!
m=size(A,1)
!
!
allocate(G(m,m),STAT=info)
if (debug) write (*,*) "mg1_fi_solve:  Allocating G, stat=",info
if (info /= 0 ) then
     errorflag=info
     return
endif

ds=.false.
method=3   
!
! setup flags
! algflag=0,1,2,3,4,5,6,7,8
!
select case(algflag)
case (0) !Natural, Basic, x(0)=0 
      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
      method = 1
case (1) !Natural, Basic, x(0)=1
      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
      method = 1
      do i=1,m 
      X0(i,i)=1.0d0
      enddo
case (2) !Natural, Shift Acceleration
      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
      method = 1
      ds=.true.
case (3) !Traditional, Basic, x(0)=0
      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
      method = 2
case (4) !Traditional, Basic, x(0)=1
      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
      method = 2
      do i=1,m
      X0(i,i)=1.0d0
      enddo
case (5) !Traditional, Shift Acceleration
      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
      method = 2
      ds=.true.
case (6) !U-based, Basic, x(0)=0
      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
case (7) !U-based, Basic, x(0)=1
      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
case (8) !U-based, Shift Acceleration
      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
      ds=.true.
case(9) !Natural, Basic, x(0) user
      method = 1
case(10) !Traditional, Basic, x(0) user
      method = 2
case(11) !U-based, Basic, x(0) user
case(12) !Natural, Shift Acceleration, x(0) user
      method = 1
      ds=.true.
case(13) !Traditional, Shift Acceleration, x(0) user
      method = 2
      ds=.true.
case(14) !U-based, Shift Acceleration, x(0) user
      ds=.true.
case default
if(debug) write(*,*) "mg1_fi_solve: strange algflag=",algflag
end select
!
! setup goals
! only G for now
!
call fi(A=A, doshift=ds, 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
write(wout,*)"drift=",drift
 call print_it

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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem M/G/1 Algorithm Invariant Subspace
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine mg1_is_solve(algflag, goalflag,errorflag,errmax,itermax)
use ponte_f_f
use smc_tools
use is_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
!
!        globals in smc_tools
! integer,parameter :: dp=kind(0.d0)
! integer :: info=0 info /=0 will catch allocation errors
!
!        globals inpwcr_interface
! subroutine pwcr(a,eps,g,err)  
! 
implicit none
integer algflag   ! input:   0=basic x(0)=0, 1=basic x(0)=1,  2=shift acceleration
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
real(dp)::drift,norm1,myeps     !
integer m
!
!
!
if (debug) write (*,*) "mg1_is_solve: entering with algflag=",algflag,&
        " goalflag=",goalflag," errorflag=",errorflag," eps=",errmax," niter=",itermax
!
finish=.false.
info =0; fdrift=0
myeps=epsilon(1.d0) * 10.0d0 ** errmax
!
!
! disallocate 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)

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

call is(A, method=algflag, eps=myeps, maxit=itermax, G=G, drift=drift, err=norm1)
write(wout,*)"G-residual=",norm1
call print_it
write(wout,*)"drift=",drift
 call print_it

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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem GI/M/1 Algorithm Cyclic Reduction
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine gim1_cr_solve(algflag, goalflag,errorflag,errmax,itermax,interpmax)
use ponte_f_f
use smc_tools
use pwcr_interface
use smc_int
! 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
!
!        globals in smc_tools
! integer,parameter :: dp=kind(0.d0)
! integer :: info=0 info /=0 will catch allocation errors
!
!        globals inpwcr_interface
! subroutine pwcr(a,eps,g,err)  
! 
implicit none
integer algflag   ! input:   0=basic,  1=shift acceleration,  2= 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
real(dp)::norm1,myeps     
integer m
integer:: nba,  i
 real(dp),dimension(:),allocatable :: vi,v 
! real(dp),dimension(:,:),allocatable :: B0,bm1
 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
!
finish=.false.
info =0; fdrift=0
myeps=epsilon(1.d0) * 10.0d0 ** errmax
!
!
! disallocate 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

!!
!write(wout,*)"********** algorithm not ready ********** work in progress **********"
!call print_it
!!

call gm1tomg1(A,AA,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,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax) 
     if (algflag==2) call  pwcr(A=AA,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax) 
     if (info /= 0 ) then
       errorflag=info
        if (info == 500) 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
     call rresidual(A,R,norm1)
     write(wout,*)"R-residual=",norm1
     call print_it
     fdrift=-fdrift
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,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (algflag==2) call  pwcr(A=AA,eps=myeps,G=G,err=norm1,maxit=itermax,intpmax=interpmax)
     if (info /= 0 ) then
       errorflag=info
        if (info == 500) 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
     call rresidual(A,r,norm1)
     write(wout,*)"R-residual=",norm1
     call print_it
     fdrift=-fdrift
write(wout,*)"drift=",fdrift
 call print_it

     if (debug) write (*,*) "gm1_cr_solve:  Cyclic Reduction  did find G"
     errorflag=info

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!PI PI PI
!!$     na=size(a,3)
!!$     m0=size(a,1)!!!!pasticcio
!!$     m=size(a,1)
!!$     maxnc=900
!!$     allocate(pi0(m0),STAT=info)
!!$     if (debug) write (*,*) "qbd_cr_solve:  Allocating pi0, stat=",info
!!$     if (info /= 0 ) then
!!$     errorflag=info
!!$     return
!!$     endif
!!$     allocate(pi(m,maxnc),STAT=info)
!!$     if (debug) write (*,*) "qbd_cr_solve:  Allocating pi, stat=",info
!!$     if (info /= 0 ) then
!!$     errorflag=info
!!$     return
!!$     endif
!!$     call mg1_pi(a=a, b0=a(:,:,2), B=A(:,:,2:na), G=G, maxnc=maxnc,epspi=1.d-14,pi0=pi0,pi=pi)!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!PI PI PI

! pasticcio per definire B
!if (allocated(B)) deallocate(B)
!if (allocated(B0)) deallocate(B0)
!if (allocated(Bm1)) deallocate(Bm1)
!allocate(b(m,m,nba-2),B0(m,m),bm1(m,m),stat=info)
!     allocate(pi0(m),STAT=info)
!     if (info /= 0 ) then
!     errorflag=info
!     return
!     endif
!     maxnc=100
!     allocate(pi(m,maxnc),STAT=info)
!     if (debug) write (*,*) "qbd_cr_solve:  Allocating pi, stat=",info
!
!     if(info/=0)return
!     b=0.d0
!     B0=0
!     nbb=size(B,3)
!     v=1.d0-sum(a(:,:,1),dim=2)
!     b0(:,1)=v
!     do i=1,nbb
!        v=v-sum(a(:,:,i+1),dim=2)
!        b(:,1,i)=v
!     end do
!
!     call  gm1_pi(A=a, B=b, B0=b0, Bm1=A(:,:,1), R=r, maxnc=maxnc, epspi=1.d-14, pi0=pi0, pi=pi)
!     call  gm1_pi(A=a, B=b, B0=b0, R=r, maxnc=maxnc, epspi=1.d-14, pi0=pi0, pi=pi)
     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
end subroutine gim1_cr_solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem GI/M/1 Algorithm Logarithmic Reduction
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine gim1_lr_solve(algflag, goalflag,errorflag,errmax,itermax)
use ponte_f_f
use smc_tools
!
!        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
!
!        globals in smc_tools
! integer,parameter :: dp=kind(0.d0)
! integer :: info=0 info /=0 will catch allocation errors
!
!        globals inpwcr_interface
! subroutine pwcr(a,eps,g,err)  
! 
implicit none
integer algflag   ! input:   0=basic,  1=shift acceleration,  2= 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
!
!
!
if (debug) write (*,*) "gim1_lr_solve: entering with algflag=",algflag,&
       " goalflag=",goalflag," errorflag=",errorflag," eps=",errmax," niter=",itermax
!
finish=.false.
info =0; fdrift=0
!
!
! disallocate 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

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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem GI/M/1 Algorithm Functional Iteration
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine gim1_fi_solve(algflag, goalflag,errorflag,errmax,itermax)
use ponte_f_f
use smc_tools
use smc_int
use fi_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
!
!        globals in smc_tools
! integer,parameter :: dp=kind(0.d0)
! integer :: info=0 info /=0 will catch allocation errors
!
!        globals inpwcr_interface
! subroutine pwcr(a,eps,g,err)  
! 
implicit none
integer algflag   ! input:   0=basic x(0)=0, 1=basic x(0)=1,  2=shift acceleration
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
real(dp)::drift,norm1,myeps     !
integer m
logical ::  ds    ! do shift
integer :: method ! 1 natural, 2 traditional, 3 U-based (default)
integer :: i,nba
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
!
finish=.false.
info =0; fdrift=0
myeps=epsilon(1.d0) * 10.0d0 ** errmax
!
!
! disallocate 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)
!!
!write(wout,*)"********** algorithm not ready ********** work in progress **********"
!call print_it
!!
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

X0=0.0d0
ds=.false.
method=3   
! setup flags
! algflag=0,1,2,3,4,5,6,7,8
!
select case(algflag)
case (0) !Natural, Basic, x(0)=0 
      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
      method = 1
case (1) !Natural, Basic, x(0)=1
     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
      method = 1
      do i=1,m 
      X0(i,i)=1.0d0
      enddo
case (2) !Natural, Shift Acceleration
     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
      method = 1
      ds=.true.
case (3) !Traditional, Basic, x(0)=0
     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
      method = 2
case (4) !Traditional, Basic, x(0)=1
     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
      method = 2
      do i=1,m
      X0(i,i)=1.0d0
      enddo
case (5) !Traditional, Shift Acceleration
     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
      method = 2
      ds=.true.
case (6) !U-based, Basic, x(0)=0
     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
case (7) !U-based, Basic, x(0)=1
     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
case (8) !U-based, Shift Acceleration
     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
      ds=.true.
case(9) !Natural, Basic, x(0) user
      method = 1
case(10) !Traditional, Basic, x(0) user
      method = 2
case(11) !U-based, Basic, x(0) user
case(12) !Natural, Shift Acceleration, x(0) user
      method = 1
      ds=.true.
case(13) !Traditional, Shift Acceleration, x(0) user
      method = 2
      ds=.true.
case(14) !U-based, Shift Acceleration, x(0) user
      ds=.true.
case default
if(debug) write(*,*) "gm1_fi_solve: strange algflag=",algflag
end select
!
! setup goals
! only G for now
!

call gm1tomg1(A,AA,v)


call fi(A=AA, doshift=ds, 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
     call rresidual(A,R,norm1)
     write(wout,*)"R-residual=",norm1
     call print_it
     fdrift=-drift
     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
end subroutine gim1_fi_solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Problem GI/M/1 Algorithm Invariant Subspace
! 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine gim1_is_solve(algflag, goalflag,errorflag,errmax,itermax)
use ponte_f_f
use smc_tools
use smc_int
use is_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
!
!        globals in smc_tools
! integer,parameter :: dp=kind(0.d0)
! integer :: info=0 info /=0 will catch allocation errors
!
!        globals inpwcr_interface
! subroutine pwcr(a,eps,g,err)  
! 
implicit none
integer algflag   ! input:   0=basic ignore !
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
real(dpr)::drift,norm1,myeps     !
integer m
integer :: nba, i
real(dpr),dimension(:,:,:),allocatable:: AA
real(dpr),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
!
finish=.false.
info =0; fdrift=0
myeps=epsilon(1.d0) * 10.0d0 ** errmax
!
!
! disallocate 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

call gm1tomg1(A,AA,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
     call rresidual(A,R,norm1)
     write(wout,*)"R-residual=",norm1
     call print_it
     fdrift=-drift
     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
!!
!! write(wout,*)"********** algorithm not ready ********** work in progress **********"
!! call print_it
!!
if (debug) write (*,*) "gim1_is_solve: exiting with algflag=",algflag,&
              " goalflag=",goalflag,"errorflag=",errorflag," finish =",finish
end subroutine gim1_is_solve
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine  qbd_pi_ponte(errorflag,maxnc,epspi,pinc,boptflag)
use ponte_f_f
use smc_tools
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
!
!        globals in smc_tools
! 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
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 (*,*) "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 smc_tools
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
!
!        globals in smc_tools
! 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 :: 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
!mg1pi(A, B0, B, Bm1, G, maxnc, epspi, pi0, pi,pinc)
! this flag is unused here ?
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 smc_tools
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
!
!        globals in smc_tools
! 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 :: 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 (*,*) "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
!mg1pi(A, B0, B, Bm1, G, maxnc, epspi, pi0, pi,pinc)
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