Least-square inversion with inexact adjoints. Method of conjugate directions: A tutorial |
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); } |
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 |