program Sparse_MxV_inv_off
use sep
use dottest 
use sparse_matmult 
use cgstep_mod
use solver_mod_mio
use scaled_identity
use sparse_mapping_off

implicit none
real, dimension (:),     allocatable :: m_inv,m_mig
real, dimension (:),       pointer   :: Hcube
real, dimension (:,:)    , pointer   :: Hsp
real, dimension (2)                  :: dot1,dot2
integer                              :: stat
logical                              :: adj,add
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 !naz,nax,nay,nhx,nhy,nz,nx,ny,nhx,nhy 
real                                 :: o1,o2,o3,o4,o5,o6,o7,o8,o9,o10 
real                                 :: d1,d2,d3,d4,d5,d6,d7,d8,d9,d10
integer                              :: esize
integer                              :: ia,cont

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 sp_cont (n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,cont)
write(0,*)'cont',cont

allocate(Hcube(n1*n2*n3*n4*n5*n6*n7*n8*n9*n10),Hsp(cont,3))

call sep_read(Hcube)                 ! Sep format Hessian matrix

call sp_mapp (Hcube,Hsp,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10) ! Mapp Heesian matrix from Sep to sparse representation

deallocate(Hcube)

allocate(m_inv(n6*n7*n8*n9*n10),m_mig(n6*n7*n8*n9*n10))

call sep_read(m_mig,"mig")  ! Target migration

call smatmult_init (Hsp)
call dot_test(smatmult_lop,n6*n7*n8*n9*n10,n6*n7*n8*n9*n10,dot1,dot2)

stat = smatmult_lop(.false.,.false.,m_mig, m_inv)

write(0,*)'max(m_mig)',maxval(m_mig)
write(0,*)'max(m_inv)',maxval(m_inv)

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


deallocate(m_inv,m_mig,Hsp)


end program Sparse_MxV_inv_off
