Actual source code: trlanczos.c

slepc-3.16.2 2022-02-01
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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} bidigonalization 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: {

 75:   PetscNew(zdata);
 76:   (*zdata)->A = svd->A;
 77:   (*zdata)->B = svd->B;
 78:   (*zdata)->AT = svd->AT;
 79:   (*zdata)->BT = svd->BT;
 80:   MatCreateVecsEmpty(svd->A,NULL,&(*zdata)->y1);
 81:   MatCreateVecsEmpty(svd->B,NULL,&(*zdata)->y2);
 82:   VecGetLocalSize((*zdata)->y1,&(*zdata)->m);
 83:   BVCreateVec(svd->U,&(*zdata)->y);
 84:   return(0);
 85: }

 87: static PetscErrorCode MatDestroy_Z(Mat Z)
 88: {
 90:   MatZData       *zdata;

 93:   MatShellGetContext(Z,&zdata);
 94:   VecDestroy(&zdata->y1);
 95:   VecDestroy(&zdata->y2);
 96:   VecDestroy(&zdata->y);
 97:   PetscFree(zdata);
 98:   return(0);
 99: }

101: static PetscErrorCode MatMult_Z(Mat Z,Vec x,Vec y)
102: {
103:   MatZData       *zdata;
104:   PetscScalar    *y_elems;

108:   MatShellGetContext(Z,&zdata);
109:   VecGetArray(y,&y_elems);
110:   VecPlaceArray(zdata->y1,y_elems);
111:   VecPlaceArray(zdata->y2,y_elems+zdata->m);

113:   MatMult(zdata->A,x,zdata->y1);
114:   MatMult(zdata->B,x,zdata->y2);

116:   VecResetArray(zdata->y1);
117:   VecResetArray(zdata->y2);
118:   VecRestoreArray(y,&y_elems);
119:   return(0);
120: }

122: static PetscErrorCode MatMultTranspose_Z(Mat Z,Vec y,Vec x)
123: {
124:   MatZData          *zdata;
125:   const PetscScalar *y_elems;
126:   PetscErrorCode    ierr;

129:   MatShellGetContext(Z,&zdata);
130:   VecGetArrayRead(y,&y_elems);
131:   VecPlaceArray(zdata->y1,y_elems);
132:   VecPlaceArray(zdata->y2,y_elems+zdata->m);

134:   MatMult(zdata->AT,zdata->y1,x);
135:   MatMultAdd(zdata->BT,zdata->y2,x,x);

137:   VecResetArray(zdata->y1);
138:   VecResetArray(zdata->y2);
139:   VecRestoreArrayRead(y,&y_elems);
140:   return(0);
141: }

143: static PetscErrorCode MatCreateVecs_Z(Mat Z,Vec *right,Vec *left)
144: {
145:   MatZData       *zdata;

149:   MatShellGetContext(Z,&zdata);
150:   if (right) {
151:     MatCreateVecs(zdata->A,right,NULL);
152:   }
153:   if (left) {
154:     VecDuplicate(zdata->y,left);
155:   }
156:   return(0);
157: }

159: PetscErrorCode SVDSetUp_TRLanczos(SVD svd)
160: {
162:   PetscInt       N,m,n,p;
163:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
164:   DSType         dstype;
165:   MatZData       *zdata;
166:   Mat            mats[2],normal;
167:   MatType        Atype;
168:   PetscBool      sametype;

171:   MatGetSize(svd->A,NULL,&N);
172:   SVDSetDimensions_Default(svd);
173:   if (svd->ncv>svd->nsv+svd->mpd) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nsv+mpd");
174:   if (!lanczos->lock && svd->mpd<svd->ncv) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
175:   if (svd->max_it==PETSC_DEFAULT) svd->max_it = PetscMax(N/svd->ncv,100);
176:   if (!lanczos->keep) lanczos->keep = 0.5;
177:   svd->leftbasis = PETSC_TRUE;
178:   SVDAllocateSolution(svd,1);
179:   dstype = DSSVD;
180:   if (svd->isgeneralized) {
181:     if (lanczos->bidiag==SVD_TRLANCZOS_GBIDIAG_UPPER || lanczos->bidiag==SVD_TRLANCZOS_GBIDIAG_LOWER) dstype = DSGSVD;
182:     SVDSetWorkVecs(svd,1,1);

184:     /* Create the matrix Z=[A;B] */
185:     MatDestroy(&lanczos->Z);
186:     MatGetLocalSize(svd->A,&m,&n);
187:     MatGetLocalSize(svd->B,&p,NULL);
188:     if (lanczos->explicitmatrix) {
189:       mats[0] = svd->A;
190:       mats[1] = svd->B;
191:       MatCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,1,NULL,mats,&lanczos->Z);
192:       MatGetType(svd->A,&Atype);
193:       PetscObjectTypeCompare((PetscObject)svd->B,Atype,&sametype);
194:       if (!sametype) Atype = MATAIJ;
195:       MatConvert(lanczos->Z,Atype,MAT_INPLACE_MATRIX,&lanczos->Z);
196:     } else {
197:       MatZCreateContext(svd,&zdata);
198:       MatCreateShell(PetscObjectComm((PetscObject)svd),m+p,n,PETSC_DECIDE,PETSC_DECIDE,zdata,&lanczos->Z);
199:       MatShellSetOperation(lanczos->Z,MATOP_MULT,(void(*)(void))MatMult_Z);
200: #if defined(PETSC_USE_COMPLEX)
201:       MatShellSetOperation(lanczos->Z,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultTranspose_Z);
202: #else
203:       MatShellSetOperation(lanczos->Z,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Z);
204: #endif
205:       MatShellSetOperation(lanczos->Z,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_Z);
206:       MatShellSetOperation(lanczos->Z,MATOP_DESTROY,(void(*)(void))MatDestroy_Z);
207:     }
208:     PetscLogObjectParent((PetscObject)svd,(PetscObject)lanczos->Z);

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

213:     if (!lanczos->ksp) { SVDTRLanczosGetKSP(svd,&lanczos->ksp); }
214:     KSPSetOperators(lanczos->ksp,lanczos->Z,normal);
215:     KSPSetUp(lanczos->ksp);
216:     MatDestroy(&normal);

218:     if (lanczos->oneside) { PetscInfo(svd,"Warning: one-side option is ignored in GSVD\n"); }
219:   }
220:   DSSetType(svd->ds,dstype);
221:   DSSetCompact(svd->ds,PETSC_TRUE);
222:   DSSetExtraRow(svd->ds,PETSC_TRUE);
223:   DSAllocate(svd->ds,svd->ncv+1);
224:   return(0);
225: }

227: static PetscErrorCode SVDOneSideTRLanczosMGS(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,BV U,PetscInt nconv,PetscInt l,PetscInt n,PetscScalar* work)
228: {
230:   PetscReal      a,b;
231:   PetscInt       i,k=nconv+l;
232:   Vec            ui,ui1,vi;

235:   BVGetColumn(V,k,&vi);
236:   BVGetColumn(U,k,&ui);
237:   MatMult(svd->A,vi,ui);
238:   BVRestoreColumn(V,k,&vi);
239:   BVRestoreColumn(U,k,&ui);
240:   if (l>0) {
241:     BVSetActiveColumns(U,nconv,n);
242:     for (i=0;i<l;i++) work[i]=beta[i+nconv];
243:     BVMultColumn(U,-1.0,1.0,k,work);
244:   }
245:   BVNormColumn(U,k,NORM_2,&a);
246:   BVScaleColumn(U,k,1.0/a);
247:   alpha[k] = a;

249:   for (i=k+1;i<n;i++) {
250:     BVGetColumn(V,i,&vi);
251:     BVGetColumn(U,i-1,&ui1);
252:     MatMult(svd->AT,ui1,vi);
253:     BVRestoreColumn(V,i,&vi);
254:     BVRestoreColumn(U,i-1,&ui1);
255:     BVOrthonormalizeColumn(V,i,PETSC_FALSE,&b,NULL);
256:     beta[i-1] = b;

258:     BVGetColumn(V,i,&vi);
259:     BVGetColumn(U,i,&ui);
260:     MatMult(svd->A,vi,ui);
261:     BVRestoreColumn(V,i,&vi);
262:     BVGetColumn(U,i-1,&ui1);
263:     VecAXPY(ui,-b,ui1);
264:     BVRestoreColumn(U,i-1,&ui1);
265:     BVRestoreColumn(U,i,&ui);
266:     BVNormColumn(U,i,NORM_2,&a);
267:     BVScaleColumn(U,i,1.0/a);
268:     alpha[i] = a;
269:   }

271:   BVGetColumn(V,n,&vi);
272:   BVGetColumn(U,n-1,&ui1);
273:   MatMult(svd->AT,ui1,vi);
274:   BVRestoreColumn(V,n,&vi);
275:   BVRestoreColumn(U,n-1,&ui1);
276:   BVOrthogonalizeColumn(V,n,NULL,&b,NULL);
277:   beta[n-1] = b;
278:   return(0);
279: }

