Actual source code: trlanczos.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    SLEPc singular value solver: "trlanczos"

 13:    Method: Thick-restart Lanczos

 15:    Algorithm:

 17:        Golub-Kahan-Lanczos bidiagonalization with thick-restart.

 19:    References:

 21:        [1] G.H. Golub and W. Kahan, "Calculating the singular values
 22:            and pseudo-inverse of a matrix", SIAM J. Numer. Anal. Ser.
 23:            B 2:205-224, 1965.

 25:        [2] V. Hernandez, J.E. Roman, and A. Tomas, "A robust and
 26:            efficient parallel SVD solver based on restarted Lanczos
 27:            bidiagonalization", Elec. Trans. Numer. Anal. 31:68-85,
 28:            2008.
 29: */

 31: #include <slepc/private/svdimpl.h>          /*I "slepcsvd.h" I*/

 33: static PetscBool  cited = PETSC_FALSE;
 34: static const char citation[] =
 35:   "@Article{slepc-svd,\n"
 36:   "   author = \"V. Hern{\\'a}ndez and J. E. Rom{\\'a}n and A. Tom{\\'a}s\",\n"
 37:   "   title = \"A robust and efficient parallel {SVD} solver based on restarted {Lanczos} bidiagonalization\",\n"
 38:   "   journal = \"Electron. Trans. Numer. Anal.\",\n"
 39:   "   volume = \"31\",\n"
 40:   "   pages = \"68--85\",\n"
 41:   "   year = \"2008\"\n"
 42:   "}\n";

 44: typedef struct {
 45:   PetscBool oneside;
 46: } SVD_TRLANCZOS;

 48: PetscErrorCode SVDSetUp_TRLanczos(SVD svd)
 49: {
 51:   PetscInt       N;

 54:   SVDMatGetSize(svd,NULL,&N);
 55:   SVDSetDimensions_Default(svd);
 56:   if (svd->ncv>svd->nsv+svd->mpd) SETERRQ(PetscObjectComm((PetscObject)svd),1,"The value of ncv must not be larger than nev+mpd");
 57:   if (!svd->max_it) svd->max_it = PetscMax(N/svd->ncv,100);
 58:   svd->leftbasis = PETSC_TRUE;
 59:   SVDAllocateSolution(svd,1);
 60:   DSSetType(svd->ds,DSSVD);
 61:   DSSetCompact(svd->ds,PETSC_TRUE);
 62:   DSAllocate(svd->ds,svd->ncv);
 63:   return(0);
 64: }

 66: static PetscErrorCode SVDOneSideTRLanczosMGS(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,BV U,PetscInt nconv,PetscInt l,PetscInt n,PetscScalar* work)
 67: {
 69:   PetscReal      a,b;
 70:   PetscInt       i,k=nconv+l;
 71:   Vec            ui,ui1,vi;

 74:   BVGetColumn(V,k,&vi);
 75:   BVGetColumn(U,k,&ui);
 76:   SVDMatMult(svd,PETSC_FALSE,vi,ui);
 77:   BVRestoreColumn(V,k,&vi);
 78:   BVRestoreColumn(U,k,&ui);
 79:   if (l>0) {
 80:     for (i=0;i<l;i++) work[i]=beta[i+nconv];
 81:     BVMultColumn(U,-1.0,1.0,k,work);
 82:   }
 83:   BVNormColumn(U,k,NORM_2,&a);
 84:   BVScaleColumn(U,k,1.0/a);
 85:   alpha[k] = a;

 87:   for (i=k+1;i<n;i++) {
 88:     BVGetColumn(V,i,&vi);
 89:     BVGetColumn(U,i-1,&ui1);
 90:     SVDMatMult(svd,PETSC_TRUE,ui1,vi);
 91:     BVRestoreColumn(V,i,&vi);
 92:     BVRestoreColumn(U,i-1,&ui1);
 93:     BVOrthogonalizeColumn(V,i,NULL,&b,NULL);
 94:     BVScaleColumn(V,i,1.0/b);
 95:     beta[i-1] = b;

 97:     BVGetColumn(V,i,&vi);
 98:     BVGetColumn(U,i,&ui);
 99:     SVDMatMult(svd,PETSC_FALSE,vi,ui);
100:     BVRestoreColumn(V,i,&vi);
101:     BVGetColumn(U,i-1,&ui1);
102:     VecAXPY(ui,-b,ui1);
103:     BVRestoreColumn(U,i-1,&ui1);
104:     BVRestoreColumn(U,i,&ui);
105:     BVNormColumn(U,i,NORM_2,&a);
106:     BVScaleColumn(U,i,1.0/a);
107:     alpha[i] = a;
108:   }

110:   BVGetColumn(V,n,&vi);
111:   BVGetColumn(U,n-1,&ui1);
112:   SVDMatMult(svd,PETSC_TRUE,ui1,vi);
113:   BVRestoreColumn(V,n,&vi);
114:   BVRestoreColumn(U,n-1,&ui1);
115:   BVOrthogonalizeColumn(V,n,NULL,&b,NULL);
116:   beta[n-1] = b;
117:   return(0);
118: }

