!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  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.2 - Nov 2006                                     !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   futils.f90   utility fortran functions to be called by C  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
! drift value needed to know if Pi can be computed
subroutine getfdrift(x)
use ponte_f_f
implicit none
real(dpr)  x
x=fdrift
end subroutine getfdrift
!
!
!
! as Most Fortrans uses a private structure to hold 
!            matrices, allocation must be made in 
!            Fortran code and not in C.
!
!Am1=A(:,:,1), A0=A(:,:,1), A1=A(:,:,2), ....
! A blocks are always square
!
subroutine alloc_matA(m,qa,stat)
use ponte_f_f
implicit none
integer m,qa,stat
if (allocated(A)) deallocate(A)
allocate(A(m,m,qa),STAT=stat)
A=0.0d0
end subroutine alloc_matA

!
! to read the content of A from C routines
! remember qa=1,2,3,..
!
subroutine read_matA(n,m,qa,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m,qa
x=A(n,m,qa)
end subroutine read_matA

!
! to write to A from C routines
! remember qa=1,2,3,..
!
subroutine write_matA(n,m,qa,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m,qa
A(n,m,qa)=x
end subroutine write_matA

!
! test program
! dumps A to Main Window NotePad
!
subroutine dump_matA()
use ponte_f_f
implicit none
integer nd,md,qad,n,m,qa
! reads dimension kdmid,jd of matx
nd=size(A,1)
md=size(A,2)
qad=size(A,3)

! writes all the matrix to notePad
do qa=1,qad
 do m=1,md
   do n=1,nd
   write(wout,*)"A(",n,m,qa,")=",A(n,m,qa)
   call print_it_
   end do
 end do
end do
end subroutine dump_matA

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!The "B Matrices" are:
!   B0,  BN1 , B1=B(:,:,1), B2=B(:,:,2) ...
!depending from the problem, some are optional
!   B0 is always required
!note  block index form 1 to qb 
!   B blocks can be rectangular
subroutine alloc_matB(n,m,qb,stat)
use ponte_f_f
implicit none
integer n,m,qb,stat
if (allocated(B)) deallocate(B)
allocate(B(n,m,qb),STAT=stat)
B=0.0d0
end subroutine alloc_matB

!
! to read the content of B from C routines
!
!
subroutine read_matB(n,m,qb,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m,qb
x=B(n,m,qb)
end subroutine read_matB

!
! to write to B from C routines
!
!
subroutine write_matB(n,m,qb,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m,qb
B(n,m,qb)=x
end subroutine write_matB

!
! test program
! dumps B to Main Window NotePad
!
subroutine dump_matB()
use ponte_f_f
implicit none
integer nd,md,qbd,n,m,qb
! reads dimension nd,md,qbd of B
nd=size(B,1)
md=size(B,2)
qbd=size(B,3)
! writes all the matrix to notePad
do qb=1,qbd
 do m=1,md
   do n=1,nd
   write(wout,*)"B(",n,m,qb,")=",B(n,m,qb)
   call print_it_
   end do
 end do
end do
end subroutine dump_matB

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! B0
!
subroutine alloc_matB0(n,stat)
use ponte_f_f
implicit none
integer n,stat
if (allocated(B0)) deallocate(B0)
allocate(B0(n,n),STAT=stat)
B0=0.0d0
end subroutine alloc_matB0

!
! to read the content of B0 from C routines
!
subroutine read_matB0(n,m,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m
x=B0(n,m)
end subroutine read_matB0

!
! to write to B0 from C routines
!
subroutine write_matB0(n,m,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m
B0(n,m)=x
end subroutine write_matB0

!
! test program
! dumps B0 to Main Window NotePad
!
subroutine dump_matB0()
use ponte_f_f
implicit none
integer nd,md,n,m
! reads dimension nd,md, of B0
nd=size(B0,1)
md=size(B0,2)
! writes all the matrix to notePad
 do m=1,md
   do n=1,nd
   write(wout,*)"B0(",n,m,")=",B0(n,m)
   call print_it_
   end do
 end do
end subroutine dump_matB0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BN1
!
subroutine alloc_matBN1(n,m,stat)
use ponte_f_f
implicit none
integer n,m,stat
if (allocated(BN1)) deallocate(BN1)
allocate(BN1(n,m),STAT=stat)
BN1=0.0d0
end subroutine alloc_matBN1

!
! to read the content of BN1 from C routines
!
subroutine read_matBN1(n,m,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m
x=BN1(n,m)
end subroutine read_matBN1

!
! to write to BN1 from C routines
!
subroutine write_matBN1(n,m,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m
BN1(n,m)=x
end subroutine write_matBN1

!
! test program
! dumps BN1 to Main Window NotePad
!
subroutine dump_matBN1()
use ponte_f_f
implicit none
integer nd,md,n,m
! reads dimension nd,md, of BN1
nd=size(BN1,1)
md=size(BN1,2)
! writes all the matrix to notePad
do m=1,md
   do n=1,nd
   write(wout,*)"BN1(",n,m,")=",BN1(n,m)
   call print_it_
   end do
end do
end subroutine dump_matBN1

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine alloc_matG(m,stat)
use ponte_f_f
implicit none
integer m,stat
if (allocated(G)) deallocate(G)
allocate(G(m,m),STAT=stat)
G=0.0d0
end subroutine alloc_matG

!
! to read the content of G from C routines
!
subroutine read_matG(n,m,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m
x=G(n,m)
end subroutine read_matG

!
! to write to G from C routines
!
subroutine write_matG(n,m,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m
G(n,m)=x
end subroutine write_matG

!
! test program
! dumps G to Main Window NotePad
!
subroutine dump_matG()
use ponte_f_f
implicit none
integer n,m,nd,md
! reads dimension nd,md of G
nd=size(G,1)
md=size(G,2)
! writes all the matrix to notePad
 do m=1,md
   do n=1,nd
   write(wout,*)"G(",n,m,")=",G(n,m)
   call print_it_
   end do
 end do
end subroutine dump_matG

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
subroutine alloc_matR(m,stat)
use ponte_f_f
implicit none
integer m,stat
if (allocated(R)) deallocate(R)
allocate(R(m,m),STAT=stat)
R=0.0d0
end subroutine alloc_matR

!
! to read the content of R from C routines
!
subroutine read_matR(n,m,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m
x=R(n,m)
end subroutine read_matR

!
! to write to R from C routines
!
subroutine write_matR(n,m,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m
R(n,m)=x
end subroutine write_matR

!
! test program
! dumps R to Main Window NotePad
!
subroutine dump_matR()
use ponte_f_f
implicit none
integer n,m,nd,md
! reads dimension nd,md of R
nd=size(R,1)
md=size(R,2)
! writes all the matrix to notePad
 do m=1,md
   do n=1,nd
   write(wout,*)"R(",n,m,")=",R(n,m)
   call print_it_
   end do
 end do
end subroutine dump_matR

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
subroutine alloc_matU(m,stat)
use ponte_f_f
implicit none
integer m,stat
if (allocated(U)) deallocate(U)
allocate(U(m,m),STAT=stat)
U=0.0d0
end subroutine alloc_matU

!
! to read the content of U from C routines
!
subroutine read_matU(n,m,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m
x=U(n,m)
end subroutine read_matU

!
! to write to U from C routines
!
subroutine write_matU(n,m,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m
U(n,m)=x
end subroutine write_matU

!
! test program
! dumps U to Main Window NotePad
!
subroutine dump_matU()
use ponte_f_f
implicit none
integer n,m,nd,md
! reads dimension nd,md of U
nd=size(U,1)
md=size(U,2)
! writes all the matrix to notePad
 do m=1,md
   do n=1,nd
   write(wout,*)"U(",n,m,")=",U(n,m)
   call print_it_
   end do
 end do
end subroutine dump_matU

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
subroutine alloc_matX0(m,stat)
use ponte_f_f
implicit none
integer m,stat
if (allocated(X0)) deallocate(X0)
allocate(X0(m,m),STAT=stat)
G=0.0d0
end subroutine alloc_matX0

!
! to read the content of X0 from C routines
!
subroutine read_matX0(n,m,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m
x=X0(n,m)
end subroutine read_matX0

!
! to write to X0 from C routines
!
subroutine write_matX0(n,m,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m
X0(n,m)=x
end subroutine write_matX0

!
! test program
! dumps X0 to Main Window NotePad
!
subroutine dump_matX0()
use ponte_f_f
implicit none
integer n,m,nd,md
! reads dimension nd,md of X0
nd=size(X0,1)
md=size(X0,2)
! writes all the matrix to notePad
 do m=1,md
   do n=1,nd
   write(wout,*)"X0(",n,m,")=",X0(n,m)
   call print_it_
   end do
 end do
end subroutine dump_matX0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
subroutine alloc_matPi0(m,stat)
use ponte_f_f
implicit none
integer m,stat
if (allocated(Pi0)) deallocate(Pi0)
allocate(Pi0(m),STAT=stat)
Pi0=0.0d0
end subroutine alloc_matPi0

!
! to read the content of Pi0 from C routines
!
subroutine read_matPi0(n,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n
x=Pi0(n)
end subroutine read_matPi0

!
! to write to Pi0 from C routines
!
subroutine write_matPi0(n,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n
Pi0(n)=x
end subroutine write_matPi0

!
! test program
! dumps Pi0 to Main Window NotePad
!
subroutine dump_matPi0()
use ponte_f_f
implicit none
integer n,nd
! reads dimension nd of Pi0
nd=size(Pi0,1)
! writes all the matrix to notePad
   do n=1,nd
   write(wout,*)"Pi0(",n,")=",Pi0(n)
   call print_it_
   end do
end subroutine dump_matPi0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
subroutine alloc_matPi(n,m,stat)
use ponte_f_f
implicit none
integer n,m,stat
if (allocated(Pi)) deallocate(Pi)
allocate(Pi(n,m),STAT=stat)
Pi=0.0d0
end subroutine alloc_matPi

!
! to read the content of Pi from C routines
!
subroutine read_matPi(n,m,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m
x=Pi(n,m)
end subroutine read_matPi

!
! to write to Pi from C routines
!
subroutine write_matPi(n,m,x)
use ponte_f_f
implicit none
real(dpr) :: x
integer n,m
Pi(n,m)=x
end subroutine write_matPi

!
! test program
! dumps Pi to Main Window NotePad
!
subroutine dump_matPi()
use ponte_f_f
implicit none
integer n,m,nd,md
! reads dimension nd,md of Pi
nd=size(Pi,1)
md=size(Pi,2)
! writes all the matrix to notePad
 do m=1,md
   do n=1,nd
   write(wout,*)"Pi(",n,m,")=",Pi(n,m)
   call print_it_
   end do
 end do
end subroutine dump_matPi

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! transfer to Fortran the Options Settings
! is_debug, is_verbose, is_timings
!
subroutine set_options(deb,ver,tim)
use ponte_f_f
implicit none
integer deb,ver,tim

if (deb == 1) then 
debug=.true.  
else 
debug=.false.  
endif

if (ver == 1) then
verb=.true.
else
verb=.false.
endif

if (tim == 1) then
timing=.true.
else
timing=.false.
endif



end subroutine set_options

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
subroutine transpose_matA
use ponte_f_f
implicit none
integer nd,md,qd,j
nd=size(A,1)
md=size(A,2)
qd=size(A,3)
! traspose works only if nd=md on the two first indexes
if (nd==md) then
do j=1,qd
     A(:,:,j)=transpose(A(:,:,j))
  end do
end if
end subroutine transpose_matA

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! to get A stochatic
!
subroutine normalize_matA
use ponte_f_f
implicit none
integer nd,md,j
real(dpr) :: ss
nd=size(A,1)
md=size(A,2)

! normalize works only if nd=md 
if (nd==md) then
do j=1,md
  ss =  sum(A(j,:,:)) 
  A(j,:,:)=A(j,:,:)/ss
end do

endif
end subroutine normalize_matA