281: /*
282:   Custom CGS orthogonalization, preprocess after first orthogonalization
283: */
284: static PetscErrorCode SVDOrthogonalizeCGS(BV V,PetscInt i,PetscScalar* h,PetscReal a,BVOrthogRefineType refine,PetscReal eta,PetscReal *norm)
285: {
287:   PetscReal      sum,onorm;
288:   PetscScalar    dot;
289:   PetscInt       j;

292:   switch (refine) {
293:   case BV_ORTHOG_REFINE_NEVER:
294:     BVNormColumn(V,i,NORM_2,norm);
295:     break;
296:   case BV_ORTHOG_REFINE_ALWAYS:
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:     break;
302:   case BV_ORTHOG_REFINE_IFNEEDED:
303:     dot = h[i];
304:     onorm = PetscSqrtReal(PetscRealPart(dot)) / a;
305:     sum = 0.0;
306:     for (j=0;j<i;j++) {
307:       sum += PetscRealPart(h[j] * PetscConj(h[j]));
308:     }
309:     *norm = PetscRealPart(dot)/(a*a) - sum;
310:     if (*norm>0.0) *norm = PetscSqrtReal(*norm);
311:     else {
312:       BVNormColumn(V,i,NORM_2,norm);
313:     }
314:     if (*norm < eta*onorm) {
315:       BVSetActiveColumns(V,0,i);
316:       BVDotColumn(V,i,h);
317:       BVMultColumn(V,-1.0,1.0,i,h);
318:       BVNormColumn(V,i,NORM_2,norm);
319:     }
320:     break;
321:   }
322:   return(0);
323: }

325: static PetscErrorCode SVDOneSideTRLanczosCGS(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,BV U,PetscInt nconv,PetscInt l,PetscInt n,PetscScalar* work)
326: {
327:   PetscErrorCode     ierr;
328:   PetscReal          a,b,eta;
329:   PetscInt           i,j,k=nconv+l;
330:   Vec                ui,ui1,vi;
331:   BVOrthogRefineType refine;

334:   BVGetColumn(V,k,&vi);
335:   BVGetColumn(U,k,&ui);
336:   MatMult(svd->A,vi,ui);
337:   BVRestoreColumn(V,k,&vi);
338:   BVRestoreColumn(U,k,&ui);
339:   if (l>0) {
340:     BVSetActiveColumns(U,nconv,n);
341:     for (i=0;i<l;i++) work[i]=beta[i+nconv];
342:     BVMultColumn(U,-1.0,1.0,k,work);
343:   }
344:   BVGetOrthogonalization(V,NULL,&refine,&eta,NULL);

346:   for (i=k+1;i<n;i++) {
347:     BVGetColumn(V,i,&vi);
348:     BVGetColumn(U,i-1,&ui1);
349:     MatMult(svd->AT,ui1,vi);
350:     BVRestoreColumn(V,i,&vi);
351:     BVRestoreColumn(U,i-1,&ui1);
352:     BVNormColumnBegin(U,i-1,NORM_2,&a);
353:     if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
354:       BVSetActiveColumns(V,0,i+1);
355:       BVGetColumn(V,i,&vi);
356:       BVDotVecBegin(V,vi,work);
357:     } else {
358:       BVSetActiveColumns(V,0,i);
359:       BVDotColumnBegin(V,i,work);
360:     }
361:     BVNormColumnEnd(U,i-1,NORM_2,&a);
362:     if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
363:       BVDotVecEnd(V,vi,work);
364:       BVRestoreColumn(V,i,&vi);
365:       BVSetActiveColumns(V,0,i);
366:     } else {
367:       BVDotColumnEnd(V,i,work);
368:     }

370:     BVScaleColumn(U,i-1,1.0/a);
371:     for (j=0;j<i;j++) work[j] = work[j] / a;
372:     BVMultColumn(V,-1.0,1.0/a,i,work);
373:     SVDOrthogonalizeCGS(V,i,work,a,refine,eta,&b);
374:     BVScaleColumn(V,i,1.0/b);
375:     if (PetscAbsReal(b)<10*PETSC_MACHINE_EPSILON) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_PLIB,"Recurrence generated a zero vector; use a two-sided variant");

377:     BVGetColumn(V,i,&vi);
378:     BVGetColumn(U,i,&ui);
379:     BVGetColumn(U,i-1,&ui1);
380:     MatMult(svd->A,vi,ui);
381:     VecAXPY(ui,-b,ui1);
382:     BVRestoreColumn(V,i,&vi);
383:     BVRestoreColumn(U,i,&ui);
384:     BVRestoreColumn(U,i-1,&ui1);

386:     alpha[i-1] = a;
387:     beta[i-1] = b;
388:   }

390:   BVGetColumn(V,n,&vi);
391:   BVGetColumn(U,n-1,&ui1);
392:   MatMult(svd->AT,ui1,vi);
393:   BVRestoreColumn(V,n,&vi);
394:   BVRestoreColumn(U,n-1,&ui1);

396:   BVNormColumnBegin(svd->U,n-1,NORM_2,&a);
397:   if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
398:     BVSetActiveColumns(V,0,n+1);
399:     BVGetColumn(V,n,&vi);
400:     BVDotVecBegin(V,vi,work);
401:   } else {
402:     BVSetActiveColumns(V,0,n);
403:     BVDotColumnBegin(V,n,work);
404:   }
405:   BVNormColumnEnd(svd->U,n-1,NORM_2,&a);
406:   if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
407:     BVDotVecEnd(V,vi,work);
408:     BVRestoreColumn(V,n,&vi);
409:   } else {
410:     BVDotColumnEnd(V,n,work);
411:   }

413:   BVScaleColumn(U,n-1,1.0/a);
414:   for (j=0;j<n;j++) work[j] = work[j] / a;
415:   BVMultColumn(V,-1.0,1.0/a,n,work);
416:   SVDOrthogonalizeCGS(V,n,work,a,refine,eta,&b);
417:   BVSetActiveColumns(V,nconv,n);
418:   alpha[n-1] = a;
419:   beta[n-1] = b;
420:   return(0);
421: }

423: PetscErrorCode SVDSolve_TRLanczos(SVD svd)
424: {
426:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
427:   PetscReal      *alpha,*beta;
428:   PetscScalar    *swork=NULL,*w;
429:   PetscInt       i,k,l,nv,ld;
430:   Mat            U,V;
431:   PetscBool      breakdown=PETSC_FALSE;
432:   BVOrthogType   orthog;

435:   PetscCitationsRegister(citation,&cited);
436:   /* allocate working space */
437:   DSGetLeadingDimension(svd->ds,&ld);
438:   BVGetOrthogonalization(svd->V,&orthog,NULL,NULL,NULL);
439:   PetscMalloc1(ld,&w);
440:   if (lanczos->oneside) {
441:     PetscMalloc1(svd->ncv+1,&swork);
442:   }

444:   /* normalize start vector */
445:   if (!svd->nini) {
446:     BVSetRandomColumn(svd->V,0);
447:     BVOrthonormalizeColumn(svd->V,0,PETSC_TRUE,NULL,NULL);
448:   }

450:   l = 0;
451:   while (svd->reason == SVD_CONVERGED_ITERATING) {
452:     svd->its++;

454:     /* inner loop */
455:     nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
456:     DSGetArrayReal(svd->ds,DS_MAT_T,&alpha);
457:     beta = alpha + ld;
458:     if (lanczos->oneside) {
459:       if (orthog == BV_ORTHOG_MGS) {
460:         SVDOneSideTRLanczosMGS(svd,alpha,beta,svd->V,svd->U,svd->nconv,l,nv,swork);
461:       } else {
462:         SVDOneSideTRLanczosCGS(svd,alpha,beta,svd->V,svd->U,svd->nconv,l,nv,swork);
463:       }
464:     } else {
465:       SVDTwoSideLanczos(svd,alpha,beta,svd->V,svd->U,svd->nconv+l,&nv,&breakdown);
466:     }
467:     DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha);
468:     BVScaleColumn(svd->V,nv,1.0/beta[nv-1]);
469:     BVSetActiveColumns(svd->V,svd->nconv,nv);
470:     BVSetActiveColumns(svd->U,svd->nconv,nv);

472:     /* solve projected problem */
473:     DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l);
474:     DSSVDSetDimensions(svd->ds,nv);
475:     if (l==0) {
476:       DSSetState(svd->ds,DS_STATE_INTERMEDIATE);
477:     } else {
478:       DSSetState(svd->ds,DS_STATE_RAW);
479:     }
480:     DSSolve(svd->ds,w,NULL);
481:     DSSort(svd->ds,w,NULL,NULL,NULL,NULL);
482:     DSUpdateExtraRow(svd->ds);
483:     DSSynchronize(svd->ds,w,NULL);
484:     for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);

486:     /* check convergence */
487:     SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,&k);
488:     (*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx);

490:     /* update l */
491:     if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
492:     else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
493:     if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
494:     if (l) { PetscInfo1(svd,"Preparing to restart keeping l=%D vectors\n",l); }

