!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  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                                   !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  examplesf.f90                                              !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! setup input arrays for every Example Problem
! to be called from callbacks.c only
!
subroutine examples_setup (num, m, n, qa, qb, stat,alpha, beta )
use ponte_f_f
! logical debug global 
implicit none
integer num   !Example number
integer m     !array first dimension
integer n     !array second dimension
integer qa     !Aq mxn arrays
integer qb    !Bqb mxn arrays 0= no Bqb arrays
integer stat  !status output 0=ok, else =no memory 
real(dp) :: alpha,beta ! optional parameter for some examples
!
integer i,j,h
real(dp) s
real(dp),allocatable,dimension(:)::v
real(dp):: binomial

if (debug) write(*,*) "examples_setup: entering with num=",num," m=",m," n=",n," qa=",qa," qb=",qb," alpha=",alpha," beta=",beta

! make space in memory
if (allocated(A)) deallocate(A)
if (allocated(B)) deallocate(B)
if (allocated(B0)) deallocate(B0)
if (allocated(BN1)) deallocate(BN1)
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)

stat=1000

! qa = number of A blocks 
! A blocks are square m x m
! qb = number of B blocks
! B0 is n x n
! BN1 is m x n in QBD and MG1 and is n x m in GIM1
! B blocks are n x m in QBD and MG1 and are m x n in GIM1
! alpha beta parameters
! this example setup is insufficient for creating practical examples.
! something better must be done in next version
!




! example 2 begins
if (num == 2) then 
! n,qa, not used here ... square and qa=3
stat = 0

call alloc_matA(m,3,stat)
if (debug) write (*,*) "examples_setup: Allocating A, stat=",stat
if (stat /= 0) return 
call alloc_matB(m,m,1,stat)
if (debug) write (*,*) "examples_setup: Allocating B, stat=",stat
if (stat /= 0) return
call alloc_matB0(m,stat)
if (debug) write (*,*) "examples_setup: Allocating B0, stat=",stat
if (stat /= 0) return
call alloc_matBN1(m,m,stat)
if (debug) write (*,*) "examples_setup: Allocating BN1, stat=",stat
if (stat /= 0) return


call alloc_matG(m,stat)
if (debug) write (*,*) "examples_setup: Allocating G, stat=",stat
if (stat /= 0 ) return
if (debug) write (*,*) "examples_setup: Now computing the Matrices"
do i=1,m
do j=1,m
!A(i,j,1)=1.d0+abs(i-j)
A(i,j,1)=m+6.d0-abs(i-j)
end do
end do
A(:,:,2)=0.d0
A(:,:,3)=A(:,:,1)
A(1,1,1)=A(1,1,1)+alpha
G=A(:,:,1)+A(:,:,2)+A(:,:,3)
do i=1,m
s=sum(G(i,:))
if (s == 0.0) then
stat=2
return
end if
A(i,:,1)=A(i,:,1)/s
A(i,:,2)=A(i,:,2)/s
A(i,:,3)=A(i,:,3)/s
end do
if (debug) write (*,*) "examples_setup: Dellocating G if allocated"
if (allocated(G)) deallocate(G)

! default B for this example - no user changes on B
B0=A(:,:,1)+A(:,:,2)
BN1=A(:,:,1)
B(:,:,1)=A(:,:,3)

if (debug) write (*,*) "examples_setup: All done for example 1"

endif
! example 2 ends
!
! example 1 begins
! Am1 rank 1
if (num == 1) then
! n,qa,qb,B not used here
stat = 0

call alloc_matA(m,3,stat)
if (debug) write (*,*) "examples_setup: Allocating A, stat=",stat
if (stat /= 0) return
call alloc_matG(m,stat)
if (debug) write (*,*) "examples_setup: Allocating G, stat=",stat
if (stat /= 0 ) return

call alloc_matB(m,m,1,stat)
if (debug) write (*,*) "examples_setup: Allocating B, stat=",stat
if (stat /= 0) return
call alloc_matB0(m,stat)
if (debug) write (*,*) "examples_setup: Allocating B0, stat=",stat
if (stat /= 0) return
call alloc_matBN1(m,m,stat)
if (debug) write (*,*) "examples_setup: Allocating BN1, stat=",stat
if (stat /= 0) return

if (debug) write (*,*) "examples_setup: Now computing the Matrices"
!
do i=1,m
   A(:,i,1)=0.8d0*i
end do
do i=1,m
   do j=1,m
      A(i,j,2)=(i+j)/(2.d0*m)
      A(i,j,3)=abs(i-j)/(4.d0*m)
   end do
end do

! normalization
G=A(:,:,1)+A(:,:,2)+A(:,:,3)
do i=1,m
s=sum(G(i,:))
if (s == 0.0) then
stat=2
return
end if
A(i,:,1)=A(i,:,1)/s
A(i,:,2)=A(i,:,2)/s
A(i,:,3)=A(i,:,3)/s
end do

