Actual source code: trlanczos.c

slepc-3.17.2 2022-08-09
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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>

 33: static PetscBool  cited = PETSC_FALSE,citedg = 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";
 43: static const char citationg[] =
 44:   "@Article{slepc-gsvd,\n"
 45:   "   author = \"F. Alvarruiz and C. Campos and J. E. Roman\",\n"
 46:   "   title = \"Thick-restarted {Lanczos} bidiagonalization methods for the {GSVD}\",\n"
 47:   "   note = \"In preparation\",\n"
 48:   "   year = \"2021\"\n"
 49:   "}\n";

 51: typedef struct {
 52:   /* user parameters */
 53:   PetscBool           oneside;   /* one-sided variant */
 54:   PetscReal           keep;      /* restart parameter */
 55:   PetscBool           lock;      /* locking/non-locking variant */
 56:   KSP                 ksp;       /* solver for least-squares problem in GSVD */
 57:   SVDTRLanczosGBidiag bidiag;    /* bidiagonalization variant for GSVD */
 58:   PetscBool           explicitmatrix;
 59:   /* auxiliary variables */
 60:   Mat                 Z;         /* aux matrix for GSVD, Z=[A;B] */
 61: } SVD_TRLANCZOS;

 63: /* Context for shell matrix [A; B] */
 64: typedef struct {
 65:   Mat      A,B,AT,BT;
 66:   Vec      y1,y2,y;
 67:   PetscInt m;
 68: } MatZData;

 70: static PetscErrorCode MatZCreateContext(SVD svd,MatZData **zdata)
 71: {
 72:   PetscNew(zdata);
 73:   (*zdata)->A = svd->A;
 74:   (*zdata)->B = svd->B;
 75:   (*zdata)->AT = svd->AT;
 76:   (*zdata)->BT = svd->BT;
 77:   MatCreateVecsEmpty(svd->A,NULL,&(*zdata)->y1);
 78:   MatCreateVecsEmpty(svd->B,NULL,&(*zdata)->y2);
 79:   VecGetLocalSize((*zdata)->y1,&(*zdata)->m);
 80:   BVCreateVec(svd->U,&(*zdata)->y);
 81:   PetscFunctionReturn(0);
 82: }

 84: static PetscErrorCode MatDestroy_Z(Mat Z)
 85: {
 86:   MatZData       *zdata;

 88:   MatShellGetContext(Z,&zdata);
 89:   VecDestroy(&zdata->y1);
 90:   VecDestroy(&zdata->y2);
 91:   VecDestroy(&zdata->y);
 92:   PetscFree(zdata);
 93:   PetscFunctionReturn(0);
 94: }

 96: static PetscErrorCode MatMult_Z(Mat Z,Vec x,Vec y)
 97: {
 98:   MatZData       *zdata;
 99:   PetscScalar    *y_elems;

101:   MatShellGetContext(Z,&zdata);
102:   VecGetArray(y,&y_elems);
103:   VecPlaceArray(zdata->y1,y_elems);
104:   VecPlaceArray(zdata->y2,y_elems+zdata->m);

106:   MatMult(zdata->A,x,zdata->y1);
107:   MatMult(zdata->B,x,zdata->y2);

109:   VecResetArray(zdata->y1);
110:   VecResetArray(zdata->y2);
111:   VecRestoreArray(y,&y_elems);
112:   PetscFunctionReturn(0);
113: }

115: static PetscErrorCode MatMultTranspose_Z(Mat Z,Vec y,Vec x)
116: {
117:   MatZData          *zdata;
118:   const PetscScalar *y_elems;

120:   MatShellGetContext(Z,&zdata);
121:   VecGetArrayRead(y,&y_elems);
122:   VecPlaceArray(zdata->y1,y_elems);
123:   VecPlaceArray(zdata->y2,y_elems+zdata->m);

125:   MatMult(zdata->AT,zdata->y1,x);
126:   MatMultAdd(zdata->BT,zdata->y2,x,x);

128:   VecResetArray(zdata->y1);
129:   VecResetArray(zdata->y2);
130:   VecRestoreArrayRead(y,&y_elems);
131:   PetscFunctionReturn(0);
132: }

134: static PetscErrorCode MatCreateVecs_Z(Mat Z,Vec *right,Vec *left)
135: {
136:   MatZData       *zdata;

138:   MatShellGetContext(Z,&zdata);
139:   if (right) MatCreateVecs(zdata->A,right,NULL);
140:   if (left) VecDuplicate(zdata->y,left);
141:   PetscFunctionReturn(0);
142: }

144: PetscErrorCode SVDSetUp_TRLanczos(SVD svd)
145: {
146:   PetscInt       N,m,n,p;
147:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
148:   DSType         dstype;
149:   MatZData       *zdata;
150:   Mat            mats[2],normal;
151:   MatType        Atype;
152:   PetscBool      sametype;

154:   MatGetSize(svd->A,NULL,&N);
155:   SVDSetDimensions_Default(svd);
158:   if (svd->max_it==PETSC_DEFAULT) svd->max_it = PetscMax(N/svd->ncv,100);
159:   if (!lanczos->keep) lanczos->keep = 0.5;
160:   svd->leftbasis = PETSC_TRUE;
161:   SVDAllocateSolution(svd,1);
162:   dstype = DSSVD;
163:   if (svd->isgeneralized) {
164:     if (lanczos->bidiag==SVD_TRLANCZOS_GBIDIAG_UPPER || lanczos->bidiag==SVD_TRLANCZOS_GBIDIAG_LOWER) dstype = DSGSVD;
165:     SVDSetWorkVecs(svd,1,1);

167:     if (svd->conv==SVD_CONV_ABS) {  /* Residual norms are multiplied by matrix norms */
168:       if (!svd->nrma) MatNorm(svd->OP,NORM_INFINITY,&svd->nrma);
169:       if (!svd->nrmb) MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb);
170:     }

172:     /* Create the matrix Z=[A;B] */
173:     MatDestroy(&lanczos->Z);
174:     MatGetLocalSize(svd->A,&m,&n);
175:     MatGetLocalSize(svd->B,&p,NULL);
176:     if (lanczos->explicitmatrix) {
177:       mats[0] = svd->A;
178:       mats[1] = svd->B;
179:       MatCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,1,NULL,mats,&lanczos->Z);
180:       MatGetType(svd->A,&Atype);
181:       PetscObjectTypeCompare((PetscObject)svd->B,Atype,&sametype);
182:       if (!sametype) Atype = MATAIJ;
183:       MatConvert(lanczos->Z,Atype,MAT_INPLACE_MATRIX,&lanczos->Z);
184:     } else {
185:       MatZCreateContext(svd,&zdata);
186:       MatCreateShell(PetscObjectComm((PetscObject)svd),m+p,n,PETSC_DECIDE,PETSC_DECIDE,zdata,&lanczos->Z);
187:       MatShellSetOperation(lanczos->Z,MATOP_MULT,(void(*)(void))MatMult_Z);
188: #if defined(PETSC_USE_COMPLEX)
189:       MatShellSetOperation(lanczos->Z,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultTranspose_Z);
190: #else
191:       MatShellSetOperation(lanczos->Z,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Z);
192: #endif
193:       MatShellSetOperation(lanczos->Z,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_Z);
194:       MatShellSetOperation(lanczos->Z,MATOP_DESTROY,(void(*)(void))MatDestroy_Z);
195:     }
196:     PetscLogObjectParent((PetscObject)svd,(PetscObject)lanczos->Z);

198:     /* create normal equations matrix, to build the preconditioner in LSQR */
199:     MatCreateNormalHermitian(lanczos->Z,&normal);

201:     if (!lanczos->ksp) SVDTRLanczosGetKSP(svd,&lanczos->ksp);
202:     SVD_KSPSetOperators(lanczos->ksp,lanczos->Z,normal);
203:     KSPSetUp(lanczos->ksp);
204:     MatDestroy(&normal);

206:     if (lanczos->oneside) PetscInfo(svd,"Warning: one-side option is ignored in GSVD\n");
207:   }
208:   DSSetType(svd->ds,dstype);
209:   DSSetCompact(svd->ds,PETSC_TRUE);
210:   DSSetExtraRow(svd->ds,PETSC_TRUE);
211:   DSAllocate(svd->ds,svd->ncv+1);
212:   PetscFunctionReturn(0);
213: }

215: static PetscErrorCode SVDOneSideTRLanczosMGS(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,BV U,PetscInt nconv,PetscInt l,PetscInt n,PetscScalar* work)
216: {
217:   PetscReal      a,b;
218:   PetscInt       i,k=nconv+l;
219:   Vec            ui,ui1,vi;

221:   BVGetColumn(V,k,&vi);
222:   BVGetColumn(U,k,&ui);
223:   MatMult(svd->A,vi,ui);
224:   BVRestoreColumn(V,k,&vi);
225:   BVRestoreColumn(U,k,&ui);
226:   if (l>0) {
227:     BVSetActiveColumns(U,nconv,n);
228:     for (i=0;i<l;i++) work[i]=beta[i+nconv];
229:     BVMultColumn(U,-1.0,1.0,k,work);
230:   }
231:   BVNormColumn(U,k,NORM_2,&a);
232:   BVScaleColumn(U,k,1.0/a);
233:   alpha[k] = a;

235:   for (i=k+1;i<n;i++) {
236:     BVGetColumn(V,i,&vi);
237:     BVGetColumn(U,i-1,&ui1);
238:     MatMult(svd->AT,ui1,vi);
239:     BVRestoreColumn(V,i,&vi);
240:     BVRestoreColumn(U,i-1,&ui1);
241:     BVOrthonormalizeColumn(V,i,PETSC_FALSE,&b,NULL);
242:     beta[i-1] = b;

244:     BVGetColumn(V,i,&vi);
245:     BVGetColumn(U,i,&ui);
246:     MatMult(svd->A,vi,ui);
247:     BVRestoreColumn(V,i,&vi);
248:     BVGetColumn(U,i-1,&ui1);
249:     VecAXPY(ui,-b,ui1);
250:     BVRestoreColumn(U,i-1,&ui1);
251:     BVRestoreColumn(U,i,&ui);
252:     BVNormColumn(U,i,NORM_2,&a);
253:     BVScaleColumn(U,i,1.0/a);
254:     alpha[i] = a;
255:   }

257:   BVGetColumn(V,n,&vi);
258:   BVGetColumn(U,n-1,&ui1);
259:   MatMult(svd->AT,ui1,vi);
260:   BVRestoreColumn(V,n,&vi);
261:   BVRestoreColumn(U,n-1,&ui1);
262:   BVOrthogonalizeColumn(V,n,NULL,&b,NULL);
263:   beta[n-1] = b;
264:   PetscFunctionReturn(0);
265: }

