Module oc_lsqrpre_mod

1.

Purpose: preconditioned LSQR solver

$$\fbox {$ \displaystyle \left [\mathcal W\mathcal L\P \atop\e... ...ht ] \bold p\approx \left [\mathcal W\bold D\atop0\right ] $}$$

		 $\bold p= 0$ 
		 $\left [\Ud\atop\Um\right ] = \left [\mathcal W\bold D\atop0\right ]$ $\b = \sqrt{\parallel \bold U\parallel^2 }$		$\bold U=\frac{1}{\b}\bold U$		 $\bold v= \left [\P^*\mathcal L^*\mathcal W^*\quad\epsilon I\right ] \left [\Ud\atop\Um\right ]$				$\alpha= \sqrt{\parallel \bold v\parallel^2 }$		$\bold v=\frac{1}{\
a}\bold v$		 $\bold w= \bold v$ 
		 $\bar\rho =\alpha$		$\bar\phi = \b$ 
		 iterate { 
		 		  $\bold U= -\alpha\bold U$ $\left [\Ud\atop\Um\right ] =\left [\Ud\atop\Um\right ] +\left [\mathcal W\mathcal L\P\atop\epsilon I\right ]\bold v$		 		  $\b = \sqrt{\parallel \bold U\parallel^2 }$ $\bold U=\frac{1}{\b}\bold U$
		 		  $\bold v= -\b \bold v$ $\bold v=\bold v+\left [\P^*\mathcal L^*\mathcal W^*\quad\epsilon I\right ]\left [\Ud\atop\Um\right ]$		 		  $\alpha= \sqrt{\parallel \bold v\parallel^2 }$ $\bold v=\frac{1}{\alpha}\bold v$
		 		  $\rho = \sqrt{\bar\rho^2 + \b^2}$ 
		 		  $c=\frac{\bar\rho}{\rho}$ $\bar\rho=-c\alpha$ 
		 		  $s=\frac{\b}{\rho}$ $\theta = s\alpha$ 
		 		  $\phi = c\bar\phi$ $\bar\phi=s\bar\phi$ 
		 		  $t_1=\frac{\phi}{\rho}$ $t_2=-\frac{\theta}{\rho}$
		 		  $\bold p= \bold p+ t_1 \bold w$ 
		 		  $\bold w= \bold v+ t_2 \bold w$ 
		 		 }
		 $\bold m=\P\bold p$

2.

Functions and subroutines

(a)

subroutine oc_lsqrpre_init(niter,eps,maxmem,verb,movie)

Purpose: initialize the preconditioned LSQR solver

• niter: iterations number
• eps: scaling factor
• verb: verbose flag (optional)
• movie: movie output flag (optional)

(b)

subroutine oc_lsqrpre( L,P, x_,yy_,npre,W,op1 ...op9)

Purpose: preconditioned LSQR solver

• L: out-of-core linear operator
• P: out-of-core preconditioning operator
• npre: dimension of the preconditioning output
• W: out-of-core weighting operator

E  

m.H:
	Spike n1=50 n2=30 nsp=6  	k1=17,34,27,20,33,41 l1=17,34,27,20,33,41   	k2=1,2,15,18,29,30   l2=1,2,15,18,29,30     	mag=-1,1,1,-1,1,-1 >$@

d.H: m.H $B/OCsimple
	< m.H Window >d.H
	OCsimple operation=1 model=m.H data=d.H maxmem=600 >$n

view1: m.H d.H
	< m.H Merge axis=3 space=n d.H | Window | 	Grey pclip=100 gainpanel=a eout=1 > j.T
	< j.T Grey color=g bias=0.5 | Tube

i.H: d.H $B/OCsimple
	< m.H Window | Scale rscale=0. >$@
	OCsimple operation=2 model=i.H data=d.H  	maxmem=600 niter=20 eps=0.02 verb=y >$n

view2: d.H i.H m.H
	< m.H Merge axis=3 space=n d.H i.H| Window | 	Grey pclip=100 gainpanel=a eout=1 > j.T
	< j.T Grey color=g bias=0.5  	titles=model:data:inversion | Tube

dot: m.H d.H $B/OCsimple
	OCsimple operation=0 model=m.H data=d.H maxmem=600 >$n

F  

program OCsimple
  use sep
  use oc_global_mod
  use oc_file_mod
  use oc_dottest_mod
  use oc_cgstep_mod
  use oc_solver_mod
  use oc_laplacian_mod

implicit none
  logical             :: verb
  character(len=128)  :: x_,yy_,t_, name
  integer             :: maxmem, stat, niter, nf, operation
  type(fileinfo)      :: file
  type(cfilter)       :: aa
  real                :: eps

call sep_init()
  call from_param("operation",operation,0)
  call from_param("maxmem",maxmem)
  call from_param("verb",verb,.false.)
  call from_param("nf",nf,5)
  call from_param("niter",niter,10)
  call from_param("eps",eps,1.0)
  x_= "model"; call auxinout(x_ )
  yy_="data" ; call auxinout(yy_)
  name="test.H"; t_=oc_clone(F, x_,name,maxmem)
  call sep_close()
  !! operator init
  call oc_allocatefile(file, x_, maxmem)
  call oc_infofile(file)
  do while(2*nf+1 > file%n(1))
     nf=nf-1
  end do
  call oc_deallocatefile(file)
  call oc_laplacian_init(x_,nf,10,0.0,maxmem)

select case(operation)
  case(0) !! dot product test
     call oc_dottest_init(maxmem=maxmem)
     call oc_dottest(oc_laplacian, x_,yy_)
  case(1) !! simple forward operator
     stat = oc_laplacian(F,F,x_,yy_)
  case(2) !! inversion
     call oc_solver_init(niter,maxmem,verb)
     call oc_solver(oc_laplacian,oc_cgstep,x_,yy_)
  case default
     call seperr("missing operation")
  end select

call exit (0)
end program OCsimple