!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  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                                   !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
!
!module smc_tools
  ! Default RC simple shift to 0     if drift<0 
  !                         to infty if drift>0 , 
  ! option no shift

  ! main_subroutine: se compare solo G calcola solo G
  !                  G e R calcola G e R
  !                  U
  !                  pi
  !       output ulteriore: drift, residuo

  ! raffinamento con Newton se residuo e' grande (opzionale?) 
  ! (da valutare con esperimenti)
!  integer,parameter :: dp=kind(0.d0)
!  integer :: info=0
!end module smc_tools

module smc_int
interface
function drftmg1(a)
integer,parameter::dp=kind(0.d0)
  real(dp),dimension(:,:,:)::a
  real(dp)::drftmg1
end function drftmg1

end interface
  interface 
     subroutine drft(am1, a0, a1, v, drift)
       integer, parameter :: dp=kind(0.d0)
       real(dp),intent(in),dimension(:,:) :: am1, a0, a1
       real(dp), intent(out)              :: drift
       real(dp),dimension(:),intent(out)  :: v
     end subroutine drft
  end interface
  interface 
     subroutine gth(a,p)
       integer, parameter :: dp=kind(0.d0)
       real(dp),dimension(:,:),intent(in) :: a
       real(dp),dimension(:),intent(out)  :: p
     end subroutine gth
  end interface
  interface 
     subroutine crqbd(Am1, A0, A1, G, R, U, doshift, dogth, drift, nerror, maxit)
       integer, parameter :: dp=kind(0.d0)
       real(dp),dimension(:,:),intent(in)  :: Am1, A0, A1
       real(dp),dimension(:,:),intent(out) :: G
       real(dp),dimension(:,:),intent(out),optional ::  R, U
       real(dp),intent(out),optional ::drift
       logical, optional,intent(in) :: dogth
       integer, optional,intent(in) :: doshift
       integer,optional,intent(in) :: nerror
       integer, optional,intent(in) :: maxit
     end subroutine crqbd
  end interface
  interface 
     subroutine lrqbd(Am1, A0, A1,  G, R, U, doshift, dogth, drift, nerror, maxit)
       integer, parameter :: dp=kind(0.d0)
       real(dp),dimension(:,:),intent(in) :: Am1, A0, A1
       real(dp),dimension(:,:),intent(out),optional :: G, R, U
       real(dp),intent(out),optional :: drift
       logical, optional,intent(in) :: dogth
       integer, optional,intent(in) :: doshift
       integer,optional,intent(in) :: nerror
       integer, optional,intent(in) :: maxit
     end subroutine lrqbd
  end interface
  interface
     subroutine bgth(v,a,b)
       integer, parameter :: dp=kind(0.d0)
       real(dp),intent(in), dimension(:,:) :: a
       real(dp),intent(inout), dimension(:,:) ::  b
       real(dp),intent(inout),dimension(:)   :: v
     end subroutine bgth
  end interface
  interface
     SUBROUTINE shift(a,ds,tau,uu,vv,drift,ta)
       integer, parameter :: dp=kind(0.d0)
       REAL(dp), DIMENSION(:,:,:):: a, ta
       REAL(dp), DIMENSION(:) :: uu,vv
       INTEGER :: ds
       REAL(dp)                                :: drift,tau
     end SUBROUTINE shift
  end interface
  
interface
subroutine gm1tomg1(A, AA, v)
  integer,parameter::dp=kind(0.d0)
  real(dp),dimension(:,:,:):: A,AA
  real(dp),dimension(:):: v
end subroutine gm1tomg1
end interface

interface
   subroutine gm1tomg1b(A, AA,tau, w)
     integer,parameter::dp=kind(0.d0)
     real(dp),dimension(:,:,:):: A,AA
     real(dp),dimension(:):: w
     real(dp):: tau 
   end subroutine gm1tomg1b
end interface

interface
   subroutine mg1tomg1(A, AA,tau, w)
     integer,parameter::dp=kind(0.d0)
     real(dp),dimension(:,:,:):: A,AA
     real(dp),dimension(:):: w
     real(dp):: tau 
   end subroutine mg1tomg1
end interface

interface
   subroutine rresidual(A,r,norm1)
     integer,parameter::dp=kind(0.d0)
     real(dp),dimension(:,:,:):: A
     real(dp),dimension(:,:):: R
     real(dp)::norm1
   end subroutine rresidual
end interface

interface
     subroutine gtou(A0,A1,G,U)
       integer,parameter:: dp=kind(0.d0)
       real(dp),dimension(:,:) :: a0,a1,g,u
     end subroutine gtou
end interface

interface
     subroutine utor(A1,u,r)
       integer,parameter:: dp=kind(0.d0)
       real(dp),dimension(:,:) :: a1,u,r
     end subroutine utor
end interface

interface
subroutine qbdrres(Am1,A0,A1,r,res) 
  integer,parameter:: dp=kind(0.d0)
  real(dp),dimension(:,:) :: am1,a0,a1,r
  real(dp)::res
end subroutine qbdrres
end interface

interface
subroutine qbdures(Am1,A0,A1,u,res) 
  integer,parameter:: dp=kind(0.d0)
  real(dp),dimension(:,:) :: am1,a0,a1,u
  real(dp)::res
end subroutine qbdures
end interface

interface
subroutine compute_tau(a,z,u,v)
! Computes the solution different from 1 of the equation
! rho(S(z))-z=0, by means of the deflated Newton method
! z=z-nwt/(1-nwt/(z-1)), nwt=(rho(S(z))-z)/(v^TS(z)'u-1)
! where S(z)u=rho u, v^TS(z)=rho v^T
! S(z)=sum_i A_i z^i
  implicit none
  integer,parameter::dp=kind(0.d0)
  integer :: n,m
  real(dp),intent(in),dimension(:,:,:)::a
  real(dp),intent(out) :: z
  real(dp),dimension(:)::u,v

end subroutine compute_tau
end interface

interface
subroutine compute_tauqbd(am1,a0,a1,tau)
  integer,parameter:: dp=kind(0.d0)
    real(dp),intent(in),dimension(:,:) :: am1,a0,a1
!    real(dp),intent(in)               :: drift
    real(dp),intent(out)               :: tau
  end subroutine compute_tauqbd
end interface

interface
subroutine taushift(a, tau, u, v, drift, b)
  integer,parameter:: dp=kind(0.d0)
  real(dp),dimension(:,:,:) :: a,b
  real(dp),dimension(:) :: u,v
 real(dp)::tau,drift
end subroutine taushift
end interface

interface
subroutine mg_startup(a,ds,dual,tau,u,v,drft)
  integer,parameter:: dp=kind(0.d0)
  real(dp),dimension(:,:,:) :: a
  real(dp) :: tau, drft
  real(dp),dimension(:) :: u,v
  integer :: ds, dual
end subroutine mg_startup
end interface

end module smc_int