496:     if (svd->reason == SVD_CONVERGED_ITERATING) {
497:       if (breakdown || k==nv) {
498:         /* Start a new bidiagonalization */
499:         PetscInfo1(svd,"Breakdown in bidiagonalization (it=%D)\n",svd->its);
500:         if (k<svd->nsv) {
501:           BVSetRandomColumn(svd->V,k);
502:           BVOrthonormalizeColumn(svd->V,k,PETSC_FALSE,NULL,&breakdown);
503:           if (breakdown) {
504:             svd->reason = SVD_DIVERGED_BREAKDOWN;
505:             PetscInfo(svd,"Unable to generate more start vectors\n");
506:           }
507:         }
508:       } else {
509:         DSTruncate(svd->ds,k+l,PETSC_FALSE);
510:       }
511:     }

513:     /* compute converged singular vectors and restart vectors */
514:     DSGetMat(svd->ds,DS_MAT_V,&V);
515:     BVMultInPlace(svd->V,V,svd->nconv,k+l);
516:     MatDestroy(&V);
517:     DSGetMat(svd->ds,DS_MAT_U,&U);
518:     BVMultInPlace(svd->U,U,svd->nconv,k+l);
519:     MatDestroy(&U);

521:     /* copy the last vector to be the next initial vector */
522:     if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) {
523:       BVCopyColumn(svd->V,nv,k+l);
524:     }

526:     svd->nconv = k;
527:     SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv);
528:   }

530:   /* orthonormalize U columns in one side method */
531:   if (lanczos->oneside) {
532:     for (i=0;i<svd->nconv;i++) {
533:       BVOrthonormalizeColumn(svd->U,i,PETSC_FALSE,NULL,NULL);
534:     }
535:   }

537:   /* free working space */
538:   PetscFree(w);
539:   if (swork) { PetscFree(swork); }
540:   DSTruncate(svd->ds,svd->nconv,PETSC_TRUE);
541:   return(0);
542: }

544: static PetscErrorCode SVDTwoSideLanczosGSingle(SVD svd,PetscReal *alpha,PetscReal *beta,Mat Z,BV V,BV U,KSP ksp,PetscInt k,PetscInt *n,PetscBool *breakdown)
545: {
546:   PetscErrorCode    ierr;
547:   PetscInt          i,j,m;
548:   const PetscScalar *carr;
549:   PetscScalar       *arr;
550:   Vec               u,v,ut=svd->workl[0],x=svd->workr[0],v1;
551:   PetscBool         lindep=PETSC_FALSE;

554:   MatCreateVecsEmpty(svd->A,NULL,&v1);
555:   BVGetColumn(V,k,&v);
556:   BVGetColumn(U,k,&u);

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

567:   /* Solve least squares problem */
568:   KSPSolve(ksp,ut,x);

570:   MatMult(Z,x,v);

572:   BVRestoreColumn(U,k,&u);
573:   BVRestoreColumn(V,k,&v);
574:   BVOrthonormalizeColumn(V,k,PETSC_FALSE,alpha+k,&lindep);
575:   if (lindep) {
576:     *n = k;
577:     if (breakdown) *breakdown = lindep;
578:     return(0);
579:   }

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

583:     /* Compute vector i of BV U */
584:     BVGetColumn(V,i-1,&v);
585:     VecGetArray(v,&arr);
586:     VecPlaceArray(v1,arr);
587:     VecRestoreArray(v,&arr);
588:     BVRestoreColumn(V,i-1,&v);
589:     BVInsertVec(U,i,v1);
590:     VecResetArray(v1);
591:     BVOrthonormalizeColumn(U,i,PETSC_FALSE,beta+i-1,&lindep);
592:     if (lindep) {
593:       *n = i;
594:       break;
595:     }

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

599:     BVGetColumn(V,i,&v);
600:     BVGetColumn(U,i,&u);

602:     /* Form ut=[u;0] */
603:     VecGetArrayRead(u,&carr);
604:     VecGetArray(ut,&arr);
605:     for (j=0; j<m; j++) arr[j] = carr[j];
606:     VecRestoreArrayRead(u,&carr);
607:     VecRestoreArray(ut,&arr);

609:     /* Solve least squares problem */
610:     KSPSolve(ksp,ut,x);

612:     MatMult(Z,x,v);

614:     BVRestoreColumn(U,i,&u);
615:     BVRestoreColumn(V,i,&v);
616:     BVOrthonormalizeColumn(V,i,PETSC_FALSE,alpha+i,&lindep);
617:     if (lindep) {
618:       *n = i;
619:       break;
620:     }
621:   }

623:   /* Compute vector n of BV U */
624:   if (!lindep) {
625:     BVGetColumn(V,*n-1,&v);
626:     VecGetArray(v,&arr);
627:     VecPlaceArray(v1,arr);
628:     VecRestoreArray(v,&arr);
629:     BVRestoreColumn(V,*n-1,&v);
630:     BVInsertVec(U,*n,v1);
631:     VecResetArray(v1);
632:     BVOrthonormalizeColumn(U,*n,PETSC_FALSE,beta+*n-1,&lindep);
633:   }
634:   if (breakdown) *breakdown = lindep;
635:   VecDestroy(&v1);
636:   return(0);
637: }

639: /* solve generalized problem with single bidiagonalization of Q_A */
640: PetscErrorCode SVDSolve_TRLanczosGSingle(SVD svd,BV U1,BV V)
641: {
643:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
644:   PetscReal      *alpha,*beta;
645:   PetscScalar    *w;
646:   PetscInt       i,k,l,nv,ld;
647:   Mat            U,VV;
648:   PetscBool      breakdown=PETSC_FALSE;

651:   DSGetLeadingDimension(svd->ds,&ld);
652:   PetscMalloc1(ld,&w);

654:   /* normalize start vector */
655:   if (!svd->nini) {
656:     BVSetRandomColumn(U1,0);
657:     BVOrthonormalizeColumn(U1,0,PETSC_TRUE,NULL,NULL);
658:   }

660:   l = 0;
661:   while (svd->reason == SVD_CONVERGED_ITERATING) {
662:     svd->its++;

664:     /* inner loop */
665:     nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
666:     DSGetArrayReal(svd->ds,DS_MAT_T,&alpha);
667:     beta = alpha + ld;
668:     SVDTwoSideLanczosGSingle(svd,alpha,beta,lanczos->Z,V,U1,lanczos->ksp,svd->nconv+l,&nv,&breakdown);
669:     DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha);
670:     BVSetActiveColumns(V,svd->nconv,nv);
671:     BVSetActiveColumns(U1,svd->nconv,nv);

673:     /* solve projected problem */
674:     DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l);
675:     DSSVDSetDimensions(svd->ds,nv);
676:     if (l==0) {
677:       DSSetState(svd->ds,DS_STATE_INTERMEDIATE);
678:     } else {
679:       DSSetState(svd->ds,DS_STATE_RAW);
680:     }
681:     DSSolve(svd->ds,w,NULL);
682:     DSSort(svd->ds,w,NULL,NULL,NULL,NULL);
683:     DSUpdateExtraRow(svd->ds);
684:     DSSynchronize(svd->ds,w,NULL);
685:     for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);

687:     /* check convergence */
688:     SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,&k);
689:     (*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx);

691:     /* update l */
692:     if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
693:     else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
694:     if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
695:     if (l) { PetscInfo1(svd,"Preparing to restart keeping l=%D vectors\n",l); }

697:     if (svd->reason == SVD_CONVERGED_ITERATING) {
698:       if (breakdown || k==nv) {
699:         /* Start a new bidiagonalization */
700:         PetscInfo1(svd,"Breakdown in bidiagonalization (it=%D)\n",svd->its);
701:         if (k<svd->nsv) {
702:           BVSetRandomColumn(U1,k);
703:           BVOrthonormalizeColumn(U1,k,PETSC_FALSE,NULL,&breakdown);
704:           if (breakdown) {
705:             svd->reason = SVD_DIVERGED_BREAKDOWN;
706:             PetscInfo(svd,"Unable to generate more start vectors\n");
707:           }
708:         }
709:       } else {
710:         DSTruncate(svd->ds,k+l,PETSC_FALSE);
711:       }
712:     }

714:     /* compute converged singular vectors and restart vectors */
715:     DSGetMat(svd->ds,DS_MAT_U,&U);
716:     BVMultInPlace(V,U,svd->nconv,k+l);
717:     MatDestroy(&U);
718:     DSGetMat(svd->ds,DS_MAT_V,&VV);
719:     BVMultInPlace(U1,VV,svd->nconv,k+l);
720:     MatDestroy(&VV);

722:     /* copy the last vector to be the next initial vector */
723:     if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) {
724:       BVCopyColumn(U1,nv,k+l);
725:     }

727:     svd->nconv = k;
728:     SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv);
729:   }

731:   PetscFree(w);
732:   return(0);
733: }

