Least-square inversion with inexact adjoints. Method of conjugate directions: A tutorial

Next: EXAMPLES Up: Fomel: Conjugate directions Previous: Short memory of the

# PROGRAM

The program in Table 1 implements one iteration of the conjugate-direction method. It is based upon Jon Claerbout's cgstep() program (, ) and uses an analogous naming convention. Vectors in the data space are denoted by double letters.

 void sf_cdstep(bool forget /* restart flag */, int nx /* model size */, int ny /* data size */, float* x /* current model [nx] */, const float* g /* gradient [nx] */, float* rr /* data residual [ny] */, const float* gg /* conjugate gradient [ny] */) /*< Step of conjugate-direction iteration. The data residual is rr = A x - dat >*/ { float *s, *si, *ss; double alpha, beta; int i, n, ix, iy; s = sf_floatalloc(nx+ny); ss = s+nx; for (ix=0; ix < nx; ix++) { s[ix] = g[ix]; } for (iy=0; iy < ny; iy++) { ss[iy] = gg[iy]; } sf_llist_rewind (steps); n = sf_llist_depth (steps); for (i=0; i < n; i++) { sf_llist_down (steps, &si, &beta); alpha = - cblas_dsdot(ny, gg, 1, si+nx, 1) / beta; cblas_saxpy(nx+ny,alpha,si,1,s,1); } beta = cblas_dsdot(ny, s+nx, 1, s+nx, 1); if (beta < DBL_EPSILON) return; sf_llist_add (steps, s, beta); if (forget) sf_llist_chop (steps); alpha = - cblas_dsdot(ny, rr, 1, ss, 1) / beta; cblas_saxpy(nx,alpha,s, 1,x, 1); cblas_saxpy(ny,alpha,ss,1,rr,1); } 

Table 1. The source of this program is RSF/api/c/cdstep.c

In addition to the previous steps and their conjugate counterparts (array s), the program stores the squared norms (variable beta) to avoid recomputation. For practical reasons, the number of remembered iterations can actually be smaller than the total number of iterations.

 Least-square inversion with inexact adjoints. Method of conjugate directions: A tutorial

Next: EXAMPLES Up: Fomel: Conjugate directions Previous: Short memory of the

2013-03-03