! default B for this example - no user changes on B
B0=A(:,:,1)+A(:,:,2)
BN1=A(:,:,1)
B(:,:,1)=A(:,:,3)

if (allocated(G)) deallocate(G)
if (debug) write (*,*) "examples_setup: Dellocating G if allocated"
if (debug) write (*,*) "examples_setup: All done for example 2"
endif
! example 1  ends
!
!  example 3 begins
if (num == 3) then
n=m-1

! tridiagonal matrices
! qa,qb,not used here
stat = 0

call alloc_matA(m,3,stat)
if (debug) write (*,*) "examples_setup: Allocating A, stat=",stat
if (stat /= 0) return
call alloc_matG(m,stat)
if (debug) write (*,*) "examples_setup: Allocating G, stat=",stat
if (stat /= 0 ) return

call alloc_matB(n,m,1,stat)
if (debug) write (*,*) "examples_setup: Allocating B, stat=",stat
if (stat /= 0) return
call alloc_matB0(n,stat)
if (debug) write (*,*) "examples_setup: Allocating B0, stat=",stat
if (stat /= 0) return
call alloc_matBN1(m,n,stat)
if (debug) write (*,*) "examples_setup: Allocating BN1, stat=",stat
if (stat /= 0) return

if (debug) write (*,*) "examples_setup: Now computing the Matrices"
!
A(:,:,1)=0.d0
A(:,:,2)=0.d0
A(:,:,3)=0.d0
do i=1,m
   do j=1,m
      if (abs(i-j).le.1) then
         A(i,j,1)=abs(2.d0*i-3.d0*j)
         A(i,j,2)=(i+j)/(2.d0*m)
         A(i,j,3)=abs(i-j)/(4.d0*m)
      end if
   end do
end do

! normalization
G=A(:,:,1)+A(:,:,2)+A(:,:,3)
do i=1,m
s=sum(G(i,:))
if (s == 0.0) then
stat=2
return
end if
A(i,:,1)=A(i,:,1)/s
A(i,:,2)=A(i,:,2)/s
A(i,:,3)=A(i,:,3)/s
end do

! default B for this example - no user changes on B
B0=A(1:n,1:n,1)
BN1=A(1:m,2:m,1)
B(1:n,1:m,1)=A(1:n,1:m,3)
do i=1,n
   s=sum(B0(i,:))+sum(B(i,:,1))
   B0(i,:)=B0(i,:)/s
   B(i,:,1)=B(i,:,1)/s
end do

   BN1(:,1)=BN1(:,1)+A(:,1,1)
!
!


if (allocated(G)) deallocate(G)
if (debug) write (*,*) "examples_setup: Dellocating G if allocated"
if (debug) write (*,*) "examples_setup: All done for example 3"
endif
! example 3 ends
!

!
! example 4 begins
! system is pos. recurrent
! in this example m=n=qa=21, qb,B not used here
! a=a(m+1,m+1,m+1)
! default alpha=0.001; d=12
! 
if (num==4) then
stat = 0
qb=qa-2
call alloc_matA(m,qa,stat)
if (debug) write (*,*) "examples_setup: Allocating A, stat=",stat
if (stat /= 0) return
call alloc_matG(m,stat)
if (debug) write (*,*) "examples_setup: Allocating G, stat=",stat
if (stat /= 0 ) return
call alloc_matR(m,stat)
if (debug) write (*,*) "examples_setup: Allocating R, stat=",stat
if (stat /= 0 ) return

call alloc_matB(m,m,qb,stat)
if (debug) write (*,*) "examples_setup: Allocating A, stat=",stat
if (stat /= 0) return
call alloc_matB0(m,stat)
if (debug) write (*,*) "examples_setup: Allocating G, stat=",stat
if (stat /= 0 ) return

if (debug) write (*,*) "examples_setup: Now computing the Matrices"
h=m-1
A=0.d0
G=0.0d0
G(1,1)=1.d0-h*alpha
  G(1,2)=h*alpha
  do i=2, h
     G(i,i-1)=(i-1.d0)*alpha
     G(i,i+1)=(h-i+1.d0)*alpha
     G(i,i)=1.d0-G(i,i-1)-G(i,i+1);
  end do
  G(h+1,h)=alpha*h
  G(h+1,h+1)=1.d0-alpha*h
  do j=0,h
     R=0.d0
     do i=0,h
        if (i >= j) then
           R(i+1,i+1)=binomial(i,j)*(1.d0/beta)**j*((beta-1.d0)/beta)**(i-j);
        end if
     end do
     if (sum(matmul(R,G))> 10**(-100)) then
        A(1:m,1:m,j+1)=matmul(R,G)
     end if
  end do