120: /*
121:   Custom CGS orthogonalization, preprocess after first orthogonalization
122: */
123: static PetscErrorCode SVDOrthogonalizeCGS(BV V,PetscInt i,PetscScalar* h,PetscReal a,BVOrthogRefineType refine,PetscReal eta,PetscReal *norm)
124: {
126:   PetscReal      sum,onorm;
127:   PetscScalar    dot;
128:   PetscInt       j;

131:   switch (refine) {
132:   case BV_ORTHOG_REFINE_NEVER:
133:     BVNormColumn(V,i,NORM_2,norm);
134:     break;
135:   case BV_ORTHOG_REFINE_ALWAYS:
136:     BVSetActiveColumns(V,0,i);
137:     BVDotColumn(V,i,h);
138:     BVMultColumn(V,-1.0,1.0,i,h);
139:     BVNormColumn(V,i,NORM_2,norm);
140:     break;
141:   case BV_ORTHOG_REFINE_IFNEEDED:
142:     dot = h[i];
143:     onorm = PetscSqrtReal(PetscRealPart(dot)) / a;
144:     sum = 0.0;
145:     for (j=0;j<i;j++) {
146:       sum += PetscRealPart(h[j] * PetscConj(h[j]));
147:     }
148:     *norm = PetscRealPart(dot)/(a*a) - sum;
149:     if (*norm>0.0) *norm = PetscSqrtReal(*norm);
150:     else {
151:       BVNormColumn(V,i,NORM_2,norm);
152:     }
153:     if (*norm < eta*onorm) {
154:       BVSetActiveColumns(V,0,i);
155:       BVDotColumn(V,i,h);
156:       BVMultColumn(V,-1.0,1.0,i,h);
157:       BVNormColumn(V,i,NORM_2,norm);
158:     }
159:     break;
160:   }
161:   return(0);
162: }

164: static PetscErrorCode SVDOneSideTRLanczosCGS(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,BV U,PetscInt nconv,PetscInt l,PetscInt n,PetscScalar* work)
165: {
166:   PetscErrorCode     ierr;
167:   PetscReal          a,b,eta;
168:   PetscInt           i,j,k=nconv+l;
169:   Vec                ui,ui1,vi;
170:   BVOrthogRefineType refine;

173:   BVGetColumn(V,k,&vi);
174:   BVGetColumn(U,k,&ui);
175:   SVDMatMult(svd,PETSC_FALSE,vi,ui);
176:   BVRestoreColumn(V,k,&vi);
177:   BVRestoreColumn(U,k,&ui);
178:   if (l>0) {
179:     for (i=0;i<l;i++) work[i]=beta[i+nconv];
180:     BVMultColumn(U,-1.0,1.0,k,work);
181:   }
182:   BVGetOrthogonalization(V,NULL,&refine,&eta,NULL);

184:   for (i=k+1;i<n;i++) {
185:     BVGetColumn(V,i,&vi);
186:     BVGetColumn(U,i-1,&ui1);
187:     SVDMatMult(svd,PETSC_TRUE,ui1,vi);
188:     BVRestoreColumn(V,i,&vi);
189:     BVRestoreColumn(U,i-1,&ui1);
190:     BVNormColumnBegin(U,i-1,NORM_2,&a);
191:     if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
192:       BVSetActiveColumns(V,0,i+1);
193:       BVGetColumn(V,i,&vi);
194:       BVDotVecBegin(V,vi,work);
195:     } else {
196:       BVSetActiveColumns(V,0,i);
197:       BVDotColumnBegin(V,i,work);
198:     }
199:     BVNormColumnEnd(U,i-1,NORM_2,&a);
200:     if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
201:       BVDotVecEnd(V,vi,work);
202:       BVRestoreColumn(V,i,&vi);
203:       BVSetActiveColumns(V,0,i);
204:     } else {
205:       BVDotColumnEnd(V,i,work);
206:     }