735: /* Move generalized left singular vectors (0..nconv) from U1 and U2 to its final destination svd->U (single variant) */
736: PETSC_STATIC_INLINE PetscErrorCode SVDLeftSingularVectors_Single(SVD svd,BV U1,BV U2)
737: {
738:   PetscErrorCode    ierr;
739:   PetscInt          i,k,m,p;
740:   Vec               u,u1,u2;
741:   PetscScalar       *ua,*u2a;
742:   const PetscScalar *u1a;
743:   PetscReal         s;

746:   MatGetLocalSize(svd->A,&m,NULL);
747:   MatGetLocalSize(svd->B,&p,NULL);
748:   for (i=0;i<svd->nconv;i++) {
749:     BVGetColumn(U1,i,&u1);
750:     BVGetColumn(U2,i,&u2);
751:     BVGetColumn(svd->U,i,&u);
752:     VecGetArrayRead(u1,&u1a);
753:     VecGetArray(u,&ua);
754:     VecGetArray(u2,&u2a);
755:     /* Copy column from U1 to upper part of u */
756:     for (k=0;k<m;k++) ua[k] = u1a[k];
757:     /* Copy column from lower part of U to U2. Orthogonalize column in U2 and copy back to U */
758:     for (k=0;k<p;k++) u2a[k] = ua[m+k];
759:     VecRestoreArray(u2,&u2a);
760:     BVRestoreColumn(U2,i,&u2);
761:     BVOrthonormalizeColumn(U2,i,PETSC_FALSE,&s,NULL);
762:     BVGetColumn(U2,i,&u2);
763:     VecGetArray(u2,&u2a);
764:     for (k=0;k<p;k++) ua[m+k] = u2a[k];
765:     /* Update singular value */
766:     svd->sigma[i] /= s;
767:     VecRestoreArrayRead(u1,&u1a);
768:     VecRestoreArray(u,&ua);
769:     VecRestoreArray(u2,&u2a);
770:     BVRestoreColumn(U1,i,&u1);
771:     BVRestoreColumn(U2,i,&u2);
772:     BVRestoreColumn(svd->U,i,&u);
773:   }
774:   return(0);
775: }

777: 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)
778: {
779:   PetscErrorCode    ierr;
780:   PetscInt          i,j,m,p;
781:   const PetscScalar *carr;
782:   PetscScalar       *arr,*u2arr;
783:   Vec               u,v,ut=svd->workl[0],x=svd->workr[0],v1,u2;
784:   PetscBool         lindep=PETSC_FALSE,lindep1=PETSC_FALSE,lindep2=PETSC_FALSE;

787:   MatCreateVecsEmpty(svd->A,NULL,&v1);
788:   MatGetLocalSize(svd->A,&m,NULL);
789:   MatGetLocalSize(svd->B,&p,NULL);

791:   for (i=k; i<*n; i++) {
792:     /* Compute vector i of BV U1 */
793:     BVGetColumn(V,i,&v);
794:     VecGetArrayRead(v,&carr);
795:     VecPlaceArray(v1,carr);
796:     BVInsertVec(U1,i,v1);
797:     VecResetArray(v1);
798:     BVOrthonormalizeColumn(U1,i,PETSC_FALSE,alpha+i,&lindep1);

800:     /* Compute vector i of BV U2 */
801:     BVGetColumn(U2,i,&u2);
802:     VecGetArray(u2,&u2arr);
803:     if (i%2) {
804:       for (j=0; j<p; j++) u2arr[j] = -carr[m+j];
805:     } else {
806:       for (j=0; j<p; j++) u2arr[j] = carr[m+j];
807:     }
808:     VecRestoreArray(u2,&u2arr);
809:     BVRestoreColumn(U2,i,&u2);
810:     VecRestoreArrayRead(v,&carr);
811:     BVRestoreColumn(V,i,&v);
812:     BVOrthonormalizeColumn(U2,i,PETSC_FALSE,alphah+i,&lindep2);
813:     if (i%2) alphah[i] = -alphah[i];
814:     if (lindep1 || lindep2) {
815:       lindep = PETSC_TRUE;
816:       *n = i;
817:       break;
818:     }

820:     /* Compute vector i+1 of BV V */
821:     BVGetColumn(V,i+1,&v);
822:     /* Form ut=[u;0] */
823:     BVGetColumn(U1,i,&u);
824:     VecZeroEntries(ut);
825:     VecGetArrayRead(u,&carr);
826:     VecGetArray(ut,&arr);
827:     for (j=0; j<m; j++) arr[j] = carr[j];
828:     VecRestoreArrayRead(u,&carr);
829:     VecRestoreArray(ut,&arr);
830:     /* Solve least squares problem */
831:     KSPSolve(ksp,ut,x);
832:     MatMult(Z,x,v);
833:     BVRestoreColumn(U1,i,&u);
834:     BVRestoreColumn(V,i+1,&v);
835:     BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,beta+i,&lindep);
836:     betah[i] = -alpha[i]*beta[i]/alphah[i];
837:     if (lindep) {
838:       *n = i;
839:       break;
840:     }
841:   }
842:   if (breakdown) *breakdown = lindep;
843:   VecDestroy(&v1);
844:   return(0);
845: }

847: /* generate random initial vector in column k for joint upper-upper bidiagonalization */
848: PETSC_STATIC_INLINE PetscErrorCode SVDInitialVectorGUpper(SVD svd,BV V,BV U1,PetscInt k,PetscBool *breakdown)
849: {
850:   PetscErrorCode    ierr;
851:   SVD_TRLANCZOS     *lanczos = (SVD_TRLANCZOS*)svd->data;
852:   Vec               u,v,ut=svd->workl[0],x=svd->workr[0];
853:   PetscInt          m,j;
854:   PetscRandom       rand;
855:   PetscScalar       *arr;
856:   const PetscScalar *carr;

859:   BVCreateVec(U1,&u);
860:   VecGetLocalSize(u,&m);
861:   BVGetRandomContext(U1,&rand);
862:   VecSetRandom(u,rand);
863:   /* Form ut=[u;0] */
864:   VecZeroEntries(ut);
865:   VecGetArrayRead(u,&carr);
866:   VecGetArray(ut,&arr);
867:   for (j=0; j<m; j++) arr[j] = carr[j];
868:   VecRestoreArrayRead(u,&carr);
869:   VecRestoreArray(ut,&arr);
870:   VecDestroy(&u);
871:   /* Solve least squares problem and premultiply the result by Z */
872:   KSPSolve(lanczos->ksp,ut,x);
873:   BVGetColumn(V,k,&v);
874:   MatMult(lanczos->Z,x,v);
875:   BVRestoreColumn(V,k,&v);
876:   if (breakdown) { BVOrthonormalizeColumn(V,k,PETSC_FALSE,NULL,breakdown); }
877:   else { BVOrthonormalizeColumn(V,k,PETSC_TRUE,NULL,NULL); }
878:   return(0);
879: }

881: /* solve generalized problem with joint upper-upper bidiagonalization */
882: PetscErrorCode SVDSolve_TRLanczosGUpper(SVD svd,BV U1,BV U2,BV V)
883: {
885:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
886:   PetscReal      *alpha,*beta,*alphah,*betah;
887:   PetscScalar    *w;
888:   PetscInt       i,k,l,nv,ld;
889:   Mat            U,Vmat,X;
890:   PetscBool      breakdown=PETSC_FALSE;

893:   DSGetLeadingDimension(svd->ds,&ld);
894:   PetscMalloc1(ld,&w);

896:   /* normalize start vector */
897:   if (!svd->nini) {
898:     SVDInitialVectorGUpper(svd,V,U1,0,NULL);
899:   }

901:   l = 0;
902:   while (svd->reason == SVD_CONVERGED_ITERATING) {
903:     svd->its++;

905:     /* inner loop */
906:     nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
907:     DSGetArrayReal(svd->ds,DS_MAT_T,&alpha);
908:     DSGetArrayReal(svd->ds,DS_MAT_D,&alphah);
909:     beta = alpha + ld;
910:     betah = alpha + 2*ld;
911:     SVDTwoSideLanczosGUpper(svd,alpha,beta,alphah,betah,lanczos->Z,U1,U2,V,lanczos->ksp,svd->nconv+l,&nv,&breakdown);
912:     DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha);
913:     DSRestoreArrayReal(svd->ds,DS_MAT_D,&alphah);
914:     BVSetActiveColumns(V,svd->nconv,nv);
915:     BVSetActiveColumns(U1,svd->nconv,nv);
916:     BVSetActiveColumns(U2,svd->nconv,nv);

918:     /* solve projected problem */
919:     DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l);
920:     DSGSVDSetDimensions(svd->ds,nv,nv);
921:     if (l==0) {
922:       DSSetState(svd->ds,DS_STATE_INTERMEDIATE);
923:     } else {
924:       DSSetState(svd->ds,DS_STATE_RAW);
925:     }
926:     DSSolve(svd->ds,w,NULL);
927:     DSSort(svd->ds,w,NULL,NULL,NULL,NULL);
928:     DSUpdateExtraRow(svd->ds);
929:     DSSynchronize(svd->ds,w,NULL);
930:     for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);

932:     /* check convergence */
933:     SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,&k);
934:     (*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx);

936:     /* update l */
937:     if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
938:     else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
939:     if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
940:     if (l) { PetscInfo1(svd,"Preparing to restart keeping l=%D vectors\n",l); }