267: /*
268:   Custom CGS orthogonalization, preprocess after first orthogonalization
269: */
270: static PetscErrorCode SVDOrthogonalizeCGS(BV V,PetscInt i,PetscScalar* h,PetscReal a,BVOrthogRefineType refine,PetscReal eta,PetscReal *norm)
271: {
272:   PetscReal      sum,onorm;
273:   PetscScalar    dot;
274:   PetscInt       j;

276:   switch (refine) {
277:   case BV_ORTHOG_REFINE_NEVER:
278:     BVNormColumn(V,i,NORM_2,norm);
279:     break;
280:   case BV_ORTHOG_REFINE_ALWAYS:
281:     BVSetActiveColumns(V,0,i);
282:     BVDotColumn(V,i,h);
283:     BVMultColumn(V,-1.0,1.0,i,h);
284:     BVNormColumn(V,i,NORM_2,norm);
285:     break;
286:   case BV_ORTHOG_REFINE_IFNEEDED:
287:     dot = h[i];
288:     onorm = PetscSqrtReal(PetscRealPart(dot)) / a;
289:     sum = 0.0;
290:     for (j=0;j<i;j++) {
291:       sum += PetscRealPart(h[j] * PetscConj(h[j]));
292:     }
293:     *norm = PetscRealPart(dot)/(a*a) - sum;
294:     if (*norm>0.0) *norm = PetscSqrtReal(*norm);
295:     else BVNormColumn(V,i,NORM_2,norm);
296:     if (*norm < eta*onorm) {
297:       BVSetActiveColumns(V,0,i);
298:       BVDotColumn(V,i,h);
299:       BVMultColumn(V,-1.0,1.0,i,h);
300:       BVNormColumn(V,i,NORM_2,norm);
301:     }
302:     break;
303:   }
304:   PetscFunctionReturn(0);
305: }

307: static PetscErrorCode SVDOneSideTRLanczosCGS(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,BV U,PetscInt nconv,PetscInt l,PetscInt n,PetscScalar* work)
308: {
309:   PetscReal          a,b,eta;
310:   PetscInt           i,j,k=nconv+l;
311:   Vec                ui,ui1,vi;
312:   BVOrthogRefineType refine;

314:   BVGetColumn(V,k,&vi);
315:   BVGetColumn(U,k,&ui);
316:   MatMult(svd->A,vi,ui);
317:   BVRestoreColumn(V,k,&vi);
318:   BVRestoreColumn(U,k,&ui);
319:   if (l>0) {
320:     BVSetActiveColumns(U,nconv,n);
321:     for (i=0;i<l;i++) work[i]=beta[i+nconv];
322:     BVMultColumn(U,-1.0,1.0,k,work);
323:   }
324:   BVGetOrthogonalization(V,NULL,&refine,&eta,NULL);

326:   for (i=k+1;i<n;i++) {
327:     BVGetColumn(V,i,&vi);
328:     BVGetColumn(U,i-1,&ui1);
329:     MatMult(svd->AT,ui1,vi);
330:     BVRestoreColumn(V,i,&vi);
331:     BVRestoreColumn(U,i-1,&ui1);
332:     BVNormColumnBegin(U,i-1,NORM_2,&a);
333:     if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
334:       BVSetActiveColumns(V,0,i+1);
335:       BVGetColumn(V,i,&vi);
336:       BVDotVecBegin(V,vi,work);
337:     } else {
338:       BVSetActiveColumns(V,0,i);
339:       BVDotColumnBegin(V,i,work);
340:     }
341:     BVNormColumnEnd(U,i-1,NORM_2,&a);
342:     if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
343:       BVDotVecEnd(V,vi,work);
344:       BVRestoreColumn(V,i,&vi);
345:       BVSetActiveColumns(V,0,i);
346:     } else BVDotColumnEnd(V,i,work);

348:     BVScaleColumn(U,i-1,1.0/a);
349:     for (j=0;j<i;j++) work[j] = work[j] / a;
350:     BVMultColumn(V,-1.0,1.0/a,i,work);
351:     SVDOrthogonalizeCGS(V,i,work,a,refine,eta,&b);
352:     BVScaleColumn(V,i,1.0/b);

355:     BVGetColumn(V,i,&vi);
356:     BVGetColumn(U,i,&ui);
357:     BVGetColumn(U,i-1,&ui1);
358:     MatMult(svd->A,vi,ui);
359:     VecAXPY(ui,-b,ui1);
360:     BVRestoreColumn(V,i,&vi);
361:     BVRestoreColumn(U,i,&ui);
362:     BVRestoreColumn(U,i-1,&ui1);

364:     alpha[i-1] = a;
365:     beta[i-1] = b;
366:   }

368:   BVGetColumn(V,n,&vi);
369:   BVGetColumn(U,n-1,&ui1);
370:   MatMult(svd->AT,ui1,vi);
371:   BVRestoreColumn(V,n,&vi);
372:   BVRestoreColumn(U,n-1,&ui1);

374:   BVNormColumnBegin(svd->U,n-1,NORM_2,&a);
375:   if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
376:     BVSetActiveColumns(V,0,n+1);
377:     BVGetColumn(V,n,&vi);
378:     BVDotVecBegin(V,vi,work);
379:   } else {
380:     BVSetActiveColumns(V,0,n);
381:     BVDotColumnBegin(V,n,work);
382:   }
383:   BVNormColumnEnd(svd->U,n-1,NORM_2,&a);
384:   if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
385:     BVDotVecEnd(V,vi,work);
386:     BVRestoreColumn(V,n,&vi);
387:   } else BVDotColumnEnd(V,n,work);

389:   BVScaleColumn(U,n-1,1.0/a);
390:   for (j=0;j<n;j++) work[j] = work[j] / a;
391:   BVMultColumn(V,-1.0,1.0/a,n,work);
392:   SVDOrthogonalizeCGS(V,n,work,a,refine,eta,&b);
393:   BVSetActiveColumns(V,nconv,n);
394:   alpha[n-1] = a;
395:   beta[n-1] = b;
396:   PetscFunctionReturn(0);
397: }

399: PetscErrorCode SVDSolve_TRLanczos(SVD svd)
400: {
401:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
402:   PetscReal      *alpha,*beta;
403:   PetscScalar    *swork=NULL,*w;
404:   PetscInt       i,k,l,nv,ld;
405:   Mat            U,V;
406:   PetscBool      breakdown=PETSC_FALSE;
407:   BVOrthogType   orthog;

409:   PetscCitationsRegister(citation,&cited);
410:   /* allocate working space */
411:   DSGetLeadingDimension(svd->ds,&ld);
412:   BVGetOrthogonalization(svd->V,&orthog,NULL,NULL,NULL);
413:   PetscMalloc1(ld,&w);
414:   if (lanczos->oneside) PetscMalloc1(svd->ncv+1,&swork);

416:   /* normalize start vector */
417:   if (!svd->nini) {
418:     BVSetRandomColumn(svd->V,0);
419:     BVOrthonormalizeColumn(svd->V,0,PETSC_TRUE,NULL,NULL);
420:   }

422:   l = 0;
423:   while (svd->reason == SVD_CONVERGED_ITERATING) {
424:     svd->its++;

426:     /* inner loop */
427:     nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
428:     DSGetArrayReal(svd->ds,DS_MAT_T,&alpha);
429:     beta = alpha + ld;
430:     if (lanczos->oneside) {
431:       if (orthog == BV_ORTHOG_MGS) SVDOneSideTRLanczosMGS(svd,alpha,beta,svd->V,svd->U,svd->nconv,l,nv,swork);
432:       else SVDOneSideTRLanczosCGS(svd,alpha,beta,svd->V,svd->U,svd->nconv,l,nv,swork);
433:     } else SVDTwoSideLanczos(svd,alpha,beta,svd->V,svd->U,svd->nconv+l,&nv,&breakdown);
434:     DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha);
435:     BVScaleColumn(svd->V,nv,1.0/beta[nv-1]);
436:     BVSetActiveColumns(svd->V,svd->nconv,nv);
437:     BVSetActiveColumns(svd->U,svd->nconv,nv);

439:     /* solve projected problem */
440:     DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l);
441:     DSSVDSetDimensions(svd->ds,nv);
442:     DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE);
443:     DSSolve(svd->ds,w,NULL);
444:     DSSort(svd->ds,w,NULL,NULL,NULL,NULL);
445:     DSUpdateExtraRow(svd->ds);
446:     DSSynchronize(svd->ds,w,NULL);
447:     for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);

449:     /* check convergence */
450:     SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,1.0,&k);
451:     (*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx);

453:     /* update l */
454:     if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
455:     else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
456:     if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
457:     if (l) PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l);

