Actual source code: lsqr.c

petsc-master 2018-05-24
Report Typos and Errors

  2: /* lourens.vanzanen@shell.com contributed the standard error estimates of the solution, Jul 25, 2006 */
  3: /* Bas van't Hof contributed the preconditioned aspects Feb 10, 2010 */

  5: #define SWAP(a,b,c) { c = a; a = b; b = c; }

  7:  #include <petsc/private/kspimpl.h>

  9: typedef struct {
 10:   PetscInt  nwork_n,nwork_m;
 11:   Vec       *vwork_m;   /* work vectors of length m, where the system is size m x n */
 12:   Vec       *vwork_n;   /* work vectors of length n */
 13:   Vec       se;         /* Optional standard error vector */
 14:   PetscBool se_flg;     /* flag for -ksp_lsqr_set_standard_error */
 15:   PetscReal arnorm;     /* Norm of the vector A.r */
 16:   PetscReal anorm;      /* Frobenius norm of the matrix A */
 17:   PetscReal rhs_norm;   /* Norm of the right hand side */
 18: } KSP_LSQR;

 20: static PetscErrorCode  VecSquare(Vec v)
 21: {
 23:   PetscScalar    *x;
 24:   PetscInt       i, n;

 27:   VecGetLocalSize(v, &n);
 28:   VecGetArray(v, &x);
 29:   for (i = 0; i < n; i++) x[i] *= PetscConj(x[i]);
 30:   VecRestoreArray(v, &x);
 31:   return(0);
 32: }

 34: static PetscErrorCode KSPSetUp_LSQR(KSP ksp)
 35: {
 37:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;
 38:   PetscBool      nopreconditioner;

 41:   PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);
 42:   /*  nopreconditioner =PETSC_FALSE; */

 44:   lsqr->nwork_m = 2;
 45:   if (lsqr->vwork_m) {
 46:     VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
 47:   }
 48:   if (nopreconditioner) lsqr->nwork_n = 4;
 49:   else lsqr->nwork_n = 5;

 51:   if (lsqr->vwork_n) {
 52:     VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
 53:   }
 54:   KSPCreateVecs(ksp,lsqr->nwork_n,&lsqr->vwork_n,lsqr->nwork_m,&lsqr->vwork_m);
 55:   if (lsqr->se_flg && !lsqr->se) {
 56:     /* lsqr->se is not set by user, get it from pmat */
 57:     Vec *se;
 58:     KSPCreateVecs(ksp,1,&se,0,NULL);
 59:     lsqr->se = *se;
 60:     PetscFree(se);
 61:   }
 62:   return(0);
 63: }

 65: static PetscErrorCode KSPSolve_LSQR(KSP ksp)
 66: {
 68:   PetscInt       i,size1,size2;
 69:   PetscScalar    rho,rhobar,phi,phibar,theta,c,s,tmp,tau;
 70:   PetscReal      beta,alpha,rnorm;
 71:   Vec            X,B,V,V1,U,U1,TMP,W,W2,SE,Z = NULL;
 72:   Mat            Amat,Pmat;
 73:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;
 74:   PetscBool      diagonalscale,nopreconditioner;

 77:   PCGetDiagonalScale(ksp->pc,&diagonalscale);
 78:   if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);

 80:   PCGetOperators(ksp->pc,&Amat,&Pmat);
 81:   PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);

 83:   /*  nopreconditioner =PETSC_FALSE; */
 84:   /* Calculate norm of right hand side */
 85:   VecNorm(ksp->vec_rhs,NORM_2,&lsqr->rhs_norm);

 87:   /* mark norm of matrix with negative number to indicate it has not yet been computed */
 88:   lsqr->anorm = -1.0;

 90:   /* vectors of length m, where system size is mxn */
 91:   B  = ksp->vec_rhs;
 92:   U  = lsqr->vwork_m[0];
 93:   U1 = lsqr->vwork_m[1];

 95:   /* vectors of length n */
 96:   X  = ksp->vec_sol;
 97:   W  = lsqr->vwork_n[0];
 98:   V  = lsqr->vwork_n[1];
 99:   V1 = lsqr->vwork_n[2];
100:   W2 = lsqr->vwork_n[3];
101:   if (!nopreconditioner) Z = lsqr->vwork_n[4];