942:     if (svd->reason == SVD_CONVERGED_ITERATING) {
943:       if (breakdown || k==nv) {
944:         /* Start a new bidiagonalization */
945:         PetscInfo1(svd,"Breakdown in bidiagonalization (it=%D)\n",svd->its);
946:         if (k<svd->nsv) {
947:           SVDInitialVectorGUpper(svd,V,U1,k,&breakdown);
948:           if (breakdown) {
949:             svd->reason = SVD_DIVERGED_BREAKDOWN;
950:             PetscInfo(svd,"Unable to generate more start vectors\n");
951:           }
952:         }
953:       } else {
954:         DSTruncate(svd->ds,k+l,PETSC_FALSE);
955:       }
956:     }
957:     /* compute converged singular vectors and restart vectors */
958:     DSGetMat(svd->ds,DS_MAT_X,&X);
959:     BVMultInPlace(V,X,svd->nconv,k+l);
960:     MatDestroy(&X);
961:     DSGetMat(svd->ds,DS_MAT_U,&U);
962:     BVMultInPlace(U1,U,svd->nconv,k+l);
963:     MatDestroy(&U);
964:     DSGetMat(svd->ds,DS_MAT_V,&Vmat);
965:     BVMultInPlace(U2,Vmat,svd->nconv,k+l);
966:     MatDestroy(&Vmat);

968:     /* copy the last vector to be the next initial vector */
969:     if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) {
970:       BVCopyColumn(V,nv,k+l);
971:     }

973:     svd->nconv = k;
974:     SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv);
975:   }

977:   PetscFree(w);
978:   return(0);
979: }

981: /* Move generalized left singular vectors (0..nconv) from U1 and U2 to its final destination svd->U (upper and lower variants) */
982: PETSC_STATIC_INLINE PetscErrorCode SVDLeftSingularVectors(SVD svd,BV U1,BV U2)
983: {
984:   PetscErrorCode    ierr;
985:   PetscInt          i,k,m,p;
986:   Vec               u,u1,u2;
987:   PetscScalar       *ua;
988:   const PetscScalar *u1a,*u2a;

991:   MatGetLocalSize(svd->A,&m,NULL);
992:   MatGetLocalSize(svd->B,&p,NULL);
993:   for (i=0;i<svd->nconv;i++) {
994:     BVGetColumn(U1,i,&u1);
995:     BVGetColumn(U2,i,&u2);
996:     BVGetColumn(svd->U,i,&u);
997:     VecGetArrayRead(u1,&u1a);
998:     VecGetArrayRead(u2,&u2a);
999:     VecGetArray(u,&ua);
1000:     /* Copy column from u1 to upper part of u */
1001:     for (k=0;k<m;k++) ua[k] = u1a[k];
1002:     /* Copy column from u2 to lower part of u */
1003:     for (k=0;k<p;k++) ua[m+k] = u2a[k];
1004:     VecRestoreArrayRead(u1,&u1a);
1005:     VecRestoreArrayRead(u2,&u2a);
1006:     VecRestoreArray(u,&ua);
1007:     BVRestoreColumn(U1,i,&u1);
1008:     BVRestoreColumn(U2,i,&u2);
1009:     BVRestoreColumn(svd->U,i,&u);
1010:   }
1011:   return(0);
1012: }

1014: 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)
1015: {
1016:   PetscErrorCode    ierr;
1017:   PetscInt          i,j,m,p;
1018:   const PetscScalar *carr;
1019:   PetscScalar       *arr,*u2arr;
1020:   Vec               u,v,ut=svd->workl[0],x=svd->workr[0],v1,u2;
1021:   PetscBool         lindep=PETSC_FALSE;

1024:   MatCreateVecsEmpty(svd->A,NULL,&v1);
1025:   MatGetLocalSize(svd->A,&m,NULL);
1026:   MatGetLocalSize(svd->B,&p,NULL);

1028:   for (i=k; i<*n; i++) {
1029:     /* Compute vector i of BV U2 */
1030:     BVGetColumn(V,i,&v);
1031:     VecGetArrayRead(v,&carr);
1032:     BVGetColumn(U2,i,&u2);
1033:     VecGetArray(u2,&u2arr);
1034:     if (i%2) {
1035:       for (j=0; j<p; j++) u2arr[j] = -carr[m+j];
1036:     } else {
1037:       for (j=0; j<p; j++) u2arr[j] = carr[m+j];
1038:     }
1039:     VecRestoreArray(u2,&u2arr);
1040:     BVRestoreColumn(U2,i,&u2);
1041:     BVOrthonormalizeColumn(U2,i,PETSC_FALSE,alphah+i,&lindep);
1042:     if (i%2) alphah[i] = -alphah[i];
1043:     if (lindep) {
1044:       BVRestoreColumn(V,i,&v);
1045:       *n = i;
1046:       break;
1047:     }

1049:     /* Compute vector i+1 of BV U1 */
1050:     VecPlaceArray(v1,carr);
1051:     BVInsertVec(U1,i+1,v1);
1052:     VecResetArray(v1);
1053:     BVOrthonormalizeColumn(U1,i+1,PETSC_FALSE,beta+i,&lindep);
1054:     VecRestoreArrayRead(v,&carr);
1055:     BVRestoreColumn(V,i,&v);
1056:     if (lindep) {
1057:       *n = i+1;
1058:       break;
1059:     }

1061:     /* Compute vector i+1 of BV V */
1062:     BVGetColumn(V,i+1,&v);
1063:     /* Form ut=[u;0] where u is column i+1 of BV U1 */
1064:     BVGetColumn(U1,i+1,&u);
1065:     VecZeroEntries(ut);
1066:     VecGetArrayRead(u,&carr);
1067:     VecGetArray(ut,&arr);
1068:     for (j=0; j<m; j++) arr[j] = carr[j];
1069:     VecRestoreArrayRead(u,&carr);
1070:     VecRestoreArray(ut,&arr);
1071:     /* Solve least squares problem */
1072:     KSPSolve(ksp,ut,x);
1073:     MatMult(Z,x,v);
1074:     BVRestoreColumn(U1,i+1,&u);
1075:     BVRestoreColumn(V,i+1,&v);
1076:     BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,alpha+i+1,&lindep);
1077:     betah[i] = -alpha[i+1]*beta[i]/alphah[i];
1078:     if (lindep) {
1079:       *n = i+1;
1080:       break;
1081:     }
1082:   }
1083:   if (breakdown) *breakdown = lindep;
1084:   VecDestroy(&v1);
1085:   return(0);
1086: }

1088: /* generate random initial vector in column k for joint lower-upper bidiagonalization */
1089: PETSC_STATIC_INLINE PetscErrorCode SVDInitialVectorGLower(SVD svd,BV V,BV U1,PetscInt k,PetscBool *breakdown)
1090: {
1091:   PetscErrorCode    ierr;
1092:   SVD_TRLANCZOS     *lanczos = (SVD_TRLANCZOS*)svd->data;
1093:   const PetscScalar *carr;
1094:   PetscScalar       *arr;
1095:   PetscReal         *alpha;
1096:   PetscInt          j,m;
1097:   Vec               u,v,ut=svd->workl[0],x=svd->workr[0];

1100:   BVSetRandomColumn(U1,k);
1101:   if (breakdown) { BVOrthonormalizeColumn(U1,k,PETSC_FALSE,NULL,breakdown); }
1102:   else { BVOrthonormalizeColumn(U1,k,PETSC_TRUE,NULL,NULL); }

1104:   if (!breakdown || !*breakdown) {
1105:     MatGetLocalSize(svd->A,&m,NULL);
1106:     /* Compute k-th vector of BV V */
1107:     BVGetColumn(V,k,&v);
1108:     /* Form ut=[u;0] where u is the 1st column of U1 */
1109:     BVGetColumn(U1,k,&u);
1110:     VecZeroEntries(ut);
1111:     VecGetArrayRead(u,&carr);
1112:     VecGetArray(ut,&arr);
1113:     for (j=0; j<m; j++) arr[j] = carr[j];
1114:     VecRestoreArrayRead(u,&carr);
1115:     VecRestoreArray(ut,&arr);
1116:     /* Solve least squares problem */
1117:     KSPSolve(lanczos->ksp,ut,x);
1118:     MatMult(lanczos->Z,x,v);
1119:     BVRestoreColumn(U1,k,&u);
1120:     BVRestoreColumn(V,k,&v);
1121:     DSGetArrayReal(svd->ds,DS_MAT_T,&alpha);
1122:     if (breakdown) { BVOrthonormalizeColumn(V,k,PETSC_FALSE,alpha+k,breakdown); }
1123:     else { BVOrthonormalizeColumn(V,k,PETSC_TRUE,alpha+k,NULL); }
1124:     DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha);
1125:   }
1126:   return(0);
1127: }