459:     if (svd->reason == SVD_CONVERGED_ITERATING) {
460:       if (PetscUnlikely(breakdown || k==nv)) {
461:         /* Start a new bidiagonalization */
462:         PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its);
463:         if (k<svd->nsv) {
464:           BVSetRandomColumn(svd->V,k);
465:           BVOrthonormalizeColumn(svd->V,k,PETSC_FALSE,NULL,&breakdown);
466:           if (breakdown) {
467:             svd->reason = SVD_DIVERGED_BREAKDOWN;
468:             PetscInfo(svd,"Unable to generate more start vectors\n");
469:           }
470:         }
471:       } else DSTruncate(svd->ds,k+l,PETSC_FALSE);
472:     }

474:     /* compute converged singular vectors and restart vectors */
475:     DSGetMat(svd->ds,DS_MAT_V,&V);
476:     BVMultInPlace(svd->V,V,svd->nconv,k+l);
477:     MatDestroy(&V);
478:     DSGetMat(svd->ds,DS_MAT_U,&U);
479:     BVMultInPlace(svd->U,U,svd->nconv,k+l);
480:     MatDestroy(&U);

482:     /* copy the last vector to be the next initial vector */
483:     if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) BVCopyColumn(svd->V,nv,k+l);

485:     svd->nconv = k;
486:     SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv);
487:   }

489:   /* orthonormalize U columns in one side method */
490:   if (lanczos->oneside) {
491:     for (i=0;i<svd->nconv;i++) BVOrthonormalizeColumn(svd->U,i,PETSC_FALSE,NULL,NULL);
492:   }

494:   /* free working space */
495:   PetscFree(w);
496:   if (swork) PetscFree(swork);
497:   DSTruncate(svd->ds,svd->nconv,PETSC_TRUE);
498:   PetscFunctionReturn(0);
499: }

501: static PetscErrorCode SVDTwoSideLanczosGSingle(SVD svd,PetscReal *alpha,PetscReal *beta,Mat Z,BV V,BV U,KSP ksp,PetscInt k,PetscInt *n,PetscBool *breakdown)
502: {
503:   PetscInt          i,j,m;
504:   const PetscScalar *carr;
505:   PetscScalar       *arr;
506:   Vec               u,v,ut=svd->workl[0],x=svd->workr[0],v1;
507:   PetscBool         lindep=PETSC_FALSE;

509:   MatCreateVecsEmpty(svd->A,NULL,&v1);
510:   BVGetColumn(V,k,&v);
511:   BVGetColumn(U,k,&u);

513:   /* Form ut=[u;0] */
514:   VecZeroEntries(ut);
515:   VecGetLocalSize(u,&m);
516:   VecGetArrayRead(u,&carr);
517:   VecGetArray(ut,&arr);
518:   for (j=0; j<m; j++) arr[j] = carr[j];
519:   VecRestoreArrayRead(u,&carr);
520:   VecRestoreArray(ut,&arr);

522:   /* Solve least squares problem */
523:   KSPSolve(ksp,ut,x);

525:   MatMult(Z,x,v);

527:   BVRestoreColumn(U,k,&u);
528:   BVRestoreColumn(V,k,&v);
529:   BVOrthonormalizeColumn(V,k,PETSC_FALSE,alpha+k,&lindep);
530:   if (PetscUnlikely(lindep)) {
531:     *n = k;
532:     if (breakdown) *breakdown = lindep;
533:     PetscFunctionReturn(0);
534:   }

536:   for (i=k+1; i<*n; i++) {

538:     /* Compute vector i of BV U */
539:     BVGetColumn(V,i-1,&v);
540:     VecGetArray(v,&arr);
541:     VecPlaceArray(v1,arr);
542:     VecRestoreArray(v,&arr);
543:     BVRestoreColumn(V,i-1,&v);
544:     BVInsertVec(U,i,v1);
545:     VecResetArray(v1);
546:     BVOrthonormalizeColumn(U,i,PETSC_FALSE,beta+i-1,&lindep);
547:     if (PetscUnlikely(lindep)) {
548:       *n = i;
549:       break;
550:     }

552:     /* Compute vector i of BV V */

554:     BVGetColumn(V,i,&v);
555:     BVGetColumn(U,i,&u);

557:     /* Form ut=[u;0] */
558:     VecGetArrayRead(u,&carr);
559:     VecGetArray(ut,&arr);
560:     for (j=0; j<m; j++) arr[j] = carr[j];
561:     VecRestoreArrayRead(u,&carr);
562:     VecRestoreArray(ut,&arr);

564:     /* Solve least squares problem */
565:     KSPSolve(ksp,ut,x);

567:     MatMult(Z,x,v);

569:     BVRestoreColumn(U,i,&u);
570:     BVRestoreColumn(V,i,&v);
571:     BVOrthonormalizeColumn(V,i,PETSC_FALSE,alpha+i,&lindep);
572:     if (PetscUnlikely(lindep)) {
573:       *n = i;
574:       break;
575:     }
576:   }

578:   /* Compute vector n of BV U */
579:   if (!lindep) {
580:     BVGetColumn(V,*n-1,&v);
581:     VecGetArray(v,&arr);
582:     VecPlaceArray(v1,arr);
583:     VecRestoreArray(v,&arr);
584:     BVRestoreColumn(V,*n-1,&v);
585:     BVInsertVec(U,*n,v1);
586:     VecResetArray(v1);
587:     BVOrthonormalizeColumn(U,*n,PETSC_FALSE,beta+*n-1,&lindep);
588:   }
589:   if (breakdown) *breakdown = lindep;
590:   VecDestroy(&v1);
591:   PetscFunctionReturn(0);
592: }

594: /* solve generalized problem with single bidiagonalization of Q_A */
595: PetscErrorCode SVDSolve_TRLanczosGSingle(SVD svd,BV U1,BV V)
596: {
597:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
598:   PetscReal      *alpha,*beta,normr;
599:   PetscScalar    *w;
600:   PetscInt       i,k,l,nv,ld;
601:   Mat            U,VV;
602:   PetscBool      breakdown=PETSC_FALSE;

604:   DSGetLeadingDimension(svd->ds,&ld);
605:   PetscMalloc1(ld,&w);
606:   normr = (svd->conv==SVD_CONV_ABS)? PetscMax(svd->nrma,svd->nrmb): 1.0;

608:   /* normalize start vector */
609:   if (!svd->ninil) {
610:     BVSetRandomColumn(U1,0);
611:     BVOrthonormalizeColumn(U1,0,PETSC_TRUE,NULL,NULL);
612:   }

614:   l = 0;
615:   while (svd->reason == SVD_CONVERGED_ITERATING) {
616:     svd->its++;

618:     /* inner loop */
619:     nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
620:     DSGetArrayReal(svd->ds,DS_MAT_T,&alpha);
621:     beta = alpha + ld;
622:     SVDTwoSideLanczosGSingle(svd,alpha,beta,lanczos->Z,V,U1,lanczos->ksp,svd->nconv+l,&nv,&breakdown);
623:     DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha);
624:     BVSetActiveColumns(V,svd->nconv,nv);
625:     BVSetActiveColumns(U1,svd->nconv,nv);

627:     /* solve projected problem */
628:     DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l);
629:     DSSVDSetDimensions(svd->ds,nv);
630:     DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE);
631:     DSSolve(svd->ds,w,NULL);
632:     DSSort(svd->ds,w,NULL,NULL,NULL,NULL);
633:     DSUpdateExtraRow(svd->ds);
634:     DSSynchronize(svd->ds,w,NULL);
635:     for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);

637:     /* check convergence */
638:     SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,normr,&k);
639:     (*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx);

641:     /* update l */
642:     if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
643:     else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
644:     if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
645:     if (l) PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l);

647:     if (svd->reason == SVD_CONVERGED_ITERATING) {
648:       if (PetscUnlikely(breakdown || k==nv)) {
649:         /* Start a new bidiagonalization */
650:         PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its);
651:         if (k<svd->nsv) {
652:           BVSetRandomColumn(U1,k);
653:           BVOrthonormalizeColumn(U1,k,PETSC_FALSE,NULL,&breakdown);
654:           if (breakdown) {
655:             svd->reason = SVD_DIVERGED_BREAKDOWN;
656:             PetscInfo(svd,"Unable to generate more start vectors\n");
657:           }
658:         }
659:       } else DSTruncate(svd->ds,k+l,PETSC_FALSE);
660:     }

662:     /* compute converged singular vectors and restart vectors */
663:     DSGetMat(svd->ds,DS_MAT_U,&U);
664:     BVMultInPlace(V,U,svd->nconv,k+l);
665:     MatDestroy(&U);
666:     DSGetMat(svd->ds,DS_MAT_V,&VV);
667:     BVMultInPlace(U1,VV,svd->nconv,k+l);
668:     MatDestroy(&VV);

670:     /* copy the last vector to be the next initial vector */
671:     if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) BVCopyColumn(U1,nv,k+l);

673:     svd->nconv = k;
674:     SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv);
675:   }

677:   PetscFree(w);
678:   PetscFunctionReturn(0);
679: }