do h=1,qa
   do i=1,m
      do j=1,m
         if (A(i,j,h)<0.d0) A(i,j,h)=0.d0
      end do
   end do
end do
 b0=a(:,:,1)+a(:,:,2)
 B=A(:,:,3:qa)

if (allocated(G)) deallocate(G)
if (debug) write (*,*) "examples_setup: Dellocating G if allocated"
if (allocated(R)) deallocate(R)
if (debug) write (*,*) "examples_setup: Dellocating R if allocated"


if (debug) write (*,*) "examples_setup: All done for example 4 - no BN1 here"
endif
! example 4 ends
!


! example 5 begins
! system is pos. recurrent if .....
! in this example m=n=qa=21, qb,B not used here
! a=a(m+1,m+1,m+1)
! default alpha=0.001; d=12
! 
if (num==5) then
stat = 0

call alloc_matA(m,qa,stat)
if (debug) write (*,*) "examples_setup: Allocating A, stat=",stat
if (stat /= 0) return
call alloc_matG(m,stat)
if (debug) write (*,*) "examples_setup: Allocating G, stat=",stat
if (stat /= 0 ) return
call alloc_matR(m,stat)
if (debug) write (*,*) "examples_setup: Allocating R, stat=",stat
if (stat /= 0 ) return

call alloc_matB(n,m,qb,stat)
if (debug) write (*,*) "examples_setup: Allocating A, stat=",stat
if (stat /= 0) return
call alloc_matB0(n,stat)
if (debug) write (*,*) "examples_setup: Allocating G, stat=",stat
if (stat /= 0 ) return
call alloc_matBN1(m,n,stat)
if (debug) write (*,*) "examples_setup: Allocating BN1, stat=",stat
if (stat /= 0) return

if (debug) write (*,*) "examples_setup: Now computing the Matrices"
A=0.d0
A(1,3,5)=1.d0;A(2,1,5)=1.d0;A(3,2,5)=1.d0
A(1,1,4)=1.d0;A(2,2,4)=2.d0;A(3,3,4)=3.d0
A(:,:,3)=A(:,:,4)+A(:,:,5)
A(:,:,2)=A(:,:,3)+A(:,:,4)
A(:,:,1)=A(:,:,3)+A(:,:,2)+2.01d0*A(:,:,5)
G=sum(A,dim=3)
do i=1,m
s=sum(G(i,:))
A(i,:,1)=A(i,:,1)/s
A(i,:,2)=A(i,:,2)/s
A(i,:,3)=A(i,:,3)/s
A(i,:,4)=A(i,:,4)/s
A(i,:,5)=A(i,:,5)/s
end do
b=1.d0*(1.d0/(qb*m+n))
b0=1.d0*(1.d0/(qb*m+n))
!bn1=sum(a(1,:,1))
do i=1,m
bn1(i,:)=sum(a(i,:,1))/2.d0
end do
if (allocated(G)) deallocate(G)
if (debug) write (*,*) "examples_setup: Dellocating G if allocated"
if (allocated(R)) deallocate(R)
if (debug) write (*,*) "examples_setup: Dellocating R if allocated"
if (debug) write (*,*) "examples_setup: All done for example 5"
endif
! example 5 ends
!



! example 6 begins
if (num==6) then !GM1 transient
stat = 0

m=3
qa=6

call alloc_matA(m,qa,stat)
if (debug) write (*,*) "examples_setup: Allocating A, stat=",stat
if (stat /= 0) return
a=0.d0
A(1,:,1)=1.d0/4.d0
A(2,:,1)=1.d0/4.d0; a(2,2,1)=3.d0/16.d0
A(3,2,1)=1.d0/4; A(3,3,1)=1.d0/2;

A(1,1,2)=1.d0/8.d0; A(2,2,2)=1.d0/16.d0; A(3,3,2)=1.d0/8.d0

A(1,2,6)=1.d0/8.d0; A(2,1,6)=1.d0/8.d0; A(2,3,6)=1.d0/8.d0; A(3,2,6)=1.d0/8.d0
if (debug) write (*,*) "examples_setup: All done for example 3 - no B here"
endif
! example 6 ends
!

!
! example 7 begins
if (num==7) then !GM1 positive recurrent
stat = 0

call alloc_matA(m,qa,stat)
if (debug) write (*,*) "examples_setup: Allocating A, stat=",stat
if (stat /= 0) return
a=0.d0

A(1,1,1)=1.d0/8.d0;A(2,2,1)=1.d0/4.d0;A(3,3,1)=3.d0/8.d0
A(1,2,2)=1.d0/8.d0;A(2,1,2)=1.d0/8.d0;A(2,3,2)=1.d0/8.d0;A(3,2,2)=1.d0/8.d0
A(1,:,6)=1.d0/4.d0
A(2,1,6)=1.d0/4.d0;A(2,3,6)=1.d0/4.d0;A(3,2,6)=1.d0/4.d0;A(3,3,6)=1.d0/4.d0