1129: /* solve generalized problem with joint lower-upper bidiagonalization */
1130: PetscErrorCode SVDSolve_TRLanczosGLower(SVD svd,BV U1,BV U2,BV V)
1131: {
1133:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
1134:   PetscReal      *alpha,*beta,*alphah,*betah;
1135:   PetscScalar    *w;
1136:   PetscInt       i,k,l,nv,ld;
1137:   Mat            U,Vmat,X;
1138:   PetscBool      breakdown=PETSC_FALSE;

1141:   DSGetLeadingDimension(svd->ds,&ld);
1142:   PetscMalloc1(ld,&w);

1144:   /* normalize start vector */
1145:   if (!svd->nini) {
1146:     SVDInitialVectorGLower(svd,V,U1,0,NULL);
1147:   }

1149:   l = 0;
1150:   while (svd->reason == SVD_CONVERGED_ITERATING) {
1151:     svd->its++;

1153:     /* inner loop */
1154:     nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
1155:     DSGetArrayReal(svd->ds,DS_MAT_T,&alpha);
1156:     DSGetArrayReal(svd->ds,DS_MAT_D,&alphah);
1157:     beta = alpha + ld;
1158:     betah = alpha + 2*ld;
1159:     SVDTwoSideLanczosGLower(svd,alpha,beta,alphah,betah,lanczos->Z,U1,U2,V,lanczos->ksp,svd->nconv+l,&nv,&breakdown);
1160:     DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha);
1161:     DSRestoreArrayReal(svd->ds,DS_MAT_D,&alphah);
1162:     BVSetActiveColumns(V,svd->nconv,nv);
1163:     BVSetActiveColumns(U1,svd->nconv,nv+1);
1164:     BVSetActiveColumns(U2,svd->nconv,nv);

1166:     /* solve projected problem */
1167:     DSSetDimensions(svd->ds,nv+1,svd->nconv,svd->nconv+l);
1168:     DSGSVDSetDimensions(svd->ds,nv,nv);
1169:     if (l==0) {
1170:       DSSetState(svd->ds,DS_STATE_INTERMEDIATE);
1171:     } else {
1172:       DSSetState(svd->ds,DS_STATE_RAW);
1173:     }
1174:     DSSolve(svd->ds,w,NULL);
1175:     DSSort(svd->ds,w,NULL,NULL,NULL,NULL);
1176:     DSUpdateExtraRow(svd->ds);
1177:     DSSynchronize(svd->ds,w,NULL);
1178:     for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);

1180:     /* check convergence */
1181:     SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,&k);
1182:     (*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx);

1184:     /* update l */
1185:     if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
1186:     else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
1187:     if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
1188:     if (l) { PetscInfo1(svd,"Preparing to restart keeping l=%D vectors\n",l); }

1190:     if (svd->reason == SVD_CONVERGED_ITERATING) {
1191:       if (breakdown || k==nv) {
1192:         /* Start a new bidiagonalization */
1193:         PetscInfo1(svd,"Breakdown in bidiagonalization (it=%D)\n",svd->its);
1194:         if (k<svd->nsv) {
1195:           SVDInitialVectorGLower(svd,V,U1,k,&breakdown);
1196:           if (breakdown) {
1197:             svd->reason = SVD_DIVERGED_BREAKDOWN;
1198:             PetscInfo(svd,"Unable to generate more start vectors\n");
1199:           }
1200:         }
1201:       } else {
1202:         DSTruncate(svd->ds,k+l,PETSC_FALSE);
1203:       }
1204:     }

1206:     /* compute converged singular vectors and restart vectors */
1207:     DSGetMat(svd->ds,DS_MAT_X,&X);
1208:     BVMultInPlace(V,X,svd->nconv,k+l);
1209:     MatDestroy(&X);
1210:     DSGetMat(svd->ds,DS_MAT_U,&U);
1211:     BVMultInPlace(U1,U,svd->nconv,k+l+1);
1212:     MatDestroy(&U);
1213:     DSGetMat(svd->ds,DS_MAT_V,&Vmat);
1214:     BVMultInPlace(U2,Vmat,svd->nconv,k+l);
1215:     MatDestroy(&Vmat);

1217:     /* copy the last vector to be the next initial vector */
1218:     if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) {
1219:       BVCopyColumn(V,nv,k+l);
1220:     }

1222:     svd->nconv = k;
1223:     SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv);
1224:   }

1226:   PetscFree(w);
1227:   return(0);
1228: }

1230: PetscErrorCode SVDSolve_TRLanczos_GSVD(SVD svd)
1231: {
1233:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
1234:   PetscInt       k,m,p;
1235:   BV             U1,U2;
1236:   BVType         type;
1237:   Mat            U,V;

1240:   PetscCitationsRegister(citationg,&citedg);

1242:   MatGetLocalSize(svd->A,&m,NULL);
1243:   MatGetLocalSize(svd->B,&p,NULL);

1245:   /* Create BV for U1 */
1246:   BVCreate(PetscObjectComm((PetscObject)svd),&U1);
1247:   PetscLogObjectParent((PetscObject)svd,(PetscObject)U1);
1248:   BVGetType(svd->U,&type);
1249:   BVSetType(U1,type);
1250:   BVGetSizes(svd->U,NULL,NULL,&k);
1251:   BVSetSizes(U1,m,PETSC_DECIDE,k);

1253:   /* Create BV for U2 */
1254:   BVCreate(PetscObjectComm((PetscObject)svd),&U2);
1255:   PetscLogObjectParent((PetscObject)svd,(PetscObject)U2);
1256:   BVSetType(U2,type);
1257:   BVSetSizes(U2,p,PETSC_DECIDE,k);

1259:   switch (lanczos->bidiag) {
1260:     case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1261:       SVDSolve_TRLanczosGSingle(svd,U1,svd->U);
1262:       break;
1263:     case SVD_TRLANCZOS_GBIDIAG_UPPER:
1264:       SVDSolve_TRLanczosGUpper(svd,U1,U2,svd->U);
1265:       break;
1266:     case SVD_TRLANCZOS_GBIDIAG_LOWER:
1267:       SVDSolve_TRLanczosGLower(svd,U1,U2,svd->U);
1268:       break;
1269:   }

1271:   /* Compute converged right singular vectors */
1272:   BVSetActiveColumns(svd->U,0,svd->nconv);
1273:   BVSetActiveColumns(svd->V,0,svd->nconv);
1274:   BVGetMat(svd->U,&U);
1275:   BVGetMat(svd->V,&V);
1276:   KSPMatSolve(lanczos->ksp,U,V);
1277:   BVRestoreMat(svd->U,&U);
1278:   BVRestoreMat(svd->V,&V);

1280:   /* Finish computing left singular vectors and move them to its place */
1281:   switch (lanczos->bidiag) {
1282:     case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1283:       SVDLeftSingularVectors_Single(svd,U1,U2);
1284:       break;
1285:     case SVD_TRLANCZOS_GBIDIAG_UPPER:
1286:     case SVD_TRLANCZOS_GBIDIAG_LOWER:
1287:       SVDLeftSingularVectors(svd,U1,U2);
1288:       break;
1289:   }

1291:   BVDestroy(&U2);
1292:   BVDestroy(&U1);
1293:   DSTruncate(svd->ds,svd->nconv,PETSC_TRUE);
1294:   return(0);
1295: }

1297: PetscErrorCode SVDSetFromOptions_TRLanczos(PetscOptionItems *PetscOptionsObject,SVD svd)
1298: {
1299:   PetscErrorCode      ierr;
1300:   SVD_TRLANCZOS       *lanczos = (SVD_TRLANCZOS*)svd->data;
1301:   PetscBool           flg,val,lock;
1302:   PetscReal           keep;
1303:   SVDTRLanczosGBidiag bidiag;

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

1308:     PetscOptionsBool("-svd_trlanczos_oneside","Use one-side reorthogonalization","SVDTRLanczosSetOneSide",lanczos->oneside,&val,&flg);
1309:     if (flg) { SVDTRLanczosSetOneSide(svd,val); }

1311:     PetscOptionsReal("-svd_trlanczos_restart","Proportion of vectors kept after restart","SVDTRLanczosSetRestart",0.5,&keep,&flg);
1312:     if (flg) { SVDTRLanczosSetRestart(svd,keep); }

1314:     PetscOptionsBool("-svd_trlanczos_locking","Choose between locking and non-locking variants","SVDTRLanczosSetLocking",PETSC_TRUE,&lock,&flg);
1315:     if (flg) { SVDTRLanczosSetLocking(svd,lock); }

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

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

1323:   PetscOptionsTail();

1325:   if (svd->OPb) {
1326:     if (!lanczos->ksp) { SVDTRLanczosGetKSP(svd,&lanczos->ksp); }
1327:     KSPSetFromOptions(lanczos->ksp);
1328:   }
1329:   return(0);
1330: }

1332: static PetscErrorCode SVDTRLanczosSetOneSide_TRLanczos(SVD svd,PetscBool oneside)
1333: {
1334:   SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;

1337:   if (lanczos->oneside != oneside) {
1338:     lanczos->oneside = oneside;
1339:     svd->state = SVD_STATE_INITIAL;
1340:   }
1341:   return(0);
1342: }

1344: /*@
1345:    SVDTRLanczosSetOneSide - Indicate if the variant of the Lanczos method
1346:    to be used is one-sided or two-sided.

1348:    Logically Collective on svd

1350:    Input Parameters:
1351: +  svd     - singular value solver
1352: -  oneside - boolean flag indicating if the method is one-sided or not

1354:    Options Database Key:
1355: .  -svd_trlanczos_oneside <boolean> - Indicates the boolean flag

1357:    Note:
1358:    By default, a two-sided variant is selected, which is sometimes slightly
1359:    more robust. However, the one-sided variant is faster because it avoids
1360:    the orthogonalization associated to left singular vectors.

1362:    Level: advanced

1364: .seealso: SVDLanczosSetOneSide()
1365: @*/
1366: PetscErrorCode SVDTRLanczosSetOneSide(SVD svd,PetscBool oneside)
1367: {

1373:   PetscTryMethod(svd,"SVDTRLanczosSetOneSide_C",(SVD,PetscBool),(svd,oneside));
1374:   return(0);
1375: }