681: /* Move generalized left singular vectors (0..nconv) from U1 and U2 to its final destination svd->U (single variant) */
682: static inline PetscErrorCode SVDLeftSingularVectors_Single(SVD svd,BV U1,BV U2)
683: {
684:   PetscInt          i,k,m,p;
685:   Vec               u,u1,u2;
686:   PetscScalar       *ua,*u2a;
687:   const PetscScalar *u1a;
688:   PetscReal         s;

690:   MatGetLocalSize(svd->A,&m,NULL);
691:   MatGetLocalSize(svd->B,&p,NULL);
692:   for (i=0;i<svd->nconv;i++) {
693:     BVGetColumn(U1,i,&u1);
694:     BVGetColumn(U2,i,&u2);
695:     BVGetColumn(svd->U,i,&u);
696:     VecGetArrayRead(u1,&u1a);
697:     VecGetArray(u,&ua);
698:     VecGetArray(u2,&u2a);
699:     /* Copy column from U1 to upper part of u */
700:     for (k=0;k<m;k++) ua[k] = u1a[k];
701:     /* Copy column from lower part of U to U2. Orthogonalize column in U2 and copy back to U */
702:     for (k=0;k<p;k++) u2a[k] = ua[m+k];
703:     VecRestoreArray(u2,&u2a);
704:     BVRestoreColumn(U2,i,&u2);
705:     BVOrthonormalizeColumn(U2,i,PETSC_FALSE,&s,NULL);
706:     BVGetColumn(U2,i,&u2);
707:     VecGetArray(u2,&u2a);
708:     for (k=0;k<p;k++) ua[m+k] = u2a[k];
709:     /* Update singular value */
710:     svd->sigma[i] /= s;
711:     VecRestoreArrayRead(u1,&u1a);
712:     VecRestoreArray(u,&ua);
713:     VecRestoreArray(u2,&u2a);
714:     BVRestoreColumn(U1,i,&u1);
715:     BVRestoreColumn(U2,i,&u2);
716:     BVRestoreColumn(svd->U,i,&u);
717:   }
718:   PetscFunctionReturn(0);
719: }

721: static PetscErrorCode SVDTwoSideLanczosGUpper(SVD svd,PetscReal *alpha,PetscReal *beta,PetscReal *alphah,PetscReal *betah,Mat Z,BV U1,BV U2,BV V,KSP ksp,PetscInt k,PetscInt *n,PetscBool *breakdown)
722: {
723:   PetscInt          i,j,m,p;
724:   const PetscScalar *carr;
725:   PetscScalar       *arr,*u2arr;
726:   Vec               u,v,ut=svd->workl[0],x=svd->workr[0],v1,u2;
727:   PetscBool         lindep=PETSC_FALSE,lindep1=PETSC_FALSE,lindep2=PETSC_FALSE;

729:   MatCreateVecsEmpty(svd->A,NULL,&v1);
730:   MatGetLocalSize(svd->A,&m,NULL);
731:   MatGetLocalSize(svd->B,&p,NULL);

733:   for (i=k; i<*n; i++) {
734:     /* Compute vector i of BV U1 */
735:     BVGetColumn(V,i,&v);
736:     VecGetArrayRead(v,&carr);
737:     VecPlaceArray(v1,carr);
738:     BVInsertVec(U1,i,v1);
739:     VecResetArray(v1);
740:     BVOrthonormalizeColumn(U1,i,PETSC_FALSE,alpha+i,&lindep1);

742:     /* Compute vector i of BV U2 */
743:     BVGetColumn(U2,i,&u2);
744:     VecGetArray(u2,&u2arr);
745:     if (i%2) {
746:       for (j=0; j<p; j++) u2arr[j] = -carr[m+j];
747:     } else {
748:       for (j=0; j<p; j++) u2arr[j] = carr[m+j];
749:     }
750:     VecRestoreArray(u2,&u2arr);
751:     BVRestoreColumn(U2,i,&u2);
752:     VecRestoreArrayRead(v,&carr);
753:     BVRestoreColumn(V,i,&v);
754:     BVOrthonormalizeColumn(U2,i,PETSC_FALSE,alphah+i,&lindep2);
755:     if (i%2) alphah[i] = -alphah[i];
756:     if (PetscUnlikely(lindep1 || lindep2)) {
757:       lindep = PETSC_TRUE;
758:       *n = i;
759:       break;
760:     }

762:     /* Compute vector i+1 of BV V */
763:     BVGetColumn(V,i+1,&v);
764:     /* Form ut=[u;0] */
765:     BVGetColumn(U1,i,&u);
766:     VecZeroEntries(ut);
767:     VecGetArrayRead(u,&carr);
768:     VecGetArray(ut,&arr);
769:     for (j=0; j<m; j++) arr[j] = carr[j];
770:     VecRestoreArrayRead(u,&carr);
771:     VecRestoreArray(ut,&arr);
772:     /* Solve least squares problem */
773:     KSPSolve(ksp,ut,x);
774:     MatMult(Z,x,v);
775:     BVRestoreColumn(U1,i,&u);
776:     BVRestoreColumn(V,i+1,&v);
777:     BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,beta+i,&lindep);
778:     betah[i] = -alpha[i]*beta[i]/alphah[i];
779:     if (PetscUnlikely(lindep)) {
780:       *n = i;
781:       break;
782:     }
783:   }
784:   if (breakdown) *breakdown = lindep;
785:   VecDestroy(&v1);
786:   PetscFunctionReturn(0);
787: }

789: /* generate random initial vector in column k for joint upper-upper bidiagonalization */
790: static inline PetscErrorCode SVDInitialVectorGUpper(SVD svd,BV V,BV U1,PetscInt k,PetscBool *breakdown)
791: {
792:   SVD_TRLANCZOS     *lanczos = (SVD_TRLANCZOS*)svd->data;
793:   Vec               u,v,ut=svd->workl[0],x=svd->workr[0];
794:   PetscInt          m,j;
795:   PetscScalar       *arr;
796:   const PetscScalar *carr;

798:   /* Form ut=[u;0] where u is the k-th column of U1 */
799:   VecZeroEntries(ut);
800:   BVGetColumn(U1,k,&u);
801:   VecGetLocalSize(u,&m);
802:   VecGetArrayRead(u,&carr);
803:   VecGetArray(ut,&arr);
804:   for (j=0; j<m; j++) arr[j] = carr[j];
805:   VecRestoreArrayRead(u,&carr);
806:   VecRestoreArray(ut,&arr);
807:   BVRestoreColumn(U1,k,&u);
808:   /* Solve least squares problem Z*x=ut for x. Then set v=Z*x */
809:   KSPSolve(lanczos->ksp,ut,x);
810:   BVGetColumn(V,k,&v);
811:   MatMult(lanczos->Z,x,v);
812:   BVRestoreColumn(V,k,&v);
813:   if (breakdown) BVOrthonormalizeColumn(V,k,PETSC_FALSE,NULL,breakdown);
814:   else BVOrthonormalizeColumn(V,k,PETSC_TRUE,NULL,NULL);
815:   PetscFunctionReturn(0);
816: }

818: /* solve generalized problem with joint upper-upper bidiagonalization */
819: PetscErrorCode SVDSolve_TRLanczosGUpper(SVD svd,BV U1,BV U2,BV V)
820: {
821:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
822:   PetscReal      *alpha,*beta,*alphah,*betah,normr;
823:   PetscScalar    *w;
824:   PetscInt       i,k,l,nv,ld;
825:   Mat            U,Vmat,X;
826:   PetscBool      breakdown=PETSC_FALSE;

828:   DSGetLeadingDimension(svd->ds,&ld);
829:   PetscMalloc1(ld,&w);
830:   normr = (svd->conv==SVD_CONV_ABS)? PetscMax(svd->nrma,svd->nrmb): 1.0;

832:   /* normalize start vector */
833:   if (!svd->ninil) BVSetRandomColumn(U1,0);
834:   SVDInitialVectorGUpper(svd,V,U1,0,NULL);

836:   l = 0;
837:   while (svd->reason == SVD_CONVERGED_ITERATING) {
838:     svd->its++;

840:     /* inner loop */
841:     nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
842:     DSGetArrayReal(svd->ds,DS_MAT_T,&alpha);
843:     DSGetArrayReal(svd->ds,DS_MAT_D,&alphah);
844:     beta = alpha + ld;
845:     betah = alpha + 2*ld;
846:     SVDTwoSideLanczosGUpper(svd,alpha,beta,alphah,betah,lanczos->Z,U1,U2,V,lanczos->ksp,svd->nconv+l,&nv,&breakdown);
847:     DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha);
848:     DSRestoreArrayReal(svd->ds,DS_MAT_D,&alphah);
849:     BVSetActiveColumns(V,svd->nconv,nv);
850:     BVSetActiveColumns(U1,svd->nconv,nv);
851:     BVSetActiveColumns(U2,svd->nconv,nv);

853:     /* solve projected problem */
854:     DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l);
855:     DSGSVDSetDimensions(svd->ds,nv,nv);
856:     DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE);
857:     DSSolve(svd->ds,w,NULL);
858:     DSSort(svd->ds,w,NULL,NULL,NULL,NULL);
859:     DSUpdateExtraRow(svd->ds);
860:     DSSynchronize(svd->ds,w,NULL);
861:     for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);

863:     /* check convergence */
864:     SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,normr,&k);
865:     (*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx);

867:     /* update l */
868:     if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
869:     else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
870:     if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
871:     if (l) PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l);

873:     if (svd->reason == SVD_CONVERGED_ITERATING) {
874:       if (PetscUnlikely(breakdown || k==nv)) {
875:         /* Start a new bidiagonalization */
876:         PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its);
877:         if (k<svd->nsv) {
878:           BVSetRandomColumn(U1,k);
879:           SVDInitialVectorGUpper(svd,V,U1,k,&breakdown);
880:           if (breakdown) {
881:             svd->reason = SVD_DIVERGED_BREAKDOWN;
882:             PetscInfo(svd,"Unable to generate more start vectors\n");
883:           }
884:         }
885:       } else DSTruncate(svd->ds,k+l,PETSC_FALSE);
886:     }
887:     /* compute converged singular vectors and restart vectors */
888:     DSGetMat(svd->ds,DS_MAT_X,&X);
889:     BVMultInPlace(V,X,svd->nconv,k+l);
890:     MatDestroy(&X);
891:     DSGetMat(svd->ds,DS_MAT_U,&U);
892:     BVMultInPlace(U1,U,svd->nconv,k+l);
893:     MatDestroy(&U);
894:     DSGetMat(svd->ds,DS_MAT_V,&Vmat);
895:     BVMultInPlace(U2,Vmat,svd->nconv,k+l);
896:     MatDestroy(&Vmat);

