program Sparse_MxV_inv_off_fly_dp
use sep
use DataModule, only : dp
use sparse_matmult_fly_sym_dp 
use prec_dp
use solver_dp

implicit none
real (dp), dimension (:),allocatable :: m_inv,m_mig,xo
real (dp), dimension (:),  pointer   :: prec
real (dp), dimension (:),  pointer   :: Hcube
real     , dimension (:),allocatable :: auxi,auxi1
real (dp), dimension (:,:),allocatable :: remov,mmov
real     , dimension (:,:),allocatable :: auxi3
integer                              :: stat,n
logical                              :: adj,add, initial
integer                              :: niter,na,np  ! iterations
real                                 :: eps,alpha ! scaling parameter
real (dp)                            :: maxi
integer                              :: n1,n2,n3,n4,n5,n6,n7,n8,n9,n10 !nhx,nhy,naz,nax,nay,nhx,nhy,nz,nx,ny 
integer(kind=8)                      :: bn1,bn2,bn3,bn4,bn5,bn6,bn7,bn8,bn9,bn10 !nhx,nhy,naz,nax,nay,nhx,nhy,nz,nx,ny
real                                 :: o1,o2,o3,o4,o5,o6,o7,o8,o9,o10 
real                                 :: d1,d2,d3,d4,d5,d6,d7,d8,d9,d10
integer                              :: esize,i
integer                              :: ia
integer, external                    :: sreed
integer(kind=8)                      :: block

call sep_init()

call from_history('n1',n1); call from_history('d1',d1); call from_history('o1',o1)
call from_history('n2',n2); call from_history('d2',d2); call from_history('o2',o2)
call from_history('n3',n3); call from_history('d3',d3); call from_history('o3',o3)
call from_history('n4',n4); call from_history('d4',d4); call from_history('o4',o4)
call from_history('n5',n5); call from_history('d5',d5); call from_history('o5',o5)
call from_history('n6',n6); call from_history('d6',d6); call from_history('o6',o6)
call from_history('n7',n7); call from_history('d7',d7); call from_history('o7',o7)
call from_history('n8',n8); call from_history('d8',d8); call from_history('o8',o8)
call from_history('n9',n9); call from_history('d9',d9); call from_history('o9',o9)
call from_history('n10',n10); call from_history('d10',d10); call from_history('o10',o10)
call from_history("esize",esize)

call from_param('eps',eps,.00001)
call from_param('niter',niter,10)
call from_param('na',na,1)
call from_param('alpha',alpha,1.)
call from_param('initial',initial,.false.)
bn1=n1;bn2=n2;bn3=n3;bn4=n4;bn5=n5;bn6=n6;bn7=n7;bn8=n8;bn9=n9;bn10=n10
block=bn1*bn2*bn3*bn4*bn5*bn6*bn7*bn8*bn9*bn10
n=n6*n7*n8*n9*n10
write(0,*)"initial model",initial
write(0,*)"n1,n2,n3,n4,n5,n6,n7,n8,n9,n10",n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
write(0,*)"memory allocation";
write(0,*)"Hcube",dp*block/1024/1024/1024,"G"

allocate(Hcube(block),auxi(n1*n2*n3*n4*n5))

! Sep format Hessian matrix
do i=1,n 
if  (sreed("in",auxi,4*n1*n2*n3*n4*n5) /= 4*n1*n2*n3*n4*n5) &
    call seperr("Error reading Hessian")
    Hcube(1+(i-1)*n1*n2*n3*n4*n5:i*n1*n2*n3*n4*n5)=auxi
enddo
maxi=maxval(Hcube)
write(0,*)"LtL max", maxi

allocate (m_inv(n),m_mig(n),auxi1(n))
call sep_read(auxi1,"mig")  ! Target migration
m_mig=auxi1

allocate (prec(n))
call sep_read(auxi1,"prec")
prec=auxi1

allocate(remov(n,niter))
allocate(mmov (n,niter))