208:     BVScaleColumn(U,i-1,1.0/a);
209:     for (j=0;j<i;j++) work[j] = work[j] / a;
210:     BVMultColumn(V,-1.0,1.0/a,i,work);
211:     SVDOrthogonalizeCGS(V,i,work,a,refine,eta,&b);
212:     BVScaleColumn(V,i,1.0/b);
213:     if (PetscAbsReal(b)<10*PETSC_MACHINE_EPSILON) SETERRQ(PetscObjectComm((PetscObject)svd),1,"Recurrence generated a zero vector; use a two-sided variant");

215:     BVGetColumn(V,i,&vi);
216:     BVGetColumn(U,i,&ui);
217:     BVGetColumn(U,i-1,&ui1);
218:     SVDMatMult(svd,PETSC_FALSE,vi,ui);
219:     VecAXPY(ui,-b,ui1);
220:     BVRestoreColumn(V,i,&vi);
221:     BVRestoreColumn(U,i,&ui);
222:     BVRestoreColumn(U,i-1,&ui1);

224:     alpha[i-1] = a;
225:     beta[i-1] = b;
226:   }

228:   BVGetColumn(V,n,&vi);
229:   BVGetColumn(U,n-1,&ui1);
230:   SVDMatMult(svd,PETSC_TRUE,ui1,vi);
231:   BVRestoreColumn(V,n,&vi);
232:   BVRestoreColumn(U,n-1,&ui1);

234:   BVNormColumnBegin(svd->U,n-1,NORM_2,&a);
235:   if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
236:     BVSetActiveColumns(V,0,n+1);
237:     BVGetColumn(V,n,&vi);
238:     BVDotVecBegin(V,vi,work);
239:   } else {
240:     BVSetActiveColumns(V,0,n);
241:     BVDotColumnBegin(V,n,work);
242:   }
243:   BVNormColumnEnd(svd->U,n-1,NORM_2,&a);
244:   if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
245:     BVDotVecEnd(V,vi,work);
246:     BVRestoreColumn(V,n,&vi);
247:   } else {
248:     BVDotColumnEnd(V,n,work);
249:   }

251:   BVScaleColumn(U,n-1,1.0/a);
252:   for (j=0;j<n;j++) work[j] = work[j] / a;
253:   BVMultColumn(V,-1.0,1.0/a,n,work);
254:   SVDOrthogonalizeCGS(V,n,work,a,refine,eta,&b);
255:   BVSetActiveColumns(V,nconv,n);
256:   alpha[n-1] = a;
257:   beta[n-1] = b;
258:   return(0);
259: }