898:     /* copy the last vector to be the next initial vector */
899:     if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) BVCopyColumn(V,nv,k+l);

901:     svd->nconv = k;
902:     SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv);
903:   }

905:   PetscFree(w);
906:   PetscFunctionReturn(0);
907: }

909: /* Move generalized left singular vectors (0..nconv) from U1 and U2 to its final destination svd->U (upper and lower variants) */
910: static inline PetscErrorCode SVDLeftSingularVectors(SVD svd,BV U1,BV U2)
911: {
912:   PetscInt          i,k,m,p;
913:   Vec               u,u1,u2;
914:   PetscScalar       *ua;
915:   const PetscScalar *u1a,*u2a;

917:   MatGetLocalSize(svd->A,&m,NULL);
918:   MatGetLocalSize(svd->B,&p,NULL);
919:   for (i=0;i<svd->nconv;i++) {
920:     BVGetColumn(U1,i,&u1);
921:     BVGetColumn(U2,i,&u2);
922:     BVGetColumn(svd->U,i,&u);
923:     VecGetArrayRead(u1,&u1a);
924:     VecGetArrayRead(u2,&u2a);
925:     VecGetArray(u,&ua);
926:     /* Copy column from u1 to upper part of u */
927:     for (k=0;k<m;k++) ua[k] = u1a[k];
928:     /* Copy column from u2 to lower part of u */
929:     for (k=0;k<p;k++) ua[m+k] = u2a[k];
930:     VecRestoreArrayRead(u1,&u1a);
931:     VecRestoreArrayRead(u2,&u2a);
932:     VecRestoreArray(u,&ua);
933:     BVRestoreColumn(U1,i,&u1);
934:     BVRestoreColumn(U2,i,&u2);
935:     BVRestoreColumn(svd->U,i,&u);
936:   }
937:   PetscFunctionReturn(0);
938: }

940: static PetscErrorCode SVDTwoSideLanczosGLower(SVD svd,PetscReal *alpha,PetscReal *beta,PetscReal *alphah,PetscReal *betah,Mat Z,BV U1,BV U2,BV V,KSP ksp,PetscInt k,PetscInt *n,PetscBool *breakdown)
941: {
942:   PetscInt          i,j,m,p;
943:   const PetscScalar *carr;
944:   PetscScalar       *arr,*u2arr;
945:   Vec               u,v,ut=svd->workl[0],x=svd->workr[0],v1,u2;
946:   PetscBool         lindep=PETSC_FALSE;

948:   MatCreateVecsEmpty(svd->A,NULL,&v1);
949:   MatGetLocalSize(svd->A,&m,NULL);
950:   MatGetLocalSize(svd->B,&p,NULL);

952:   for (i=k; i<*n; i++) {
953:     /* Compute vector i of BV U2 */
954:     BVGetColumn(V,i,&v);
955:     VecGetArrayRead(v,&carr);
956:     BVGetColumn(U2,i,&u2);
957:     VecGetArray(u2,&u2arr);
958:     if (i%2) {
959:       for (j=0; j<p; j++) u2arr[j] = -carr[m+j];
960:     } else {
961:       for (j=0; j<p; j++) u2arr[j] = carr[m+j];
962:     }
963:     VecRestoreArray(u2,&u2arr);
964:     BVRestoreColumn(U2,i,&u2);
965:     BVOrthonormalizeColumn(U2,i,PETSC_FALSE,alphah+i,&lindep);
966:     if (i%2) alphah[i] = -alphah[i];
967:     if (PetscUnlikely(lindep)) {
968:       BVRestoreColumn(V,i,&v);
969:       *n = i;
970:       break;
971:     }

973:     /* Compute vector i+1 of BV U1 */
974:     VecPlaceArray(v1,carr);
975:     BVInsertVec(U1,i+1,v1);
976:     VecResetArray(v1);
977:     BVOrthonormalizeColumn(U1,i+1,PETSC_FALSE,beta+i,&lindep);
978:     VecRestoreArrayRead(v,&carr);
979:     BVRestoreColumn(V,i,&v);
980:     if (PetscUnlikely(lindep)) {
981:       *n = i+1;
982:       break;
983:     }

985:     /* Compute vector i+1 of BV V */
986:     BVGetColumn(V,i+1,&v);
987:     /* Form ut=[u;0] where u is column i+1 of BV U1 */
988:     BVGetColumn(U1,i+1,&u);
989:     VecZeroEntries(ut);
990:     VecGetArrayRead(u,&carr);
991:     VecGetArray(ut,&arr);
992:     for (j=0; j<m; j++) arr[j] = carr[j];
993:     VecRestoreArrayRead(u,&carr);
994:     VecRestoreArray(ut,&arr);
995:     /* Solve least squares problem */
996:     KSPSolve(ksp,ut,x);
997:     MatMult(Z,x,v);
998:     BVRestoreColumn(U1,i+1,&u);
999:     BVRestoreColumn(V,i+1,&v);
1000:     BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,alpha+i+1,&lindep);
1001:     betah[i] = -alpha[i+1]*beta[i]/alphah[i];
1002:     if (PetscUnlikely(lindep)) {
1003:       *n = i+1;
1004:       break;
1005:     }
1006:   }
1007:   if (breakdown) *breakdown = lindep;
1008:   VecDestroy(&v1);
1009:   PetscFunctionReturn(0);
1010: }

1012: /* generate random initial vector in column k for joint lower-upper bidiagonalization */
1013: static inline PetscErrorCode SVDInitialVectorGLower(SVD svd,BV V,BV U1,BV U2,PetscInt k,PetscBool *breakdown)
1014: {
1015:   SVD_TRLANCZOS     *lanczos = (SVD_TRLANCZOS*)svd->data;
1016:   const PetscScalar *carr;
1017:   PetscScalar       *arr;
1018:   PetscReal         *alpha;
1019:   PetscInt          j,m,p;
1020:   Vec               u,uh,v,ut=svd->workl[0],x=svd->workr[0];

1022:   MatGetLocalSize(svd->A,&m,NULL);
1023:   MatGetLocalSize(svd->B,&p,NULL);
1024:   /* Form ut=[0;uh], where uh is the k-th column of U2 */
1025:   BVGetColumn(U2,k,&uh);
1026:   VecZeroEntries(ut);
1027:   VecGetArrayRead(uh,&carr);
1028:   VecGetArray(ut,&arr);
1029:   for (j=0; j<p; j++) arr[m+j] = carr[j];
1030:   VecRestoreArrayRead(uh,&carr);
1031:   VecRestoreArray(ut,&arr);
1032:   BVRestoreColumn(U2,k,&uh);
1033:   /* Solve least squares problem Z*x=ut for x. Then set ut=Z*x */
1034:   KSPSolve(lanczos->ksp,ut,x);
1035:   MatMult(lanczos->Z,x,ut);
1036:   /* Form u, column k of BV U1, as the upper part of ut and orthonormalize */
1037:   MatCreateVecsEmpty(svd->A,NULL,&u);
1038:   VecGetArrayRead(ut,&carr);
1039:   VecPlaceArray(u,carr);
1040:   BVInsertVec(U1,k,u);
1041:   VecResetArray(u);
1042:   VecRestoreArrayRead(ut,&carr);
1043:   VecDestroy(&u);
1044:   if (breakdown) BVOrthonormalizeColumn(U1,k,PETSC_FALSE,NULL,breakdown);
1045:   else BVOrthonormalizeColumn(U1,k,PETSC_TRUE,NULL,NULL);

1047:   if (!breakdown || !*breakdown) {
1048:     MatGetLocalSize(svd->A,&m,NULL);
1049:     /* Compute k-th vector of BV V */
1050:     BVGetColumn(V,k,&v);
1051:     /* Form ut=[u;0] where u is the 1st column of U1 */
1052:     BVGetColumn(U1,k,&u);
1053:     VecZeroEntries(ut);
1054:     VecGetArrayRead(u,&carr);
1055:     VecGetArray(ut,&arr);
1056:     for (j=0; j<m; j++) arr[j] = carr[j];
1057:     VecRestoreArrayRead(u,&carr);
1058:     VecRestoreArray(ut,&arr);
1059:     /* Solve least squares problem */
1060:     KSPSolve(lanczos->ksp,ut,x);
1061:     MatMult(lanczos->Z,x,v);
1062:     BVRestoreColumn(U1,k,&u);
1063:     BVRestoreColumn(V,k,&v);
1064:     DSGetArrayReal(svd->ds,DS_MAT_T,&alpha);
1065:     if (breakdown) BVOrthonormalizeColumn(V,k,PETSC_FALSE,alpha+k,breakdown);
1066:     else BVOrthonormalizeColumn(V,k,PETSC_TRUE,alpha+k,NULL);
1067:     DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha);
1068:   }
1069:   PetscFunctionReturn(0);
1070: }

