  program inv01

  use sep
  use cgstep_mod
  use adj_mod
  use solver_prc_mod
  use invhess
  use tcai1
  use igrad1

  implicit none
 
  integer					:: niter
  integer                         		:: i4,i6,i8,i9,nAop
  integer					:: n4,n6,n8,n9
  real                            		:: o4,o6,o8,o9,d4,d6,d8,d9,eps
  real,dimension(:),allocatable   		:: wgt,mig,inv

  call date_and_time(date,tempo)
  write(0,*) date,'  ',tempo

  call sep_init()
  call sep_begin_prog()
  !! get info from history file
  call from_history("n4",n4);  call from_history("n6",n6);  call from_history("n8",n8);  call from_history("n9",n9)
  call from_history("d4",d4);  call from_history("d6",d6);  call from_history("d8",d8);  call from_history("d9",d9)
  call from_history("o4",o4);  call from_history("o6",o6);  call from_history("o8",o8);  call from_history("o9",o9)

  !! get info from parameter
  call from_param("niter",niter,20);  call from_param("eps",eps,1.)

  !! write info out to history file
  call to_history("n1",n8);  call to_history("n2",n6);  call to_history("n3",n9)
  call to_history("d1",d8);  call to_history("d2",d6);  call to_history("d3",d9)
  call to_history("o1",o8);  call to_history("o2",o6);  call to_history("o3",o9)


!! allocate working space
  allocate(wgt(n4*n6*n8*n9),mig(n6*n8*n9),inv(n6*n8*n9))
!!
nAop=size(m)+size(aa)
call sreed("in",wgt,n4*n6*n8*n9*4)

do i9=1,n9
   do i8=1,n8
      do i6=1,n6
         inv=0.
	 if (mod(i3,10).eq.0) write(0,*) "CMP  ", i3, " of ", n3
	 ind1=(i9-1)*n8*n6+(i8-1)*n6+1; ind2=(i9-1)*n8*n6+(i8-1)*n6+n6
 	 ind3=(i9-1)*n8*n6*n4+(i8-1)*n4*n6+(i6-1)*n4+1;	 ind4=(i9-1)*n8*n6*n4+(i8-1)*n4*n6+(i6-1)*n4+n4
	 call invhess_init(i6,n4,wgt(ind3:ind4))
	 call solver_prc(m=inv, d=mig(ind1:ind2), &
	     Fop=avomodel_lop,Sop=igrad1_lop,stepper=cgstep,nSop=nAop,niter=niter,eps=eps,&
	     verb=.false.)
	 call invhess_close()
      end do  
   end do
end do

call srite ("out",inv,n6*n8*n9)
call cgstep_close()

deallocate(wgt,mig,inv)
call sep_end_prog()
end program inv01
