Actual source code: lsqr.c

petsc-3.5.2 2014-09-08
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>
  8: #include <../src/ksp/ksp/impls/lsqr/lsqr.h>

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

 21: extern PetscErrorCode  VecSquare(Vec);

 25: static PetscErrorCode KSPSetUp_LSQR(KSP ksp)
 26: {
 28:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;
 29:   PetscBool      nopreconditioner;

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

 35:   lsqr->nwork_m = 2;
 36:   if (lsqr->vwork_m) {
 37:     VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
 38:   }
 39:   if (nopreconditioner) lsqr->nwork_n = 4;
 40:   else lsqr->nwork_n = 5;

 42:   if (lsqr->vwork_n) {
 43:     VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
 44:   }
 45:   KSPGetVecs(ksp,lsqr->nwork_n,&lsqr->vwork_n,lsqr->nwork_m,&lsqr->vwork_m);
 46:   if (lsqr->se_flg && !lsqr->se) {
 47:     /* lsqr->se is not set by user, get it from pmat */
 48:     Vec *se;
 49:     KSPGetVecs(ksp,1,&se,0,NULL);
 50:     lsqr->se = *se;
 51:     PetscFree(se);
 52:   }
 53:   return(0);
 54: }

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

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

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

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

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

 83:   /* vectors of length m, where system size is mxn */
 84:   B  = ksp->vec_rhs;
 85:   U  = lsqr->vwork_m[0];
 86:   U1 = lsqr->vwork_m[1];

 88:   /* vectors of length n */
 89:   X  = ksp->vec_sol;
 90:   W  = lsqr->vwork_n[0];
 91:   V  = lsqr->vwork_n[1];
 92:   V1 = lsqr->vwork_n[2];
 93:   W2 = lsqr->vwork_n[3];
 94:   if (!nopreconditioner) Z = lsqr->vwork_n[4];

 96:   /* standard error vector */
 97:   SE = lsqr->se;
 98:   if (SE) {
 99:     VecGetSize(SE,&size1);
100:     VecGetSize(X,&size2);
101:     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);
102:     VecSet(SE,0.0);
103:   }

105:   /* Compute initial residual, temporarily use work vector u */
106:   if (!ksp->guess_zero) {
107:     KSP_MatMult(ksp,Amat,X,U);       /*   u <- b - Ax     */
108:     VecAYPX(U,-1.0,B);
109:   } else {
110:     VecCopy(B,U);            /*   u <- b (x is 0) */
111:   }