261: PetscErrorCode SVDSolve_TRLanczos(SVD svd)
262: {
264:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
265:   PetscReal      *alpha,*beta,lastbeta,resnorm;
266:   PetscScalar    *Q,*swork=NULL,*w;
267:   PetscInt       i,k,l,nv,ld;
268:   Mat            U,VT;
269:   PetscBool      conv;
270:   BVOrthogType   orthog;

273:   PetscCitationsRegister(citation,&cited);
274:   /* allocate working space */
275:   DSGetLeadingDimension(svd->ds,&ld);
276:   BVGetOrthogonalization(svd->V,&orthog,NULL,NULL,NULL);
277:   PetscMalloc1(ld,&w);
278:   if (lanczos->oneside) {
279:     PetscMalloc1(svd->ncv+1,&swork);
280:   }

282:   /* normalize start vector */
283:   if (!svd->nini) {
284:     BVSetRandomColumn(svd->V,0);
285:     BVOrthonormalizeColumn(svd->V,0,PETSC_TRUE,NULL,NULL);
286:   }

288:   l = 0;
289:   while (svd->reason == SVD_CONVERGED_ITERATING) {
290:     svd->its++;

292:     /* inner loop */
293:     nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
294:     BVSetActiveColumns(svd->V,svd->nconv,nv);
295:     BVSetActiveColumns(svd->U,svd->nconv,nv);
296:     DSGetArrayReal(svd->ds,DS_MAT_T,&alpha);
297:     beta = alpha + ld;
298:     if (lanczos->oneside) {
299:       if (orthog == BV_ORTHOG_MGS) {
300:         SVDOneSideTRLanczosMGS(svd,alpha,beta,svd->V,svd->U,svd->nconv,l,nv,swork);
301:       } else {
302:         SVDOneSideTRLanczosCGS(svd,alpha,beta,svd->V,svd->U,svd->nconv,l,nv,swork);
303:       }
304:     } else {
305:       SVDTwoSideLanczos(svd,alpha,beta,svd->V,svd->U,svd->nconv+l,nv);
306:     }
307:     lastbeta = beta[nv-1];
308:     DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha);
309:     BVScaleColumn(svd->V,nv,1.0/lastbeta);

311:     /* compute SVD of general matrix */
312:     DSSetDimensions(svd->ds,nv,nv,svd->nconv,svd->nconv+l);
313:     if (l==0) {
314:       DSSetState(svd->ds,DS_STATE_INTERMEDIATE);
315:     } else {
316:       DSSetState(svd->ds,DS_STATE_RAW);
317:     }
318:     DSSolve(svd->ds,w,NULL);
319:     DSSort(svd->ds,w,NULL,NULL,NULL,NULL);
320:     DSSynchronize(svd->ds,w,NULL);

322:     /* compute error estimates */
323:     k = 0;
324:     conv = PETSC_TRUE;
325:     DSGetArray(svd->ds,DS_MAT_U,&Q);
326:     DSGetArrayReal(svd->ds,DS_MAT_T,&alpha);
327:     beta = alpha + ld;
328:     for (i=svd->nconv;i<nv;i++) {
329:       svd->sigma[i] = PetscRealPart(w[i]);
330:       beta[i] = PetscRealPart(Q[nv-1+i*ld])*lastbeta;
331:       resnorm = PetscAbsReal(beta[i]);
332:       (*svd->converged)(svd,svd->sigma[i],resnorm,&svd->errest[i],svd->convergedctx);
333:       if (conv) {
334:         if (svd->errest[i] < svd->tol) k++;
335:         else conv = PETSC_FALSE;
336:       }
337:     }
338:     DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha);
339:     DSRestoreArray(svd->ds,DS_MAT_U,&Q);

341:     /* check convergence and update l */
342:     (*svd->stopping)(svd,svd->its,svd->max_it,svd->nconv+k,svd->nsv,&svd->reason,svd->stoppingctx);
343:     if (svd->reason != SVD_CONVERGED_ITERATING) l = 0;
344:     else l = PetscMax((nv-svd->nconv-k)/2,0);

346:     /* compute converged singular vectors and restart vectors */
347:     DSGetMat(svd->ds,DS_MAT_VT,&VT);
348:     BVMultInPlaceTranspose(svd->V,VT,svd->nconv,svd->nconv+k+l);
349:     MatDestroy(&VT);
350:     DSGetMat(svd->ds,DS_MAT_U,&U);
351:     BVMultInPlace(svd->U,U,svd->nconv,svd->nconv+k+l);
352:     MatDestroy(&U);

354:     /* copy the last vector to be the next initial vector */
355:     if (svd->reason == SVD_CONVERGED_ITERATING) {
356:       BVCopyColumn(svd->V,nv,svd->nconv+k+l);
357:     }

359:     svd->nconv += k;
360:     SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv);
361:   }

363:   /* orthonormalize U columns in one side method */
364:   if (lanczos->oneside) {
365:     for (i=0;i<svd->nconv;i++) {
366:       BVOrthonormalizeColumn(svd->U,i,PETSC_FALSE,NULL,NULL);
367:     }
368:   }

370:   /* free working space */
371:   PetscFree(w);
372:   if (swork) { PetscFree(swork); }
373:   return(0);
374: }