1072: /* solve generalized problem with joint lower-upper bidiagonalization */
1073: PetscErrorCode SVDSolve_TRLanczosGLower(SVD svd,BV U1,BV U2,BV V)
1074: {
1075:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
1076:   PetscReal      *alpha,*beta,*alphah,*betah,normr;
1077:   PetscScalar    *w;
1078:   PetscInt       i,k,l,nv,ld;
1079:   Mat            U,Vmat,X;
1080:   PetscBool      breakdown=PETSC_FALSE;

1082:   DSGetLeadingDimension(svd->ds,&ld);
1083:   PetscMalloc1(ld,&w);
1084:   normr = (svd->conv==SVD_CONV_ABS)? PetscMax(svd->nrma,svd->nrmb): 1.0;

1086:   /* normalize start vector */
1087:   if (!svd->ninil) BVSetRandomColumn(U2,0);
1088:   SVDInitialVectorGLower(svd,V,U1,U2,0,NULL);

1090:   l = 0;
1091:   while (svd->reason == SVD_CONVERGED_ITERATING) {
1092:     svd->its++;

1094:     /* inner loop */
1095:     nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
1096:     DSGetArrayReal(svd->ds,DS_MAT_T,&alpha);
1097:     DSGetArrayReal(svd->ds,DS_MAT_D,&alphah);
1098:     beta = alpha + ld;
1099:     betah = alpha + 2*ld;
1100:     SVDTwoSideLanczosGLower(svd,alpha,beta,alphah,betah,lanczos->Z,U1,U2,V,lanczos->ksp,svd->nconv+l,&nv,&breakdown);
1101:     DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha);
1102:     DSRestoreArrayReal(svd->ds,DS_MAT_D,&alphah);
1103:     BVSetActiveColumns(V,svd->nconv,nv);
1104:     BVSetActiveColumns(U1,svd->nconv,nv+1);
1105:     BVSetActiveColumns(U2,svd->nconv,nv);

1107:     /* solve projected problem */
1108:     DSSetDimensions(svd->ds,nv+1,svd->nconv,svd->nconv+l);
1109:     DSGSVDSetDimensions(svd->ds,nv,nv);
1110:     DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE);
1111:     DSSolve(svd->ds,w,NULL);
1112:     DSSort(svd->ds,w,NULL,NULL,NULL,NULL);
1113:     DSUpdateExtraRow(svd->ds);
1114:     DSSynchronize(svd->ds,w,NULL);
1115:     for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);

1117:     /* check convergence */
1118:     SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,normr,&k);
1119:     (*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx);

1121:     /* update l */
1122:     if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
1123:     else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
1124:     if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
1125:     if (l) PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l);

1127:     if (svd->reason == SVD_CONVERGED_ITERATING) {
1128:       if (PetscUnlikely(breakdown || k==nv)) {
1129:         /* Start a new bidiagonalization */
1130:         PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its);
1131:         if (k<svd->nsv) {
1132:           BVSetRandomColumn(U2,k);
1133:           SVDInitialVectorGLower(svd,V,U1,U2,k,&breakdown);
1134:           if (breakdown) {
1135:             svd->reason = SVD_DIVERGED_BREAKDOWN;
1136:             PetscInfo(svd,"Unable to generate more start vectors\n");
1137:           }
1138:         }
1139:       } else DSTruncate(svd->ds,k+l,PETSC_FALSE);
1140:     }

1142:     /* compute converged singular vectors and restart vectors */
1143:     DSGetMat(svd->ds,DS_MAT_X,&X);
1144:     BVMultInPlace(V,X,svd->nconv,k+l);
1145:     MatDestroy(&X);
1146:     DSGetMat(svd->ds,DS_MAT_U,&U);
1147:     BVMultInPlace(U1,U,svd->nconv,k+l+1);
1148:     MatDestroy(&U);
1149:     DSGetMat(svd->ds,DS_MAT_V,&Vmat);
1150:     BVMultInPlace(U2,Vmat,svd->nconv,k+l);
1151:     MatDestroy(&Vmat);

1153:     /* copy the last vector to be the next initial vector */
1154:     if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) BVCopyColumn(V,nv,k+l);

1156:     svd->nconv = k;
1157:     SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv);
1158:   }

1160:   PetscFree(w);
1161:   PetscFunctionReturn(0);
1162: }

1164: PetscErrorCode SVDSolve_TRLanczos_GSVD(SVD svd)
1165: {
1166:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
1167:   PetscInt       k,m,p;
1168:   PetscBool      convchg=PETSC_FALSE;
1169:   BV             U1,U2;
1170:   BVType         type;
1171:   Mat            U,V;

1173:   PetscCitationsRegister(citationg,&citedg);

1175:   if (svd->converged==SVDConvergedNorm) {  /* override temporarily since computed residual is already relative to the norms */
1176:     svd->converged = SVDConvergedAbsolute;
1177:     convchg = PETSC_TRUE;
1178:   }
1179:   MatGetLocalSize(svd->A,&m,NULL);
1180:   MatGetLocalSize(svd->B,&p,NULL);

1182:   /* Create BV for U1 */
1183:   BVCreate(PetscObjectComm((PetscObject)svd),&U1);
1184:   PetscLogObjectParent((PetscObject)svd,(PetscObject)U1);
1185:   BVGetType(svd->U,&type);
1186:   BVSetType(U1,type);
1187:   BVGetSizes(svd->U,NULL,NULL,&k);
1188:   BVSetSizes(U1,m,PETSC_DECIDE,k);

1190:   /* Create BV for U2 */
1191:   BVCreate(PetscObjectComm((PetscObject)svd),&U2);
1192:   PetscLogObjectParent((PetscObject)svd,(PetscObject)U2);
1193:   BVSetType(U2,type);
1194:   BVSetSizes(U2,p,PETSC_DECIDE,k);

1196:   /* Copy initial vectors from svd->U to U1 and U2 */
1197:   if (svd->ninil) {
1198:     Vec u, uh, nest, aux[2];
1199:     BVGetColumn(U1,0,&u);
1200:     BVGetColumn(U2,0,&uh);
1201:     aux[0] = u;
1202:     aux[1] = uh;
1203:     VecCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,aux,&nest);
1204:     BVCopyVec(svd->U,0,nest);
1205:     BVRestoreColumn(U1,0,&u);
1206:     BVRestoreColumn(U2,0,&uh);
1207:     VecDestroy(&nest);
1208:   }

1210:   switch (lanczos->bidiag) {
1211:     case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1212:       SVDSolve_TRLanczosGSingle(svd,U1,svd->U);
1213:       break;
1214:     case SVD_TRLANCZOS_GBIDIAG_UPPER:
1215:       SVDSolve_TRLanczosGUpper(svd,U1,U2,svd->U);
1216:       break;
1217:     case SVD_TRLANCZOS_GBIDIAG_LOWER:
1218:       SVDSolve_TRLanczosGLower(svd,U1,U2,svd->U);
1219:       break;
1220:   }

1222:   /* Compute converged right singular vectors */
1223:   BVSetActiveColumns(svd->U,0,svd->nconv);
1224:   BVSetActiveColumns(svd->V,0,svd->nconv);
1225:   BVGetMat(svd->U,&U);
1226:   BVGetMat(svd->V,&V);
1227:   KSPMatSolve(lanczos->ksp,U,V);
1228:   BVRestoreMat(svd->U,&U);
1229:   BVRestoreMat(svd->V,&V);

1231:   /* Finish computing left singular vectors and move them to its place */
1232:   switch (lanczos->bidiag) {
1233:     case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1234:       SVDLeftSingularVectors_Single(svd,U1,U2);
1235:       break;
1236:     case SVD_TRLANCZOS_GBIDIAG_UPPER:
1237:     case SVD_TRLANCZOS_GBIDIAG_LOWER:
1238:       SVDLeftSingularVectors(svd,U1,U2);
1239:       break;
1240:   }

1242:   BVDestroy(&U2);
1243:   BVDestroy(&U1);
1244:   DSTruncate(svd->ds,svd->nconv,PETSC_TRUE);
1245:   if (convchg) svd->converged = SVDConvergedNorm;
1246:   PetscFunctionReturn(0);
1247: }

1249: PetscErrorCode SVDSetFromOptions_TRLanczos(PetscOptionItems *PetscOptionsObject,SVD svd)
1250: {
1251:   SVD_TRLANCZOS       *lanczos = (SVD_TRLANCZOS*)svd->data;
1252:   PetscBool           flg,val,lock;
1253:   PetscReal           keep;
1254:   SVDTRLanczosGBidiag bidiag;

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

1258:     PetscOptionsBool("-svd_trlanczos_oneside","Use one-side reorthogonalization","SVDTRLanczosSetOneSide",lanczos->oneside,&val,&flg);
1259:     if (flg) SVDTRLanczosSetOneSide(svd,val);

1261:     PetscOptionsReal("-svd_trlanczos_restart","Proportion of vectors kept after restart","SVDTRLanczosSetRestart",0.5,&keep,&flg);
1262:     if (flg) SVDTRLanczosSetRestart(svd,keep);

1264:     PetscOptionsBool("-svd_trlanczos_locking","Choose between locking and non-locking variants","SVDTRLanczosSetLocking",PETSC_TRUE,&lock,&flg);
1265:     if (flg) SVDTRLanczosSetLocking(svd,lock);

1267:     PetscOptionsEnum("-svd_trlanczos_gbidiag","Bidiagonalization choice for Generalized Problem","SVDTRLanczosSetGBidiag",SVDTRLanczosGBidiags,(PetscEnum)lanczos->bidiag,(PetscEnum*)&bidiag,&flg);
1268:     if (flg) SVDTRLanczosSetGBidiag(svd,bidiag);

1270:     PetscOptionsBool("-svd_trlanczos_explicitmatrix","Build explicit matrix for KSP solver","SVDTRLanczosSetExplicitMatrix",lanczos->explicitmatrix,&val,&flg);
1271:     if (flg) SVDTRLanczosSetExplicitMatrix(svd,val);

1273:   PetscOptionsTail();

1275:   if (svd->OPb) {
1276:     if (!lanczos->ksp) SVDTRLanczosGetKSP(svd,&lanczos->ksp);
1277:     KSPSetFromOptions(lanczos->ksp);
1278:   }
1279:   PetscFunctionReturn(0);
1280: }

1282: static PetscErrorCode SVDTRLanczosSetOneSide_TRLanczos(SVD svd,PetscBool oneside)
1283: {
1284:   SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;

1286:   if (lanczos->oneside != oneside) {
1287:     lanczos->oneside = oneside;
1288:     svd->state = SVD_STATE_INITIAL;
1289:   }
1290:   PetscFunctionReturn(0);
1291: }