113:   /* Test for nothing to do */
114:   VecNorm(U,NORM_2,&rnorm);
115:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
116:   ksp->its   = 0;
117:   ksp->rnorm = rnorm;
118:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
119:   KSPLogResidualHistory(ksp,rnorm);
120:   KSPMonitor(ksp,0,rnorm);
121:   (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
122:   if (ksp->reason) return(0);

124:   beta = rnorm;
125:   VecScale(U,1.0/beta);
126:   KSP_MatMultTranspose(ksp,Amat,U,V);
127:   if (nopreconditioner) {
128:     VecNorm(V,NORM_2,&alpha);
129:   } else {
130:     PCApply(ksp->pc,V,Z);
131:     VecDotRealPart(V,Z,&alpha);
132:     if (alpha <= 0.0) {
133:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
134:       return(0);
135:     }
136:     alpha = PetscSqrtReal(alpha);
137:     VecScale(Z,1.0/alpha);
138:   }
139:   VecScale(V,1.0/alpha);

141:   if (nopreconditioner) {
142:     VecCopy(V,W);
143:   } else {
144:     VecCopy(Z,W);
145:   }

147:   lsqr->arnorm = alpha * beta;
148:   phibar       = beta;
149:   rhobar       = alpha;
150:   i            = 0;
151:   do {
152:     if (nopreconditioner) {
153:       KSP_MatMult(ksp,Amat,V,U1);
154:     } else {
155:       KSP_MatMult(ksp,Amat,Z,U1);
156:     }
157:     VecAXPY(U1,-alpha,U);
158:     VecNorm(U1,NORM_2,&beta);
159:     if (beta == 0.0) {
160:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
161:       break;
162:     }
163:     VecScale(U1,1.0/beta);

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

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

191:     if (SE) {
192:       VecCopy(W,W2);
193:       VecSquare(W2);
194:       VecScale(W2,1.0/(rho*rho));
195:       VecAXPY(SE, 1.0, W2); /* SE <- SE + (w^2/rho^2) */
196:     }
197:     if (nopreconditioner) {
198:       VecAYPX(W,-theta/rho,V1);  /* w <- v - (theta/rho) w */
199:     } else {
200:       VecAYPX(W,-theta/rho,Z);   /* w <- z - (theta/rho) w */
201:     }

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

206:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
207:     ksp->its++;
208:     ksp->rnorm = rnorm;
209:     PetscObjectSAWsGrantAccess((PetscObject)ksp);
210:     KSPLogResidualHistory(ksp,rnorm);
211:     KSPMonitor(ksp,i+1,rnorm);
212:     (*ksp->converged)(ksp,i+1,rnorm,&ksp->reason,ksp->cnvP);
213:     if (ksp->reason) break;
214:     SWAP(U1,U,TMP);
215:     SWAP(V1,V,TMP);

217:     i++;
218:   } while (i<ksp->max_it);
219:   if (i >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;

221:   /* Finish off the standard error estimates */
222:   if (SE) {
223:     tmp  = 1.0;
224:     MatGetSize(Amat,&size1,&size2);
225:     if (size1 > size2) tmp = size1 - size2;
226:     tmp  = rnorm / PetscSqrtScalar(tmp);
227:     VecSqrtAbs(SE);
228:     VecScale(SE,tmp);
229:   }
230:   return(0);
231: }


236: PetscErrorCode KSPDestroy_LSQR(KSP ksp)
237: {
238:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;

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

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

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

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

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

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

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

305: /*@C
306:    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

308:    Collective on KSP

310:    Input Parameters:
311: +  ksp   - iterative context
312: .  n     - iteration number
313: .  rnorm - 2-norm (preconditioned) residual value (may be estimated).
314: -  dummy - unused monitor context

316:    Level: intermediate

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

320: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate(), KSPMonitorDefault()
321: @*/
322: PetscErrorCode  KSPLSQRMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
323: {
325:   PetscViewer    viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp));
326:   KSP_LSQR       *lsqr  = (KSP_LSQR*)ksp->data;

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

344: PetscErrorCode KSPSetFromOptions_LSQR(KSP ksp)
345: {
347:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;
348:   char           monfilename[PETSC_MAX_PATH_LEN];
349:   PetscViewer    monviewer;
350:   PetscBool      flg;

353:   PetscOptionsHead("KSP LSQR Options");
354:   PetscOptionsName("-ksp_lsqr_set_standard_error","Set Standard Error Estimates of Solution","KSPLSQRSetStandardErrorVec",&lsqr->se_flg);
355:   PetscOptionsString("-ksp_lsqr_monitor","Monitor residual norm and norm of residual of normal equations","KSPMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
356:   if (flg) {
357:     PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ksp),monfilename,&monviewer);
358:     KSPMonitorSet(ksp,KSPLSQRMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
359:   }
360:   PetscOptionsTail();
361:   return(0);
362: }

366: PetscErrorCode KSPView_LSQR(KSP ksp,PetscViewer viewer)
367: {
368:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;
370:   PetscBool      iascii;

373:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
374:   if (iascii) {
375:     if (lsqr->se) {
376:       PetscReal rnorm;
377:       KSPLSQRGetStandardErrorVec(ksp,&lsqr->se);
378:       VecNorm(lsqr->se,NORM_2,&rnorm);
379:       PetscViewerASCIIPrintf(viewer,"  Norm of Standard Error %g, Iterations %D\n",(double)rnorm,ksp->its);
380:     }
381:   }
382:   return(0);
383: }

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

391:    Collective on KSP

393:    Input Parameters:
394: +  ksp   - iterative context
395: .  n     - iteration number
396: .  rnorm - 2-norm residual value (may be estimated)
397: -  ctx - convergence context which must be created by KSPConvergedDefaultCreate()

399:    reason is set to:
400: +   positive - if the iteration has converged;
401: .   negative - if residual norm exceeds divergence threshold;
402: -   0 - otherwise.

404:    Notes:
405:       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.

407:    Level: intermediate

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

411: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
412:           KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy(), KSPConvergedDefault()
413: @*/
414: PetscErrorCode  KSPLSQRDefaultConverged(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
415: {
417:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;

420:   KSPConvergedDefault(ksp,n,rnorm,reason,ctx);
421:   if (!n || *reason) return(0);
422:   if (lsqr->arnorm/lsqr->rhs_norm < ksp->rtol) *reason = KSP_CONVERGED_RTOL_NORMAL;
423:   if (lsqr->arnorm < ksp->abstol) *reason = KSP_CONVERGED_ATOL_NORMAL;
424:   return(0);
425: }



429: /*MC
430:      KSPLSQR - This implements LSQR

432:    Options Database Keys:
433: +   -ksp_lsqr_set_standard_error  - Set Standard Error Estimates of Solution see KSPLSQRSetStandardErrorVec()
434: .   -ksp_lsqr_monitor - Monitor residual norm and norm of residual of normal equations
435: -   see KSPSolve()

437:    Level: beginner

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

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

446:      Supports only left preconditioning.

448:    References:The original unpreconditioned algorithm can be found in Paige and Saunders, ACM Transactions on Mathematical Software, Vol 8, pp 43-71, 1982.
449:      In exact arithmetic the LSQR method (with no preconditioning) is identical to the KSPCG algorithm applied to the normal equations.
450:      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
451:      track the true norm of the residual and uses that in the convergence test.

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


457:    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()

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

461: M*/
464: PETSC_EXTERN PetscErrorCode KSPCreate_LSQR(KSP ksp)
465: {
466:   KSP_LSQR       *lsqr;

470:   PetscNewLog(ksp,&lsqr);
471:   lsqr->se     = NULL;
472:   lsqr->se_flg = PETSC_FALSE;
473:   lsqr->arnorm = 0.0;
474:   ksp->data    = (void*)lsqr;
475:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);

477:   ksp->ops->setup          = KSPSetUp_LSQR;
478:   ksp->ops->solve          = KSPSolve_LSQR;
479:   ksp->ops->destroy        = KSPDestroy_LSQR;
480:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
481:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
482:   ksp->ops->setfromoptions = KSPSetFromOptions_LSQR;
483:   ksp->ops->view           = KSPView_LSQR;
484:   ksp->converged           = KSPLSQRDefaultConverged;
485:   return(0);
486: }

490: PetscErrorCode  VecSquare(Vec v)
491: {
493:   PetscScalar    *x;
494:   PetscInt       i, n;

497:   VecGetLocalSize(v, &n);
498:   VecGetArray(v, &x);
499:   for (i = 0; i < n; i++) x[i] *= x[i];
500:   VecRestoreArray(v, &x);
501:   return(0);
502: }