376: PetscErrorCode SVDSetFromOptions_TRLanczos(PetscOptionItems *PetscOptionsObject,SVD svd)
377: {
379:   PetscBool      set,val;
380:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;

383:   PetscOptionsHead(PetscOptionsObject,"SVD TRLanczos Options");

385:     PetscOptionsBool("-svd_trlanczos_oneside","Use one-side reorthogonalization","SVDTRLanczosSetOneSide",lanczos->oneside,&val,&set);
386:     if (set) { SVDTRLanczosSetOneSide(svd,val); }

388:   PetscOptionsTail();
389:   return(0);
390: }

392: static PetscErrorCode SVDTRLanczosSetOneSide_TRLanczos(SVD svd,PetscBool oneside)
393: {
394:   SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;

397:   if (lanczos->oneside != oneside) {
398:     lanczos->oneside = oneside;
399:     svd->state = SVD_STATE_INITIAL;
400:   }
401:   return(0);
402: }

404: /*@
405:    SVDTRLanczosSetOneSide - Indicate if the variant of the Lanczos method
406:    to be used is one-sided or two-sided.

408:    Logically Collective on SVD

410:    Input Parameters:
411: +  svd     - singular value solver
412: -  oneside - boolean flag indicating if the method is one-sided or not

414:    Options Database Key:
415: .  -svd_trlanczos_oneside <boolean> - Indicates the boolean flag

417:    Note:
418:    By default, a two-sided variant is selected, which is sometimes slightly
419:    more robust. However, the one-sided variant is faster because it avoids
420:    the orthogonalization associated to left singular vectors.

422:    Level: advanced

424: .seealso: SVDLanczosSetOneSide()
425: @*/
426: PetscErrorCode SVDTRLanczosSetOneSide(SVD svd,PetscBool oneside)
427: {

433:   PetscTryMethod(svd,"SVDTRLanczosSetOneSide_C",(SVD,PetscBool),(svd,oneside));
434:   return(0);
435: }

437: static PetscErrorCode SVDTRLanczosGetOneSide_TRLanczos(SVD svd,PetscBool *oneside)
438: {
439:   SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;

442:   *oneside = lanczos->oneside;
443:   return(0);
444: }

446: /*@
447:    SVDTRLanczosGetOneSide - Gets if the variant of the Lanczos method
448:    to be used is one-sided or two-sided.

450:    Not Collective

452:    Input Parameters:
453: .  svd     - singular value solver

455:    Output Parameters:
456: .  oneside - boolean flag indicating if the method is one-sided or not

458:    Level: advanced

460: .seealso: SVDTRLanczosSetOneSide()
461: @*/
462: PetscErrorCode SVDTRLanczosGetOneSide(SVD svd,PetscBool *oneside)
463: {

469:   PetscUseMethod(svd,"SVDTRLanczosGetOneSide_C",(SVD,PetscBool*),(svd,oneside));
470:   return(0);
471: }

473: PetscErrorCode SVDDestroy_TRLanczos(SVD svd)
474: {

478:   PetscFree(svd->data);
479:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",NULL);
480:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",NULL);
481:   return(0);
482: }

484: PetscErrorCode SVDView_TRLanczos(SVD svd,PetscViewer viewer)
485: {
487:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
488:   PetscBool      isascii;

491:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
492:   if (isascii) {
493:     PetscViewerASCIIPrintf(viewer,"  %s-sided reorthogonalization\n",lanczos->oneside? "one": "two");
494:   }
495:   return(0);
496: }

498: SLEPC_EXTERN PetscErrorCode SVDCreate_TRLanczos(SVD svd)
499: {
501:   SVD_TRLANCZOS  *ctx;

504:   PetscNewLog(svd,&ctx);
505:   svd->data = (void*)ctx;

507:   svd->ops->setup          = SVDSetUp_TRLanczos;
508:   svd->ops->solve          = SVDSolve_TRLanczos;
509:   svd->ops->destroy        = SVDDestroy_TRLanczos;
510:   svd->ops->setfromoptions = SVDSetFromOptions_TRLanczos;
511:   svd->ops->view           = SVDView_TRLanczos;
512:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",SVDTRLanczosSetOneSide_TRLanczos);
513:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",SVDTRLanczosGetOneSide_TRLanczos);
514:   return(0);
515: }