program Sparse_Mxv_inv_newsolver_sym_off
use sep
use dottest 
use sparse_matmult_fly_sym 
use scaled_identity
use solver_H_mod
use chain_weight
use myarray_mod
use differential_off_mod

implicit none
real, dimension (:),     allocatable :: m_inv,m_mig,xo,resw,auxi
real, dimension(:),      pointer     :: w_out,wto
real, dimension (:,:),   allocatable :: remov,mmov
real, dimension (:),       pointer   :: Hcube
real, dimension (2)                  :: dot1,dot2
integer                              :: stat
logical                              :: adj,add, reg_off,initial,reg_damp,dot
integer                              :: niter,na,np  ! iterations
real                                 :: eps,alpha,beta_max ! scaling parameter
real                                 :: maxi
real, dimension(:),      pointer     :: err       ! residual error per iteration
real, dimension(:),      pointer     :: res       ! final residual error
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('reg_off',reg_off,.true.)
call from_param('initial',initial,.false.)
call from_param('dot',dot,.false.)
call from_param('beta_max',beta_max,1.)
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
write(0,*)"offset regularization", reg_off
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",4*block/1024/1024/1024,"G"

allocate(Hcube(block))

! Sep format Hessian matrix
do i=1,n6*n7*n8*n9*n10
write(0,*)i
if  (sreed("in",Hcube(1+(i-1)*n1*n2*n3*n4*n5:i*n1*n2*n3*n4*n5),4*n1*n2*n3*n4*n5) /= 4*n1*n2*n3*n4*n5) &
    call seperr("Error reading Hessian")
enddo

STOP

maxi=maxval(Hcube)
write(0,*)"LtL max", maxi
Hcube=(Hcube)

allocate(m_inv(n6*n7*n8*n9*n10),m_mig(n6*n7*n8*n9*n10),xo(n6*n7*n8*n9*n10),wto(n6*n7*n8*n9*n10),res(n6*n7*n8*n9*n10))
if( exist_file("remov")) allocate(remov(n6*n7*n8*n9*n10,niter))
if( exist_file("mmov")) allocate(mmov(n6*n7*n8*n9*n10,niter))
if( exist_file("resw")) allocate(resw(n6*n7*n8*n9*n10))
if( exist_file("prec")) allocate(w_out(n6*n7*n8*n9*n10),auxi(n6*n7*n8*n9*n10))
call sep_read(m_mig,"mig")  ! Target migration
if( exist_file("resw")) then
  call sep_read(resw,"resw")
  wto=resw
  write(0,*)"read weight"
else
  wto=1.  
endif 
if( exist_file("prec")) then
  call sep_read(auxi,"prec")
  w_out=auxi
endif

if( exist_file("prec").and.exist_file("resw")) then
 where ( w_out > 5. )
   wto=0.
 end where
endif

call smatmult_init (Hcube,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10)
if( exist_file("prec")) call scaled_identity_v_init(w_out)
call differential_off_init(n1,d1,o1,eps)
call  myarray_init (eps) 

if (dot) call dot_test(myarray,n6*n7*n8*n9*n10,n6*n7*n8*n9*n10,dot1,dot2)
!call dot_test(differential_off2_exec,n6*n7*n8*n9*n10,n6*n7*n8*n9*n10,dot1,dot2)
!call dot_test(smatmult_lop,n6*n7*n8*n9*n10,n6*n7*n8*n9*n10,dot1,dot2)


do ia=1,na
  if (ia==1) then
   if (initial) then
     xo=m_mig
     write(0,*)"initial true",initial
   else
     xo=0.
   endif
  else 
   if( exist_file("resw")) then
    wto=(1./sqrt(sqrt(1.+(remov(:,niter)/alpha)**2)))*resw
   else
    wto=1./sqrt(sqrt(1.+(remov(:,niter)/alpha)**2))
   endif
   xo=m_inv
  endif
if ( exist_file("resw")) then
  call scaled_identity_wto_init( wto )
end if

if( exist_file("prec"))  write(0,*)"w_out",maxval(w_out),minval(w_out)
  write(0,*)"m_mig",maxval(m_mig),minval(m_mig)
  write(0,*)"wto",maxval(wto),minval(wto)

if( exist_file("prec").and.exist_file("resw")) then
  write(0,*)"solver_precond2"
  call  solver_precond2(m_inv,m_mig,myarray,scaled_identity_v_lop,xo,niter,wto,mmov,remov,beta_max)
elseif( exist_file("prec")) then
  write(0,*)"solver_precond1"
  call  solver_precond1(m_inv,m_mig,myarray,scaled_identity_v_lop,xo,niter,mmov,remov,beta_max)
elseif  (exist_file("resw")) then
  write(0,*)"solver_tiny w"
  call  solver_tiny_w(m_inv,m_mig,chain_w,niter,wto,mmov,remov)
else
  write(0,*)"solver_tiny"
  call  solver_tiny(m_inv,m_mig,myarray,niter,mmov,remov)
endif

  np=niter*na
  if( exist_file("remov")) then
     call sep_write( remov, "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
     call sep_write( mmov, "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
  if( exist_file("wmov")) then
     call sep_write( wto, "wmov")
     call to_history ("n1",n6,"wmov")
     call to_history ("n2",n7,"wmov")
     call to_history ("n3",n8,"wmov")
     call to_history ("n4",n9,"wmov")
     call to_history ("n5",n10,"wmov")
     call to_history ("n6",na,"wmov")
     call to_history ("o1",o6,"wmov")
     call to_history ("o2",o7,"wmov")
     call to_history ("o3",o8,"wmov")
     call to_history ("o4",o9,"wmov")
     call to_history ("o5",o10,"wmov")
     call to_history ("o6",1,"wmov")
     call to_history ("d1",d6,"wmov")
     call to_history ("d2",d7,"wmov")
     call to_history ("d3",d8,"wmov")
     call to_history ("d4",d9,"wmov")
     call to_history ("d5",d10,"wmov")
     call to_history ("d6",1,"wmov")
  end if
end do


call sep_write(m_inv)

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'," ")


if( exist_file("remov")) deallocate(remov)
if( exist_file("mmov")) deallocate(mmov)
if( exist_file("resw")) deallocate(resw)
deallocate(m_inv,m_mig,xo,wto)


end program 