103:   /* standard error vector */
104:   SE = lsqr->se;
105:   if (SE) {
106:     VecGetSize(SE,&size1);
107:     VecGetSize(X,&size2);
108:     if (size1 != size2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Standard error vector (size %D) does not match solution vector (size %D)",size1,size2);
109:     VecSet(SE,0.0);
110:   }

112:   /* Compute initial residual, temporarily use work vector u */
113:   if (!ksp->guess_zero) {
114:     KSP_MatMult(ksp,Amat,X,U);       /*   u <- b - Ax     */
115:     VecAYPX(U,-1.0,B);
116:   } else {
117:     VecCopy(B,U);            /*   u <- b (x is 0) */
118:   }

120:   /* Test for nothing to do */
121:   VecNorm(U,NORM_2,&rnorm);
122:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
123:   ksp->its   = 0;
124:   ksp->rnorm = rnorm;
125:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
126:   KSPLogResidualHistory(ksp,rnorm);
127:   KSPMonitor(ksp,0,rnorm);
128:   (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
129:   if (ksp->reason) return(0);

131:   beta = rnorm;
132:   VecScale(U,1.0/beta);
133:   KSP_MatMultTranspose(ksp,Amat,U,V);
134:   if (nopreconditioner) {
135:     VecNorm(V,NORM_2,&alpha);
136:   } else {
137:     /* this is an application of the preconditioner for the normal equations; not the operator, see the manual page */
138:     PCApply(ksp->pc,V,Z);
139:     VecDotRealPart(V,Z,&alpha);
140:     if (alpha <= 0.0) {
141:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
142:       return(0);
143:     }
144:     alpha = PetscSqrtReal(alpha);
145:     VecScale(Z,1.0/alpha);
146:   }
147:   VecScale(V,1.0/alpha);

149:   if (nopreconditioner) {
150:     VecCopy(V,W);
151:   } else {
152:     VecCopy(Z,W);
153:   }

155:   lsqr->arnorm = alpha * beta;
156:   phibar       = beta;
157:   rhobar       = alpha;
158:   i            = 0;
159:   do {
160:     if (nopreconditioner) {
161:       KSP_MatMult(ksp,Amat,V,U1);
162:     } else {
163:       KSP_MatMult(ksp,Amat,Z,U1);
164:     }
165:     VecAXPY(U1,-alpha,U);
166:     VecNorm(U1,NORM_2,&beta);
167:     if (beta > 0.0) {
168:       VecScale(U1,1.0/beta); /* beta*U1 = Amat*V - alpha*U */
169:     }

171:     KSP_MatMultTranspose(ksp,Amat,U1,V1);
172:     VecAXPY(V1,-beta,V);
173:     if (nopreconditioner) {
174:       VecNorm(V1,NORM_2,&alpha);
175:     } else {
176:       PCApply(ksp->pc,V1,Z);
177:       VecDotRealPart(V1,Z,&alpha);
178:       if (alpha <= 0.0) {
179:         ksp->reason = KSP_DIVERGED_BREAKDOWN;
180:         break;
181:       }
182:       alpha = PetscSqrtReal(alpha);
183:       VecScale(Z,1.0/alpha);
184:     }
185:     VecScale(V1,1.0/alpha); /* alpha*V1 = Amat^T*U1 - beta*V */
186:     rho    = PetscSqrtScalar(rhobar*rhobar + beta*beta);
187:     c      = rhobar / rho;
188:     s      = beta / rho;
189:     theta  = s * alpha;
190:     rhobar = -c * alpha;
191:     phi    = c * phibar;
192:     phibar = s * phibar;
193:     tau    = s * phi;

195:     VecAXPY(X,phi/rho,W);  /*    x <- x + (phi/rho) w   */

197:     if (SE) {
198:       VecCopy(W,W2);
199:       VecSquare(W2);
200:       VecScale(W2,1.0/(rho*rho));
201:       VecAXPY(SE, 1.0, W2); /* SE <- SE + (w^2/rho^2) */
202:     }
203:     if (nopreconditioner) {
204:       VecAYPX(W,-theta/rho,V1);  /* w <- v - (theta/rho) w */
205:     } else {
206:       VecAYPX(W,-theta/rho,Z);   /* w <- z - (theta/rho) w */
207:     }

209:     lsqr->arnorm = alpha*PetscAbsScalar(tau);
210:     rnorm        = PetscRealPart(phibar);

212:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
213:     ksp->its++;
214:     ksp->rnorm = rnorm;
215:     PetscObjectSAWsGrantAccess((PetscObject)ksp);
216:     KSPLogResidualHistory(ksp,rnorm);
217:     KSPMonitor(ksp,i+1,rnorm);
218:     (*ksp->converged)(ksp,i+1,rnorm,&ksp->reason,ksp->cnvP);
219:     if (ksp->reason) break;
220:     SWAP(U1,U,TMP);
221:     SWAP(V1,V,TMP);

223:     i++;
224:   } while (i<ksp->max_it);
225:   if (i >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;

227:   /* Finish off the standard error estimates */
228:   if (SE) {
229:     tmp  = 1.0;
230:     MatGetSize(Amat,&size1,&size2);
231:     if (size1 > size2) tmp = size1 - size2;
232:     tmp  = rnorm / PetscSqrtScalar(tmp);
233:     VecSqrtAbs(SE);
234:     VecScale(SE,tmp);
235:   }
236:   return(0);
237: }


240: PetscErrorCode KSPDestroy_LSQR(KSP ksp)
241: {
242:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;

246:   /* Free work vectors */
247:   if (lsqr->vwork_n) {
248:     VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
249:   }
250:   if (lsqr->vwork_m) {
251:     VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
252:   }
253:   if (lsqr->se_flg) {
254:     VecDestroy(&lsqr->se);
255:   }
256:   PetscFree(ksp->data);
257:   return(0);
258: }

260: PetscErrorCode  KSPLSQRSetStandardErrorVec(KSP ksp, Vec se)
261: {
262:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;

266:   VecDestroy(&lsqr->se);
267:   lsqr->se = se;
268:   return(0);
269: }

271: PetscErrorCode  KSPLSQRGetStandardErrorVec(KSP ksp,Vec *se)
272: {
273:   KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;

276:   *se = lsqr->se;
277:   return(0);
278: }

280: PetscErrorCode  KSPLSQRGetArnorm(KSP ksp,PetscReal *arnorm, PetscReal *rhs_norm, PetscReal *anorm)
281: {
282:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;

286:   *arnorm = lsqr->arnorm;
287:   if (anorm) {
288:     if (lsqr->anorm < 0.0) {
289:       PC  pc;
290:       Mat Amat;
291:       KSPGetPC(ksp,&pc);
292:       PCGetOperators(pc,&Amat,NULL);
293:       MatNorm(Amat,NORM_FROBENIUS,&lsqr->anorm);
294:     }
295:     *anorm = lsqr->anorm;
296:   }
297:   if (rhs_norm) *rhs_norm = lsqr->rhs_norm;
298:   return(0);
299: }

301: /*@C
302:    KSPLSQRMonitorDefault - Print the residual norm at each iteration of the LSQR method and the norm of the residual of the normal equations A'*A x = A' b

304:    Collective on KSP

306:    Input Parameters:
307: +  ksp   - iterative context
308: .  n     - iteration number
309: .  rnorm - 2-norm (preconditioned) residual value (may be estimated).
310: -  dummy - viewer and format context

312:    Level: intermediate

314: .keywords: KSP, default, monitor, residual

316: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate(), KSPMonitorDefault()
317: @*/
318: PetscErrorCode  KSPLSQRMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
319: {
321:   PetscViewer    viewer = dummy->viewer;
322:   KSP_LSQR       *lsqr  = (KSP_LSQR*)ksp->data;

325:   PetscViewerPushFormat(viewer,dummy->format);
326:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
327:   if (((PetscObject)ksp)->prefix) {
328:     PetscViewerASCIIPrintf(viewer,"  Residual norm and norm of normal equations for %s solve.\n",((PetscObject)ksp)->prefix);
329:   }
330:   if (!n) {
331:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e\n",n,rnorm);
332:   } else {
333:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e Residual norm normal equations %14.12e\n",n,rnorm,lsqr->arnorm);
334:   }
335:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
336:   PetscViewerPopFormat(viewer);
337:   return(0);
338: }

340: PetscErrorCode KSPSetFromOptions_LSQR(PetscOptionItems *PetscOptionsObject,KSP ksp)
341: {
343:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;

346:   PetscOptionsHead(PetscOptionsObject,"KSP LSQR Options");
347:   PetscOptionsName("-ksp_lsqr_set_standard_error","Set Standard Error Estimates of Solution","KSPLSQRSetStandardErrorVec",&lsqr->se_flg);
348:   KSPMonitorSetFromOptions(ksp,"-ksp_lsqr_monitor","Monitor residual norm and norm of residual of normal equations","KSPMonitorSet",KSPLSQRMonitorDefault);
349:   PetscOptionsTail();
350:   return(0);
351: }

353: PetscErrorCode KSPView_LSQR(KSP ksp,PetscViewer viewer)
354: {
355:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;
357:   PetscBool      iascii;

360:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
361:   if (iascii) {
362:     if (lsqr->se) {
363:       PetscReal rnorm;
364:       KSPLSQRGetStandardErrorVec(ksp,&lsqr->se);
365:       VecNorm(lsqr->se,NORM_2,&rnorm);
366:       PetscViewerASCIIPrintf(viewer,"  Norm of Standard Error %g, Iterations %D\n",(double)rnorm,ksp->its);
367:     }
368:   }
369:   return(0);
370: }

372: /*@C
373:    KSPLSQRDefaultConverged - Determines convergence of the LSQR Krylov method. This calls KSPConvergedDefault() and if that does not determine convergence then checks
374:       convergence for the least squares problem.

376:    Collective on KSP

378:    Input Parameters:
379: +  ksp   - iterative context
380: .  n     - iteration number
381: .  rnorm - 2-norm residual value (may be estimated)
382: -  ctx - convergence context which must be created by KSPConvergedDefaultCreate()

384:    reason is set to:
385: +   positive - if the iteration has converged;
386: .   negative - if residual norm exceeds divergence threshold;
387: -   0 - otherwise.

389:    Notes:
390:       Possible convergence for the least squares problem (which is based on the residual of the normal equations) are KSP_CONVERGED_RTOL_NORMAL norm and KSP_CONVERGED_ATOL_NORMAL.

392:    Level: intermediate

394: .keywords: KSP, default, convergence, residual

396: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
397:           KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy(), KSPConvergedDefault()
398: @*/
399: PetscErrorCode  KSPLSQRDefaultConverged(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
400: {
402:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;

405:   KSPConvergedDefault(ksp,n,rnorm,reason,ctx);
406:   if (!n || *reason) return(0);
407:   if (lsqr->arnorm/lsqr->rhs_norm < ksp->rtol) *reason = KSP_CONVERGED_RTOL_NORMAL;
408:   if (lsqr->arnorm < ksp->abstol) *reason = KSP_CONVERGED_ATOL_NORMAL;
409:   return(0);
410: }

412: /*MC
413:      KSPLSQR - This implements LSQR

415:    Options Database Keys:
416: +   -ksp_lsqr_set_standard_error  - Set Standard Error Estimates of Solution see KSPLSQRSetStandardErrorVec()
417: .   -ksp_lsqr_monitor - Monitor residual norm and norm of residual of normal equations
418: -   see KSPSolve()

420:    Level: beginner

422:    Notes:
423:      Supports non-square (rectangular) matrices.

425:      This varient, when applied with no preconditioning is identical to the original algorithm in exact arithematic; however, in practice, with no preconditioning
426:      due to inexact arithematic, it can converge differently. Hence when no preconditioner is used (PCType PCNONE) it automatically reverts to the original algorithm.

428:      With the PETSc built-in preconditioners, such as ICC, one should call KSPSetOperators(ksp,A,A'*A)) since the preconditioner needs to work
429:      for the normal equations A'*A.

431:      Supports only left preconditioning.

433:    References:
434: .  1. - The original unpreconditioned algorithm can be found in Paige and Saunders, ACM Transactions on Mathematical Software, Vol 8, 1982.

436:      In exact arithmetic the LSQR method (with no preconditioning) is identical to the KSPCG algorithm applied to the normal equations.
437:      The preconditioned varient was implemented by Bas van't Hof and is essentially a left preconditioning for the Normal Equations. It appears the implementation with preconditioner
438:      track the true norm of the residual and uses that in the convergence test.

440:    Developer Notes:
441:     How is this related to the KSPCGNE implementation? One difference is that KSPCGNE applies
442:             the preconditioner transpose times the preconditioner,  so one does not need to pass A'*A as the third argument to KSPSetOperators().


445:    For least squares problems without a zero to A*x = b, there are additional convergence tests for the residual of the normal equations, A'*(b - Ax), see KSPLSQRDefaultConverged()

447: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPLSQRDefaultConverged()

449: M*/
450: PETSC_EXTERN PetscErrorCode KSPCreate_LSQR(KSP ksp)
451: {
452:   KSP_LSQR       *lsqr;
453:   void           *ctx;

457:   PetscNewLog(ksp,&lsqr);
458:   lsqr->se     = NULL;
459:   lsqr->se_flg = PETSC_FALSE;
460:   lsqr->arnorm = 0.0;
461:   ksp->data    = (void*)lsqr;
462:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);

464:   ksp->ops->setup          = KSPSetUp_LSQR;
465:   ksp->ops->solve          = KSPSolve_LSQR;
466:   ksp->ops->destroy        = KSPDestroy_LSQR;
467:   ksp->ops->setfromoptions = KSPSetFromOptions_LSQR;
468:   ksp->ops->view           = KSPView_LSQR;

470:   KSPConvergedDefaultCreate(&ctx);
471:   KSPSetConvergenceTest(ksp,KSPLSQRDefaultConverged,ctx,KSPConvergedDefaultDestroy);
472:   return(0);
473: }