1293: /*@
1294:    SVDTRLanczosSetOneSide - Indicate if the variant of the Lanczos method
1295:    to be used is one-sided or two-sided.

1297:    Logically Collective on svd

1299:    Input Parameters:
1300: +  svd     - singular value solver
1301: -  oneside - boolean flag indicating if the method is one-sided or not

1303:    Options Database Key:
1304: .  -svd_trlanczos_oneside <boolean> - Indicates the boolean flag

1306:    Note:
1307:    By default, a two-sided variant is selected, which is sometimes slightly
1308:    more robust. However, the one-sided variant is faster because it avoids
1309:    the orthogonalization associated to left singular vectors.

1311:    Level: advanced

1313: .seealso: SVDLanczosSetOneSide()
1314: @*/
1315: PetscErrorCode SVDTRLanczosSetOneSide(SVD svd,PetscBool oneside)
1316: {
1319:   PetscTryMethod(svd,"SVDTRLanczosSetOneSide_C",(SVD,PetscBool),(svd,oneside));
1320:   PetscFunctionReturn(0);
1321: }

1323: static PetscErrorCode SVDTRLanczosGetOneSide_TRLanczos(SVD svd,PetscBool *oneside)
1324: {
1325:   SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;

1327:   *oneside = lanczos->oneside;
1328:   PetscFunctionReturn(0);
1329: }

1331: /*@
1332:    SVDTRLanczosGetOneSide - Gets if the variant of the Lanczos method
1333:    to be used is one-sided or two-sided.

1335:    Not Collective

1337:    Input Parameters:
1338: .  svd     - singular value solver

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

1343:    Level: advanced

1345: .seealso: SVDTRLanczosSetOneSide()
1346: @*/
1347: PetscErrorCode SVDTRLanczosGetOneSide(SVD svd,PetscBool *oneside)
1348: {
1351:   PetscUseMethod(svd,"SVDTRLanczosGetOneSide_C",(SVD,PetscBool*),(svd,oneside));
1352:   PetscFunctionReturn(0);
1353: }

1355: static PetscErrorCode SVDTRLanczosSetGBidiag_TRLanczos(SVD svd,SVDTRLanczosGBidiag bidiag)
1356: {
1357:   SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;

1359:   switch (bidiag) {
1360:     case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1361:     case SVD_TRLANCZOS_GBIDIAG_UPPER:
1362:     case SVD_TRLANCZOS_GBIDIAG_LOWER:
1363:       if (lanczos->bidiag != bidiag) {
1364:         lanczos->bidiag = bidiag;
1365:         svd->state = SVD_STATE_INITIAL;
1366:       }
1367:       break;
1368:     default:
1369:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid bidiagonalization choice");
1370:   }
1371:   PetscFunctionReturn(0);
1372: }

1374: /*@
1375:    SVDTRLanczosSetGBidiag - Sets the bidiagonalization choice to use in
1376:    the GSVD TRLanczos solver.

1378:    Logically Collective on svd

1380:    Input Parameters:
1381: +  svd - the singular value solver
1382: -  bidiag - the bidiagonalization choice

1384:    Options Database Key:
1385: .  -svd_trlanczos_gbidiag - Sets the bidiagonalization choice (either 's' or 'juu'
1386:    or 'jlu')

1388:    Level: advanced

1390: .seealso: SVDTRLanczosGetGBidiag(), SVDTRLanczosGBidiag
1391: @*/
1392: PetscErrorCode SVDTRLanczosSetGBidiag(SVD svd,SVDTRLanczosGBidiag bidiag)
1393: {
1396:   PetscTryMethod(svd,"SVDTRLanczosSetGBidiag_C",(SVD,SVDTRLanczosGBidiag),(svd,bidiag));
1397:   PetscFunctionReturn(0);
1398: }

1400: static PetscErrorCode SVDTRLanczosGetGBidiag_TRLanczos(SVD svd,SVDTRLanczosGBidiag *bidiag)
1401: {
1402:   SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;

1404:   *bidiag = lanczos->bidiag;
1405:   PetscFunctionReturn(0);
1406: }

1408: /*@
1409:    SVDTRLanczosGetGBidiag - Gets the bidiagonalization choice used in the GSVD
1410:    TRLanczos solver.

1412:    Not Collective

1414:    Input Parameter:
1415: .  svd - the singular value solver

1417:    Output Parameter:
1418: .  bidiag - the bidiagonalization choice

1420:    Level: advanced

1422: .seealso: SVDTRLanczosSetGBidiag(), SVDTRLanczosGBidiag
1423: @*/
1424: PetscErrorCode SVDTRLanczosGetGBidiag(SVD svd,SVDTRLanczosGBidiag *bidiag)
1425: {
1428:   PetscUseMethod(svd,"SVDTRLanczosGetGBidiag_C",(SVD,SVDTRLanczosGBidiag*),(svd,bidiag));
1429:   PetscFunctionReturn(0);
1430: }

1432: static PetscErrorCode SVDTRLanczosSetKSP_TRLanczos(SVD svd,KSP ksp)
1433: {
1434:   SVD_TRLANCZOS  *ctx = (SVD_TRLANCZOS*)svd->data;

1436:   PetscObjectReference((PetscObject)ksp);
1437:   KSPDestroy(&ctx->ksp);
1438:   ctx->ksp = ksp;
1439:   PetscLogObjectParent((PetscObject)svd,(PetscObject)ctx->ksp);
1440:   svd->state = SVD_STATE_INITIAL;
1441:   PetscFunctionReturn(0);
1442: }

1444: /*@
1445:    SVDTRLanczosSetKSP - Associate a linear solver object (KSP) to the SVD solver.

1447:    Collective on svd

1449:    Input Parameters:
1450: +  svd - SVD solver
1451: -  ksp - the linear solver object

1453:    Note:
1454:    Only used for the GSVD problem.

1456:    Level: advanced

1458: .seealso: SVDTRLanczosGetKSP()
1459: @*/
1460: PetscErrorCode SVDTRLanczosSetKSP(SVD svd,KSP ksp)
1461: {
1465:   PetscTryMethod(svd,"SVDTRLanczosSetKSP_C",(SVD,KSP),(svd,ksp));
1466:   PetscFunctionReturn(0);
1467: }

1469: static PetscErrorCode SVDTRLanczosGetKSP_TRLanczos(SVD svd,KSP *ksp)
1470: {
1471:   SVD_TRLANCZOS  *ctx = (SVD_TRLANCZOS*)svd->data;
1472:   PC             pc;

1474:   if (!ctx->ksp) {
1475:     /* Create linear solver */
1476:     KSPCreate(PetscObjectComm((PetscObject)svd),&ctx->ksp);
1477:     PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)svd,1);
1478:     KSPSetOptionsPrefix(ctx->ksp,((PetscObject)svd)->prefix);
1479:     KSPAppendOptionsPrefix(ctx->ksp,"svd_trlanczos_");
1480:     PetscLogObjectParent((PetscObject)svd,(PetscObject)ctx->ksp);
1481:     PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)svd)->options);
1482:     KSPSetType(ctx->ksp,KSPLSQR);
1483:     KSPGetPC(ctx->ksp,&pc);
1484:     PCSetType(pc,PCNONE);
1485:     KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE);
1486:     KSPSetTolerances(ctx->ksp,SlepcDefaultTol(svd->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1487:   }
1488:   *ksp = ctx->ksp;
1489:   PetscFunctionReturn(0);
1490: }

1492: /*@
1493:    SVDTRLanczosGetKSP - Retrieve the linear solver object (KSP) associated with
1494:    the SVD solver.

1496:    Not Collective

1498:    Input Parameter:
1499: .  svd - SVD solver

1501:    Output Parameter:
1502: .  ksp - the linear solver object

1504:    Level: advanced

1506: .seealso: SVDTRLanczosSetKSP()
1507: @*/
1508: PetscErrorCode SVDTRLanczosGetKSP(SVD svd,KSP *ksp)
1509: {
1512:   PetscUseMethod(svd,"SVDTRLanczosGetKSP_C",(SVD,KSP*),(svd,ksp));
1513:   PetscFunctionReturn(0);
1514: }

1516: static PetscErrorCode SVDTRLanczosSetRestart_TRLanczos(SVD svd,PetscReal keep)
1517: {
1518:   SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;

1520:   if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
1521:   else {
1523:     ctx->keep = keep;
1524:   }
1525:   PetscFunctionReturn(0);
1526: }

1528: /*@
1529:    SVDTRLanczosSetRestart - Sets the restart parameter for the thick-restart
1530:    Lanczos method, in particular the proportion of basis vectors that must be
1531:    kept after restart.

1533:    Logically Collective on svd

1535:    Input Parameters:
1536: +  svd  - the singular value solver
1537: -  keep - the number of vectors to be kept at restart

1539:    Options Database Key:
1540: .  -svd_trlanczos_restart - Sets the restart parameter

1542:    Notes:
1543:    Allowed values are in the range [0.1,0.9]. The default is 0.5.

1545:    Level: advanced

1547: .seealso: SVDTRLanczosGetRestart()
1548: @*/
1549: PetscErrorCode SVDTRLanczosSetRestart(SVD svd,PetscReal keep)
1550: {
1553:   PetscTryMethod(svd,"SVDTRLanczosSetRestart_C",(SVD,PetscReal),(svd,keep));
1554:   PetscFunctionReturn(0);
1555: }

1557: static PetscErrorCode SVDTRLanczosGetRestart_TRLanczos(SVD svd,PetscReal *keep)
1558: {
1559:   SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;

1561:   *keep = ctx->keep;
1562:   PetscFunctionReturn(0);
1563: }