call smatmult_init_dp (Hcube,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10)
call prec_init_dp(prec)

  if (initial) then
    xo=m_mig
    write(0,*)"initial true",initial
  else
    xo=0.
  endif

  write(0,*)"prec",maxval(prec),minval(prec)
  write(0,*)"m_mig",maxval(m_mig),minval(m_mig)


  write(0,*)"solver_precond1"
  call  solver_precond1(m_inv,m_mig,smatmult_lop_dp,prec_lop_dp,xo,niter,mmov,remov)

  np=niter*na
  allocate(auxi3(n,niter))
  if( exist_file("remov")) then
     auxi3=remov
     call sep_write( auxi3, "remov")
     call to_history ("n1",n6,"remov")
     call to_history ("n2",n7,"remov")
     call to_history ("n3",n8,"remov")
     call to_history ("n4",n9,"remov")
     call to_history ("n5",n10,"remov")
     call to_history ("n6",np,"remov")
     call to_history ("o1",o6,"remov")
     call to_history ("o2",o7,"remov")
     call to_history ("o3",o8,"remov")
     call to_history ("o4",o9,"remov")
     call to_history ("o5",o10,"remov")
     call to_history ("o6",1,"remov")
     call to_history ("d1",d6,"remov")
     call to_history ("d2",d7,"remov")
     call to_history ("d3",d8,"remov")
     call to_history ("d4",d9,"remov")
     call to_history ("d5",d10,"remov")
     call to_history ("d6",1,"remov")
  end if
  if( exist_file("mmov")) then
     auxi3=mmov
     call sep_write( auxi3, "mmov")
     call to_history ("n1",n6,"mmov")
     call to_history ("n2",n7,"mmov")
     call to_history ("n3",n8,"mmov")
     call to_history ("n4",n9,"mmov")
     call to_history ("n5",n10,"mmov")
     call to_history ("n6",np,"mmov")
     call to_history ("o1",o6,"mmov")
     call to_history ("o2",o7,"mmov")
     call to_history ("o3",o8,"mmov")
     call to_history ("o4",o9,"mmov")
     call to_history ("o5",o10,"mmov")
     call to_history ("o6",1,"mmov")
     call to_history ("d1",d6,"mmov")
     call to_history ("d2",d7,"mmov")
     call to_history ("d3",d8,"mmov")
     call to_history ("d4",d9,"mmov")
     call to_history ("d5",d10,"mmov")
     call to_history ("d6",1,"mmov")
  end if

auxi1=m_inv
call sep_write(auxi1)

call to_history('n1',n6); call to_history('d1',d6); call to_history('o1',o6) ; call to_history('label1',"hx" )
call to_history('n2',n7); call to_history('d2',d7); call to_history('o2',o7) ; call to_history('label2',"hy")
call to_history('n3',n8); call to_history('d3',d8); call to_history('o3',o8) ; call to_history('label3',"z" )
call to_history('n4',n9); call to_history('d4',d9); call to_history('o4',o9) ; call to_history('label4',"x")
call to_history('n5',n10); call to_history('d5',d10); call to_history('o5',o10) ; call to_history('label5',"y")
call to_history('n6',1); call to_history('d6',1); call to_history('o6',0.) ; call to_history('label6'," ")
call to_history('n7',1); call to_history('d7',1); call to_history('o7',0.) ; call to_history('label7'," ")
call to_history('n8',1); call to_history('d8',1); call to_history('o8',0.) ; call to_history('label8'," ")
call to_history('n9',1); call to_history('d9',1); call to_history('o9',0.) ; call to_history('label9'," ")
call to_history('n10',1); call to_history('d10',1); call to_history('o10',0.) ; call to_history('label10'," ")

!deallocate(m_inv,m_mig,xo,auxi,auxi1,prec,Hcube)
!deallocate(mmov,remov,auxi3)


end program Sparse_MxV_inv_off_fly_dp
