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 (:),       pointer   :: Hcube
integer                              :: stat
logical                              :: adj,add 
real                                 :: maxi, mini
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)

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,*)"n1,n2,n3,n4,n5",n1,n2,n3,n4,n5
write(0,*)"n6,n7,n8,n9,n10",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,1+(i-1)*n1*n2*n3*n4*n5,i*n1*n2*n3*n4*n5
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

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

call smatmult_init (Hcube,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10)
stat = smatmult_lop (.false.,.false.)

STOP
end program 