1565: /*@
1566:    SVDTRLanczosGetRestart - Gets the restart parameter used in the thick-restart
1567:    Lanczos method.

1569:    Not Collective

1571:    Input Parameter:
1572: .  svd - the singular value solver

1574:    Output Parameter:
1575: .  keep - the restart parameter

1577:    Level: advanced

1579: .seealso: SVDTRLanczosSetRestart()
1580: @*/
1581: PetscErrorCode SVDTRLanczosGetRestart(SVD svd,PetscReal *keep)
1582: {
1585:   PetscUseMethod(svd,"SVDTRLanczosGetRestart_C",(SVD,PetscReal*),(svd,keep));
1586:   PetscFunctionReturn(0);
1587: }

1589: static PetscErrorCode SVDTRLanczosSetLocking_TRLanczos(SVD svd,PetscBool lock)
1590: {
1591:   SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;

1593:   ctx->lock = lock;
1594:   PetscFunctionReturn(0);
1595: }

1597: /*@
1598:    SVDTRLanczosSetLocking - Choose between locking and non-locking variants of
1599:    the thick-restart Lanczos method.

1601:    Logically Collective on svd

1603:    Input Parameters:
1604: +  svd  - the singular value solver
1605: -  lock - true if the locking variant must be selected

1607:    Options Database Key:
1608: .  -svd_trlanczos_locking - Sets the locking flag

1610:    Notes:
1611:    The default is to lock converged singular triplets when the method restarts.
1612:    This behaviour can be changed so that all directions are kept in the
1613:    working subspace even if already converged to working accuracy (the
1614:    non-locking variant).

1616:    Level: advanced

1618: .seealso: SVDTRLanczosGetLocking()
1619: @*/
1620: PetscErrorCode SVDTRLanczosSetLocking(SVD svd,PetscBool lock)
1621: {
1624:   PetscTryMethod(svd,"SVDTRLanczosSetLocking_C",(SVD,PetscBool),(svd,lock));
1625:   PetscFunctionReturn(0);
1626: }

1628: static PetscErrorCode SVDTRLanczosGetLocking_TRLanczos(SVD svd,PetscBool *lock)
1629: {
1630:   SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;

1632:   *lock = ctx->lock;
1633:   PetscFunctionReturn(0);
1634: }

1636: /*@
1637:    SVDTRLanczosGetLocking - Gets the locking flag used in the thick-restart
1638:    Lanczos method.

1640:    Not Collective

1642:    Input Parameter:
1643: .  svd - the singular value solver

1645:    Output Parameter:
1646: .  lock - the locking flag

1648:    Level: advanced

1650: .seealso: SVDTRLanczosSetLocking()
1651: @*/
1652: PetscErrorCode SVDTRLanczosGetLocking(SVD svd,PetscBool *lock)
1653: {
1656:   PetscUseMethod(svd,"SVDTRLanczosGetLocking_C",(SVD,PetscBool*),(svd,lock));
1657:   PetscFunctionReturn(0);
1658: }

1660: static PetscErrorCode SVDTRLanczosSetExplicitMatrix_TRLanczos(SVD svd,PetscBool explicitmat)
1661: {
1662:   SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;

1664:   if (lanczos->explicitmatrix != explicitmat) {
1665:     lanczos->explicitmatrix = explicitmat;
1666:     svd->state = SVD_STATE_INITIAL;
1667:   }
1668:   PetscFunctionReturn(0);
1669: }

1671: /*@
1672:    SVDTRLanczosSetExplicitMatrix - Indicate if the matrix Z=[A;B] must
1673:    be built explicitly.

1675:    Logically Collective on svd

1677:    Input Parameters:
1678: +  svd         - singular value solver
1679: -  explicitmat - Boolean flag indicating if Z=[A;B] is built explicitly

1681:    Options Database Key:
1682: .  -svd_trlanczos_explicitmatrix <boolean> - Indicates the boolean flag

1684:    Notes:
1685:    This option is relevant for the GSVD case only.
1686:    Z is the coefficient matrix of the KSP solver used internally.

1688:    Level: advanced

1690: .seealso: SVDTRLanczosGetExplicitMatrix()
1691: @*/
1692: PetscErrorCode SVDTRLanczosSetExplicitMatrix(SVD svd,PetscBool explicitmat)
1693: {
1696:   PetscTryMethod(svd,"SVDTRLanczosSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
1697:   PetscFunctionReturn(0);
1698: }

1700: static PetscErrorCode SVDTRLanczosGetExplicitMatrix_TRLanczos(SVD svd,PetscBool *explicitmat)
1701: {
1702:   SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;

1704:   *explicitmat = lanczos->explicitmatrix;
1705:   PetscFunctionReturn(0);
1706: }

1708: /*@
1709:    SVDTRLanczosGetExplicitMatrix - Returns the flag indicating if Z=[A;B] is built explicitly.

1711:    Not Collective

1713:    Input Parameter:
1714: .  svd  - singular value solver

1716:    Output Parameter:
1717: .  explicitmat - the mode flag

1719:    Level: advanced

1721: .seealso: SVDTRLanczosSetExplicitMatrix()
1722: @*/
1723: PetscErrorCode SVDTRLanczosGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
1724: {
1727:   PetscUseMethod(svd,"SVDTRLanczosGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
1728:   PetscFunctionReturn(0);
1729: }

1731: PetscErrorCode SVDReset_TRLanczos(SVD svd)
1732: {
1733:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;

1735:   if (svd->isgeneralized || (!svd->problem_type && svd->OPb)) {
1736:     KSPReset(lanczos->ksp);
1737:     MatDestroy(&lanczos->Z);
1738:   }
1739:   PetscFunctionReturn(0);
1740: }

1742: PetscErrorCode SVDDestroy_TRLanczos(SVD svd)
1743: {
1744:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;

1746:   if (svd->isgeneralized || (!svd->problem_type && svd->OPb)) KSPDestroy(&lanczos->ksp);
1747:   PetscFree(svd->data);
1748:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",NULL);
1749:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",NULL);
1750:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetGBidiag_C",NULL);
1751:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetGBidiag_C",NULL);
1752:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetKSP_C",NULL);
1753:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetKSP_C",NULL);
1754:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetRestart_C",NULL);
1755:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetRestart_C",NULL);
1756:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetLocking_C",NULL);
1757:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetLocking_C",NULL);
1758:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetExplicitMatrix_C",NULL);
1759:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetExplicitMatrix_C",NULL);
1760:   PetscFunctionReturn(0);
1761: }

1763: PetscErrorCode SVDView_TRLanczos(SVD svd,PetscViewer viewer)
1764: {
1765:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
1766:   PetscBool      isascii;

1768:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1769:   if (isascii) {
1770:     PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*lanczos->keep));
1771:     PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",lanczos->lock?"":"non-");
1772:     if (svd->isgeneralized) {
1773:       const char *bidiag="";

1775:       switch (lanczos->bidiag) {
1776:         case SVD_TRLANCZOS_GBIDIAG_SINGLE: bidiag = "single"; break;
1777:         case SVD_TRLANCZOS_GBIDIAG_UPPER:  bidiag = "joint upper-upper"; break;
1778:         case SVD_TRLANCZOS_GBIDIAG_LOWER:  bidiag = "joint lower-upper"; break;
1779:       }
1780:       PetscViewerASCIIPrintf(viewer,"  bidiagonalization choice: %s\n",bidiag);
1781:       PetscViewerASCIIPrintf(viewer,"  %s matrix\n",lanczos->explicitmatrix?"explicit":"implicit");
1782:       if (!lanczos->ksp) SVDTRLanczosGetKSP(svd,&lanczos->ksp);
1783:       PetscViewerASCIIPushTab(viewer);
1784:       KSPView(lanczos->ksp,viewer);
1785:       PetscViewerASCIIPopTab(viewer);
1786:     } else PetscViewerASCIIPrintf(viewer,"  %s-sided reorthogonalization\n",lanczos->oneside? "one": "two");
1787:   }
1788:   PetscFunctionReturn(0);
1789: }

1791: SLEPC_EXTERN PetscErrorCode SVDCreate_TRLanczos(SVD svd)
1792: {
1793:   SVD_TRLANCZOS  *ctx;

1795:   PetscNewLog(svd,&ctx);
1796:   svd->data = (void*)ctx;

1798:   ctx->lock   = PETSC_TRUE;
1799:   ctx->bidiag = SVD_TRLANCZOS_GBIDIAG_LOWER;

1801:   svd->ops->setup          = SVDSetUp_TRLanczos;
1802:   svd->ops->solve          = SVDSolve_TRLanczos;
1803:   svd->ops->solveg         = SVDSolve_TRLanczos_GSVD;
1804:   svd->ops->destroy        = SVDDestroy_TRLanczos;
1805:   svd->ops->reset          = SVDReset_TRLanczos;
1806:   svd->ops->setfromoptions = SVDSetFromOptions_TRLanczos;
1807:   svd->ops->view           = SVDView_TRLanczos;
1808:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",SVDTRLanczosSetOneSide_TRLanczos);
1809:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",SVDTRLanczosGetOneSide_TRLanczos);
1810:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetGBidiag_C",SVDTRLanczosSetGBidiag_TRLanczos);
1811:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetGBidiag_C",SVDTRLanczosGetGBidiag_TRLanczos);
1812:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetKSP_C",SVDTRLanczosSetKSP_TRLanczos);
1813:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetKSP_C",SVDTRLanczosGetKSP_TRLanczos);
1814:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetRestart_C",SVDTRLanczosSetRestart_TRLanczos);
1815:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetRestart_C",SVDTRLanczosGetRestart_TRLanczos);
1816:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetLocking_C",SVDTRLanczosSetLocking_TRLanczos);
1817:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetLocking_C",SVDTRLanczosGetLocking_TRLanczos);
1818:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetExplicitMatrix_C",SVDTRLanczosSetExplicitMatrix_TRLanczos);
1819:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetExplicitMatrix_C",SVDTRLanczosGetExplicitMatrix_TRLanczos);
1820:   PetscFunctionReturn(0);
1821: }