1377: static PetscErrorCode SVDTRLanczosGetOneSide_TRLanczos(SVD svd,PetscBool *oneside)
1378: {
1379:   SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;

1382:   *oneside = lanczos->oneside;
1383:   return(0);
1384: }

1386: /*@
1387:    SVDTRLanczosGetOneSide - Gets if the variant of the Lanczos method
1388:    to be used is one-sided or two-sided.

1390:    Not Collective

1392:    Input Parameters:
1393: .  svd     - singular value solver

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

1398:    Level: advanced

1400: .seealso: SVDTRLanczosSetOneSide()
1401: @*/
1402: PetscErrorCode SVDTRLanczosGetOneSide(SVD svd,PetscBool *oneside)
1403: {

1409:   PetscUseMethod(svd,"SVDTRLanczosGetOneSide_C",(SVD,PetscBool*),(svd,oneside));
1410:   return(0);
1411: }

1413: static PetscErrorCode SVDTRLanczosSetGBidiag_TRLanczos(SVD svd,SVDTRLanczosGBidiag bidiag)
1414: {
1415:   SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;

1418:   switch (bidiag) {
1419:     case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1420:     case SVD_TRLANCZOS_GBIDIAG_UPPER:
1421:     case SVD_TRLANCZOS_GBIDIAG_LOWER:
1422:       if (lanczos->bidiag != bidiag) {
1423:         lanczos->bidiag = bidiag;
1424:         svd->state = SVD_STATE_INITIAL;
1425:       }
1426:       break;
1427:     default:
1428:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid bidiagonalization choice");
1429:   }
1430:   return(0);
1431: }

1433: /*@
1434:    SVDTRLanczosSetGBidiag - Sets the bidiagonalization choice to use in
1435:    the GSVD TRLanczos solver.

1437:    Logically Collective on svd

1439:    Input Parameters:
1440: +  svd - the singular value solver
1441: -  bidiag - the bidiagonalization choice

1443:    Options Database Key:
1444: .  -svd_trlanczos_gbidiag - Sets the bidiagonalization choice (either 's' or 'juu'
1445:    or 'jlu')

1447:    Level: advanced

1449: .seealso: SVDTRLanczosGetGBidiag(), SVDTRLanczosGBidiag
1450: @*/
1451: PetscErrorCode SVDTRLanczosSetGBidiag(SVD svd,SVDTRLanczosGBidiag bidiag)
1452: {

1458:   PetscTryMethod(svd,"SVDTRLanczosSetGBidiag_C",(SVD,SVDTRLanczosGBidiag),(svd,bidiag));
1459:   return(0);
1460: }

1462: static PetscErrorCode SVDTRLanczosGetGBidiag_TRLanczos(SVD svd,SVDTRLanczosGBidiag *bidiag)
1463: {
1464:   SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;

1467:   *bidiag = lanczos->bidiag;
1468:   return(0);
1469: }

1471: /*@
1472:    SVDTRLanczosGetGBidiag - Gets the bidiagonalization choice used in the GSVD
1473:    TRLanczos solver.

1475:    Not Collective

1477:    Input Parameter:
1478: .  svd - the singular value solver

1480:    Output Parameter:
1481: .  bidiag - the bidiagonalization choice

1483:    Level: advanced

1485: .seealso: SVDTRLanczosSetGBidiag(), SVDTRLanczosGBidiag
1486: @*/
1487: PetscErrorCode SVDTRLanczosGetGBidiag(SVD svd,SVDTRLanczosGBidiag *bidiag)
1488: {

1494:   PetscUseMethod(svd,"SVDTRLanczosGetGBidiag_C",(SVD,SVDTRLanczosGBidiag*),(svd,bidiag));
1495:   return(0);
1496: }

1498: static PetscErrorCode SVDTRLanczosSetKSP_TRLanczos(SVD svd,KSP ksp)
1499: {
1501:   SVD_TRLANCZOS  *ctx = (SVD_TRLANCZOS*)svd->data;

1504:   PetscObjectReference((PetscObject)ksp);
1505:   KSPDestroy(&ctx->ksp);
1506:   ctx->ksp = ksp;
1507:   PetscLogObjectParent((PetscObject)svd,(PetscObject)ctx->ksp);
1508:   svd->state = SVD_STATE_INITIAL;
1509:   return(0);
1510: }

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

1515:    Collective on svd

1517:    Input Parameters:
1518: +  svd - SVD solver
1519: -  ksp - the linear solver object

1521:    Note:
1522:    Only used for the GSVD problem.

1524:    Level: advanced

1526: .seealso: SVDTRLanczosGetKSP()
1527: @*/
1528: PetscErrorCode SVDTRLanczosSetKSP(SVD svd,KSP ksp)
1529: {

1536:   PetscTryMethod(svd,"SVDTRLanczosSetKSP_C",(SVD,KSP),(svd,ksp));
1537:   return(0);
1538: }

1540: static PetscErrorCode SVDTRLanczosGetKSP_TRLanczos(SVD svd,KSP *ksp)
1541: {
1543:   SVD_TRLANCZOS  *ctx = (SVD_TRLANCZOS*)svd->data;
1544:   PC             pc;

1547:   if (!ctx->ksp) {
1548:     /* Create linear solver */
1549:     KSPCreate(PetscObjectComm((PetscObject)svd),&ctx->ksp);
1550:     PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)svd,1);
1551:     KSPSetOptionsPrefix(ctx->ksp,((PetscObject)svd)->prefix);
1552:     KSPAppendOptionsPrefix(ctx->ksp,"svd_trlanczos_");
1553:     PetscLogObjectParent((PetscObject)svd,(PetscObject)ctx->ksp);
1554:     PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)svd)->options);
1555:     KSPSetType(ctx->ksp,KSPLSQR);
1556:     KSPGetPC(ctx->ksp,&pc);
1557:     PCSetType(pc,PCNONE);
1558:     KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE);
1559:     KSPSetTolerances(ctx->ksp,SlepcDefaultTol(svd->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1560:   }
1561:   *ksp = ctx->ksp;
1562:   return(0);
1563: }

1565: /*@
1566:    SVDTRLanczosGetKSP - Retrieve the linear solver object (KSP) associated with
1567:    the SVD solver.

1569:    Not Collective

1571:    Input Parameter:
1572: .  svd - SVD solver

1574:    Output Parameter:
1575: .  ksp - the linear solver object

1577:    Level: advanced

1579: .seealso: SVDTRLanczosSetKSP()
1580: @*/
1581: PetscErrorCode SVDTRLanczosGetKSP(SVD svd,KSP *ksp)
1582: {

1588:   PetscUseMethod(svd,"SVDTRLanczosGetKSP_C",(SVD,KSP*),(svd,ksp));
1589:   return(0);
1590: }