if (allocated(v)) deallocate(v)
if (allocated(B)) deallocate(B)
if (allocated(B0)) deallocate(B0)
allocate(b(m,m,qb),B0(m,m),v(m),stat=stat)

if (stat /= 0 ) return
b=0.d0
B0=0.d0

v=1.d0-sum(a(:,:,1),dim=2)
b0(:,1)=v
do i=1,qb
   v=v-sum(a(:,:,i+1),dim=2)
   b(:,1,i)=v
end do


if (debug) write (*,*) "examples_setup: All done for example 8"
endif
! example 7 ends
!
! example 8
if (num==8) then !QBD
stat = 0
! make space in memory
if (allocated(A)) deallocate(A)
if (allocated(B)) deallocate(B)
call alloc_matB0(n,stat)
if (debug) write (*,*) "examples_setup: Allocating G, stat=",stat
if (stat /= 0 ) return
call alloc_matBN1(m,n,stat)
if (debug) write (*,*) "examples_setup: Allocating BN1, stat=",stat
if (stat /= 0) return

call alloc_matA(m,qa,stat)

if (debug) write (*,*) "examples_setup: Allocating A, stat=",stat
if (stat /= 0) return

if (allocated(v)) deallocate(v)
allocate(v(m),stat=stat)


a=0.d0

A(10,1,3)=1.d0/8.d0;A(11,2,3)=1.d0/4.d0;A(12,3,3)=3.d0/8.d0



A(1,2,2)=1.d0/8.d0;A(2,1,2)=1.d0/8.d0;A(2,3,2)=1.d0/8.d0;A(3,2,2)=1.d0/8.d0
A(4,5,2)=1.d0/8.d0;A(5,4,2)=1.d0/8.d0;A(5,6,2)=1.d0/8.d0;A(6,5,2)=1.d0/8.d0
A(7,8,2)=1.d0/8.d0;A(8,7,2)=1.d0/8.d0;A(8,9,2)=1.d0/8.d0;A(9,8,2)=1.d0/8.d0
A(10,11,2)=1.d0/8.d0;A(11,10,2)=1.d0/8.d0;A(11,12,2)=1.d0/8.d0;A(12,11,2)=1.d0/8.d0


A(1,4,2)=1.d0/8.d0; A(2,5,2)=1.d0/4.d0; A(3,6,2)=3.d0/8.d0; 
A(4,7,2)=1.d0/8.d0; A(5,8,2)=1.d0/4.d0; A(6,9,2)=3.d0/8.d0; 
A(7,10,2)=1.d0/8.d0; A(8,11,2)=1.d0/4.d0; A(9,12,2)=3.d0/8.d0; 

A(1,1:3,1)=1.d0/4.d0
A(2,1,1)=1.d0/4.d0;A(2,3,1)=1.d0/4.d0;
A(3,2,1)=1.d0/4.d0;A(3,3,1)=1.d0/4.d0;

A(4,4:6,1)=1.d0/4.d0
A(5,4,1)=1.d0/4.d0;A(5,6,1)=1.d0/4.d0;
A(6,5,1)=1.d0/4.d0;A(6,6,1)=1.d0/4.d0;

A(7,7:9,1)=1.d0/4.d0
A(8,7,1)=1.d0/4.d0;A(8,9,1)=1.d0/4.d0;
A(9,8,1)=1.d0/4.d0;A(9,9,1)=1.d0/4.d0;

A(10,10:12,1)=1.d0/4.d0
A(11,10,1)=1.d0/4.d0;A(11,12,1)=1.d0/4.d0;
A(12,11,1)=1.d0/4.d0;A(12,12,1)=1.d0/4.d0;

B0=A(:,:,2);BN1=A(:,:,1)
B0(1:12,1:3)=0.d0
Bn1(1:3,1:3)=0.d0
v=1.d0-sum(a(1:3,4:6,2),dim=2)
b0(1:3,1)=v

v=v-sum(a(4:6,4:6,2),dim=2)
b0(4:6,1)=v

v=v-sum(a(7:9,4:6,2),dim=2)
b0(7:9,1)=v

v=v-sum(a(10:12,4:6,2),dim=2)
b0(10:12,1)=v

v=v-sum(a(1:3,4:6,1),dim=2)
bn1(1:3,1)=v
endif
! example 8 ends

if(debug) then
write(*,*) "examples_setup: exiting with num=",num," m=",m," n=",n," qa=",qa," qb=",qb," alpha=",alpha," beta=",beta
end if
end subroutine  examples_setup