1592: static PetscErrorCode SVDTRLanczosSetRestart_TRLanczos(SVD svd,PetscReal keep)
1593: {
1594:   SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;

1597:   if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
1598:   else {
1599:     if (keep<0.1 || keep>0.9) SETERRQ1(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument %g must be in the range [0.1,0.9]",keep);
1600:     ctx->keep = keep;
1601:   }
1602:   return(0);
1603: }

1605: /*@
1606:    SVDTRLanczosSetRestart - Sets the restart parameter for the thick-restart
1607:    Lanczos method, in particular the proportion of basis vectors that must be
1608:    kept after restart.

1610:    Logically Collective on svd

1612:    Input Parameters:
1613: +  svd  - the singular value solver
1614: -  keep - the number of vectors to be kept at restart

1616:    Options Database Key:
1617: .  -svd_trlanczos_restart - Sets the restart parameter

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

1622:    Level: advanced

1624: .seealso: SVDTRLanczosGetRestart()
1625: @*/
1626: PetscErrorCode SVDTRLanczosSetRestart(SVD svd,PetscReal keep)
1627: {

1633:   PetscTryMethod(svd,"SVDTRLanczosSetRestart_C",(SVD,PetscReal),(svd,keep));
1634:   return(0);
1635: }

1637: static PetscErrorCode SVDTRLanczosGetRestart_TRLanczos(SVD svd,PetscReal *keep)
1638: {
1639:   SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;

1642:   *keep = ctx->keep;
1643:   return(0);
1644: }

1646: /*@
1647:    SVDTRLanczosGetRestart - Gets the restart parameter used in the thick-restart
1648:    Lanczos method.

1650:    Not Collective

1652:    Input Parameter:
1653: .  svd - the singular value solver

1655:    Output Parameter:
1656: .  keep - the restart parameter

1658:    Level: advanced

1660: .seealso: SVDTRLanczosSetRestart()
1661: @*/
1662: PetscErrorCode SVDTRLanczosGetRestart(SVD svd,PetscReal *keep)
1663: {

1669:   PetscUseMethod(svd,"SVDTRLanczosGetRestart_C",(SVD,PetscReal*),(svd,keep));
1670:   return(0);
1671: }

1673: static PetscErrorCode SVDTRLanczosSetLocking_TRLanczos(SVD svd,PetscBool lock)
1674: {
1675:   SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;

1678:   ctx->lock = lock;
1679:   return(0);
1680: }

1682: /*@
1683:    SVDTRLanczosSetLocking - Choose between locking and non-locking variants of
1684:    the thick-restart Lanczos method.

1686:    Logically Collective on svd

1688:    Input Parameters:
1689: +  svd  - the singular value solver
1690: -  lock - true if the locking variant must be selected

1692:    Options Database Key:
1693: .  -svd_trlanczos_locking - Sets the locking flag

1695:    Notes:
1696:    The default is to lock converged singular triplets when the method restarts.
1697:    This behaviour can be changed so that all directions are kept in the
1698:    working subspace even if already converged to working accuracy (the
1699:    non-locking variant).

1701:    Level: advanced

1703: .seealso: SVDTRLanczosGetLocking()
1704: @*/
1705: PetscErrorCode SVDTRLanczosSetLocking(SVD svd,PetscBool lock)
1706: {

1712:   PetscTryMethod(svd,"SVDTRLanczosSetLocking_C",(SVD,PetscBool),(svd,lock));
1713:   return(0);
1714: }

1716: static PetscErrorCode SVDTRLanczosGetLocking_TRLanczos(SVD svd,PetscBool *lock)
1717: {
1718:   SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;

1721:   *lock = ctx->lock;
1722:   return(0);
1723: }

1725: /*@
1726:    SVDTRLanczosGetLocking - Gets the locking flag used in the thick-restart
1727:    Lanczos method.

1729:    Not Collective

1731:    Input Parameter:
1732: .  svd - the singular value solver

1734:    Output Parameter:
1735: .  lock - the locking flag

1737:    Level: advanced

1739: .seealso: SVDTRLanczosSetLocking()
1740: @*/
1741: PetscErrorCode SVDTRLanczosGetLocking(SVD svd,PetscBool *lock)
1742: {

1748:   PetscUseMethod(svd,"SVDTRLanczosGetLocking_C",(SVD,PetscBool*),(svd,lock));
1749:   return(0);
1750: }

1752: static PetscErrorCode SVDTRLanczosSetExplicitMatrix_TRLanczos(SVD svd,PetscBool explicitmatrix)
1753: {
1754:   SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;

1757:   if (lanczos->explicitmatrix != explicitmatrix) {
1758:     lanczos->explicitmatrix = explicitmatrix;
1759:     svd->state = SVD_STATE_INITIAL;
1760:   }
1761:   return(0);
1762: }

1764: /*@
1765:    SVDTRLanczosSetExplicitMatrix - Indicate if the matrix Z=[A;B] must
1766:    be built explicitly.

1768:    Logically Collective on svd

1770:    Input Parameters:
1771: +  svd      - singular value solver
1772: -  explicit - Boolean flag indicating if Z=[A;B] is built explicitly

1774:    Options Database Key:
1775: .  -svd_trlanczos_explicitmatrix <boolean> - Indicates the boolean flag

1777:    Notes:
1778:    This option is relevant for the GSVD case only.
1779:    Z is the coefficient matrix of the KSP solver used internally.

1781:    Level: advanced

1783: .seealso: SVDTRLanczosGetExplicitMatrix()
1784: @*/
1785: PetscErrorCode SVDTRLanczosSetExplicitMatrix(SVD svd,PetscBool explicitmatrix)
1786: {

1792:   PetscTryMethod(svd,"SVDTRLanczosSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmatrix));
1793:   return(0);
1794: }

1796: static PetscErrorCode SVDTRLanczosGetExplicitMatrix_TRLanczos(SVD svd,PetscBool *explicitmatrix)
1797: {
1798:   SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;

1801:   *explicitmatrix = lanczos->explicitmatrix;
1802:   return(0);
1803: }

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

1808:    Not Collective

1810:    Input Parameter:
1811: .  svd  - singular value solver

1813:    Output Parameter:
1814: .  explicit - the mode flag

1816:    Level: advanced

1818: .seealso: SVDTRLanczosSetExplicitMatrix()
1819: @*/
1820: PetscErrorCode SVDTRLanczosGetExplicitMatrix(SVD svd,PetscBool *explicitmatrix)
1821: {

1827:   PetscUseMethod(svd,"SVDTRLanczosGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmatrix));
1828:   return(0);
1829: }

1831: PetscErrorCode SVDReset_TRLanczos(SVD svd)
1832: {
1834:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;

1837:   if (svd->isgeneralized) {
1838:     KSPReset(lanczos->ksp);
1839:     MatDestroy(&lanczos->Z);
1840:   }
1841:   return(0);
1842: }

1844: PetscErrorCode SVDDestroy_TRLanczos(SVD svd)
1845: {
1847:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;

1850:   if (svd->isgeneralized) {
1851:     KSPDestroy(&lanczos->ksp);
1852:   }
1853:   PetscFree(svd->data);
1854:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",NULL);
1855:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",NULL);
1856:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetGBidiag_C",NULL);
1857:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetGBidiag_C",NULL);
1858:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetKSP_C",NULL);
1859:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetKSP_C",NULL);
1860:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetRestart_C",NULL);
1861:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetRestart_C",NULL);
1862:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetLocking_C",NULL);
1863:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetLocking_C",NULL);
1864:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetExplicitMatrix_C",NULL);
1865:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetExplicitMatrix_C",NULL);
1866:   return(0);
1867: }

1869: PetscErrorCode SVDView_TRLanczos(SVD svd,PetscViewer viewer)
1870: {
1872:   SVD_TRLANCZOS  *lanczos = (SVD_TRLANCZOS*)svd->data;
1873:   PetscBool      isascii;

1876:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1877:   if (isascii) {
1878:     PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*lanczos->keep));
1879:     PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",lanczos->lock?"":"non-");
1880:     if (svd->isgeneralized) {
1881:       const char *bidiag="";

1883:       switch (lanczos->bidiag) {
1884:         case SVD_TRLANCZOS_GBIDIAG_SINGLE: bidiag = "single"; break;
1885:         case SVD_TRLANCZOS_GBIDIAG_UPPER:  bidiag = "joint upper-upper"; break;
1886:         case SVD_TRLANCZOS_GBIDIAG_LOWER:  bidiag = "joint lower-upper"; break;
1887:       }
1888:       PetscViewerASCIIPrintf(viewer,"  bidiagonalization choice: %s\n",bidiag);
1889:       PetscViewerASCIIPrintf(viewer,"  %s matrix\n",lanczos->explicitmatrix?"explicit":"implicit");
1890:       if (!lanczos->ksp) { SVDTRLanczosGetKSP(svd,&lanczos->ksp); }
1891:       PetscViewerASCIIPushTab(viewer);
1892:       KSPView(lanczos->ksp,viewer);
1893:       PetscViewerASCIIPopTab(viewer);
1894:     } else {
1895:       PetscViewerASCIIPrintf(viewer,"  %s-sided reorthogonalization\n",lanczos->oneside? "one": "two");
1896:     }
1897:   }
1898:   return(0);
1899: }

1901: SLEPC_EXTERN PetscErrorCode SVDCreate_TRLanczos(SVD svd)
1902: {
1904:   SVD_TRLANCZOS  *ctx;

1907:   PetscNewLog(svd,&ctx);
1908:   svd->data = (void*)ctx;

1910:   ctx->lock   = PETSC_TRUE;
1911:   ctx->bidiag = SVD_TRLANCZOS_GBIDIAG_LOWER;

1913:   svd->ops->setup          = SVDSetUp_TRLanczos;
1914:   svd->ops->solve          = SVDSolve_TRLanczos;
1915:   svd->ops->solveg         = SVDSolve_TRLanczos_GSVD;
1916:   svd->ops->destroy        = SVDDestroy_TRLanczos;
1917:   svd->ops->reset          = SVDReset_TRLanczos;
1918:   svd->ops->setfromoptions = SVDSetFromOptions_TRLanczos;
1919:   svd->ops->view           = SVDView_TRLanczos;
1920:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",SVDTRLanczosSetOneSide_TRLanczos);
1921:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",SVDTRLanczosGetOneSide_TRLanczos);
1922:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetGBidiag_C",SVDTRLanczosSetGBidiag_TRLanczos);
1923:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetGBidiag_C",SVDTRLanczosGetGBidiag_TRLanczos);
1924:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetKSP_C",SVDTRLanczosSetKSP_TRLanczos);
1925:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetKSP_C",SVDTRLanczosGetKSP_TRLanczos);
1926:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetRestart_C",SVDTRLanczosSetRestart_TRLanczos);
1927:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetRestart_C",SVDTRLanczosGetRestart_TRLanczos);
1928:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetLocking_C",SVDTRLanczosSetLocking_TRLanczos);
1929:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetLocking_C",SVDTRLanczosGetLocking_TRLanczos);
1930:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetExplicitMatrix_C",SVDTRLanczosSetExplicitMatrix_TRLanczos);
1931:   PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetExplicitMatrix_C",SVDTRLanczosGetExplicitMatrix_TRLanczos);
1932:   return(0);
1933: }