Actual source code: trlanczos.c
slepc-3.22.1 2024-10-28
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>
32: #include <slepc/private/bvimpl.h>
34: static PetscBool cited = PETSC_FALSE,citedg = PETSC_FALSE;
35: static const char citation[] =
36: "@Article{slepc-svd,\n"
37: " author = \"V. Hern{\\'a}ndez and J. E. Rom{\\'a}n and A. Tom{\\'a}s\",\n"
38: " title = \"A robust and efficient parallel {SVD} solver based on restarted {Lanczos} bidiagonalization\",\n"
39: " journal = \"Electron. Trans. Numer. Anal.\",\n"
40: " volume = \"31\",\n"
41: " pages = \"68--85\",\n"
42: " year = \"2008\"\n"
43: "}\n";
44: static const char citationg[] =
45: "@Article{slepc-gsvd,\n"
46: " author = \"F. Alvarruiz and C. Campos and J. E. Roman\",\n"
47: " title = \"Thick-restarted {Lanczos} bidiagonalization methods for the {GSVD}\",\n"
48: " journal = \"J. Comput. Appl. Math.\",\n"
49: " volume = \"440\",\n"
50: " pages = \"115506\",\n"
51: " year = \"2024\"\n"
52: "}\n";
54: typedef struct {
55: /* user parameters */
56: PetscBool oneside; /* one-sided variant */
57: PetscReal keep; /* restart parameter */
58: PetscBool lock; /* locking/non-locking variant */
59: KSP ksp; /* solver for least-squares problem in GSVD */
60: SVDTRLanczosGBidiag bidiag; /* bidiagonalization variant for GSVD */
61: PetscReal scalef; /* scale factor for matrix B */
62: PetscReal scaleth; /* scale threshold for automatic scaling */
63: PetscBool explicitmatrix;
64: /* auxiliary variables */
65: Mat Z; /* aux matrix for GSVD, Z=[A;B] */
66: } SVD_TRLANCZOS;
68: /* Context for shell matrix [A; B] */
69: typedef struct {
70: Mat A,B,AT,BT;
71: Vec y1,y2,y;
72: PetscInt m;
73: PetscReal scalef;
74: } MatZData;
76: static PetscErrorCode MatZCreateContext(SVD svd,MatZData **zdata)
77: {
78: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
80: PetscFunctionBegin;
81: PetscCall(PetscNew(zdata));
82: (*zdata)->A = svd->A;
83: (*zdata)->B = svd->B;
84: (*zdata)->AT = svd->AT;
85: (*zdata)->BT = svd->BT;
86: (*zdata)->scalef = lanczos->scalef;
87: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&(*zdata)->y1));
88: PetscCall(MatCreateVecsEmpty(svd->B,NULL,&(*zdata)->y2));
89: PetscCall(VecGetLocalSize((*zdata)->y1,&(*zdata)->m));
90: PetscCall(BVCreateVec(svd->U,&(*zdata)->y));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: /* Update scale factor for B in Z=[A;B]
95: If matrices are swapped, the scale factor is inverted.*/
96: static PetscErrorCode MatZUpdateScale(SVD svd)
97: {
98: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
99: MatZData *zdata;
100: Mat mats[2],normal;
101: MatType Atype;
102: PetscBool sametype;
103: PetscReal scalef = svd->swapped? 1.0/lanczos->scalef : lanczos->scalef;
105: PetscFunctionBegin;
106: if (lanczos->explicitmatrix) {
107: /* Destroy the matrix Z and create it again */
108: PetscCall(MatDestroy(&lanczos->Z));
109: mats[0] = svd->A;
110: if (scalef == 1.0) {
111: mats[1] = svd->B;
112: } else {
113: PetscCall(MatDuplicate(svd->B,MAT_COPY_VALUES,&mats[1]));
114: PetscCall(MatScale(mats[1],scalef));
115: }
116: PetscCall(MatCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,1,NULL,mats,&lanczos->Z));
117: PetscCall(MatGetType(svd->A,&Atype));
118: PetscCall(PetscObjectTypeCompare((PetscObject)svd->B,Atype,&sametype));
119: if (!sametype) Atype = MATAIJ;
120: PetscCall(MatConvert(lanczos->Z,Atype,MAT_INPLACE_MATRIX,&lanczos->Z));
121: if (scalef != 1.0) PetscCall(MatDestroy(&mats[1]));
122: } else {
123: PetscCall(MatShellGetContext(lanczos->Z,&zdata));
124: zdata->scalef = scalef;
125: }
127: /* create normal equations matrix, to build the preconditioner in LSQR */
128: PetscCall(MatCreateNormalHermitian(lanczos->Z,&normal));
130: if (!lanczos->ksp) PetscCall(SVDTRLanczosGetKSP(svd,&lanczos->ksp));
131: PetscCall(SVD_KSPSetOperators(lanczos->ksp,lanczos->Z,normal));
132: PetscCall(KSPSetUp(lanczos->ksp));
133: PetscCall(MatDestroy(&normal));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: static PetscErrorCode MatDestroy_Z(Mat Z)
138: {
139: MatZData *zdata;
141: PetscFunctionBegin;
142: PetscCall(MatShellGetContext(Z,&zdata));
143: PetscCall(VecDestroy(&zdata->y1));
144: PetscCall(VecDestroy(&zdata->y2));
145: PetscCall(VecDestroy(&zdata->y));
146: PetscCall(PetscFree(zdata));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: static PetscErrorCode MatMult_Z(Mat Z,Vec x,Vec y)
151: {
152: MatZData *zdata;
153: PetscScalar *y_elems;
155: PetscFunctionBegin;
156: PetscCall(MatShellGetContext(Z,&zdata));
157: PetscCall(VecGetArray(y,&y_elems));
158: PetscCall(VecPlaceArray(zdata->y1,y_elems));
159: PetscCall(VecPlaceArray(zdata->y2,y_elems+zdata->m));
161: PetscCall(MatMult(zdata->A,x,zdata->y1));
162: PetscCall(MatMult(zdata->B,x,zdata->y2));
163: PetscCall(VecScale(zdata->y2,zdata->scalef));
165: PetscCall(VecResetArray(zdata->y1));
166: PetscCall(VecResetArray(zdata->y2));
167: PetscCall(VecRestoreArray(y,&y_elems));
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: static PetscErrorCode MatMultTranspose_Z(Mat Z,Vec y,Vec x)
172: {
173: MatZData *zdata;
174: const PetscScalar *y_elems;
176: PetscFunctionBegin;
177: PetscCall(MatShellGetContext(Z,&zdata));
178: PetscCall(VecGetArrayRead(y,&y_elems));
179: PetscCall(VecPlaceArray(zdata->y1,y_elems));
180: PetscCall(VecPlaceArray(zdata->y2,y_elems+zdata->m));
182: PetscCall(MatMult(zdata->BT,zdata->y2,x));
183: PetscCall(VecScale(x,zdata->scalef));
184: PetscCall(MatMultAdd(zdata->AT,zdata->y1,x,x));
186: PetscCall(VecResetArray(zdata->y1));
187: PetscCall(VecResetArray(zdata->y2));
188: PetscCall(VecRestoreArrayRead(y,&y_elems));
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: static PetscErrorCode MatCreateVecs_Z(Mat Z,Vec *right,Vec *left)
193: {
194: MatZData *zdata;
196: PetscFunctionBegin;
197: PetscCall(MatShellGetContext(Z,&zdata));
198: if (right) PetscCall(MatCreateVecs(zdata->A,right,NULL));
199: if (left) PetscCall(VecDuplicate(zdata->y,left));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: static PetscErrorCode SVDSetUp_TRLanczos(SVD svd)
204: {
205: PetscInt M,N,P,m,n,p;
206: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
207: MatZData *zdata;
208: Mat aux;
210: PetscFunctionBegin;
211: PetscCall(MatGetSize(svd->A,&M,&N));
212: PetscCall(SVDSetDimensions_Default(svd));
213: PetscCheck(svd->ncv<=svd->nsv+svd->mpd,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nsv+mpd");
214: PetscCheck(lanczos->lock || svd->mpd>=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
215: if (svd->max_it==PETSC_DETERMINE) svd->max_it = PetscMax(N/svd->ncv,100);
216: if (!lanczos->keep) lanczos->keep = 0.5;
217: svd->leftbasis = PETSC_TRUE;
218: PetscCall(SVDAllocateSolution(svd,1));
219: if (svd->isgeneralized) {
220: PetscCall(MatGetSize(svd->B,&P,NULL));
221: if (lanczos->bidiag == SVD_TRLANCZOS_GBIDIAG_LOWER && ((svd->which==SVD_LARGEST && P<=N) || (svd->which==SVD_SMALLEST && M>N && P<=N))) {
222: SlepcSwap(svd->A,svd->B,aux);
223: SlepcSwap(svd->AT,svd->BT,aux);
224: svd->swapped = PETSC_TRUE;
225: } else svd->swapped = PETSC_FALSE;
227: PetscCall(SVDSetWorkVecs(svd,1,1));
229: if (svd->conv==SVD_CONV_ABS) { /* Residual norms are multiplied by matrix norms */
230: if (!svd->nrma) PetscCall(MatNorm(svd->A,NORM_INFINITY,&svd->nrma));
231: if (!svd->nrmb) PetscCall(MatNorm(svd->B,NORM_INFINITY,&svd->nrmb));
232: }
234: /* Create the matrix Z=[A;B] */
235: PetscCall(MatGetLocalSize(svd->A,&m,&n));
236: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
237: if (!lanczos->explicitmatrix) {
238: PetscCall(MatDestroy(&lanczos->Z));
239: PetscCall(MatZCreateContext(svd,&zdata));
240: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)svd),m+p,n,PETSC_DECIDE,PETSC_DECIDE,zdata,&lanczos->Z));
241: PetscCall(MatShellSetOperation(lanczos->Z,MATOP_MULT,(void(*)(void))MatMult_Z));
242: #if defined(PETSC_USE_COMPLEX)
243: PetscCall(MatShellSetOperation(lanczos->Z,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultTranspose_Z));
244: #else
245: PetscCall(MatShellSetOperation(lanczos->Z,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Z));
246: #endif
247: PetscCall(MatShellSetOperation(lanczos->Z,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_Z));
248: PetscCall(MatShellSetOperation(lanczos->Z,MATOP_DESTROY,(void(*)(void))MatDestroy_Z));
249: }
250: /* Explicit matrix is created here, when updating the scale */
251: PetscCall(MatZUpdateScale(svd));
253: } else if (svd->ishyperbolic) {
254: PetscCall(BV_SetMatrixDiagonal(svd->swapped?svd->V:svd->U,svd->omega,svd->OP));
255: PetscCall(SVDSetWorkVecs(svd,1,0));
256: }
257: PetscCall(DSSetCompact(svd->ds,PETSC_TRUE));
258: PetscCall(DSSetExtraRow(svd->ds,PETSC_TRUE));
259: PetscCall(DSAllocate(svd->ds,svd->ncv+1));
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: static PetscErrorCode SVDOneSideTRLanczosMGS(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,BV U,PetscInt nconv,PetscInt l,PetscInt n,PetscScalar* work)
264: {
265: PetscReal a,b;
266: PetscInt i,k=nconv+l;
267: Vec ui,ui1,vi;
269: PetscFunctionBegin;
270: PetscCall(BVGetColumn(V,k,&vi));
271: PetscCall(BVGetColumn(U,k,&ui));
272: PetscCall(MatMult(svd->A,vi,ui));
273: PetscCall(BVRestoreColumn(V,k,&vi));
274: PetscCall(BVRestoreColumn(U,k,&ui));
275: if (l>0) {
276: PetscCall(BVSetActiveColumns(U,nconv,n));
277: for (i=0;i<l;i++) work[i]=beta[i+nconv];
278: PetscCall(BVMultColumn(U,-1.0,1.0,k,work));
279: }
280: PetscCall(BVNormColumn(U,k,NORM_2,&a));
281: PetscCall(BVScaleColumn(U,k,1.0/a));
282: alpha[k] = a;
284: for (i=k+1;i<n;i++) {
285: PetscCall(BVGetColumn(V,i,&vi));
286: PetscCall(BVGetColumn(U,i-1,&ui1));
287: PetscCall(MatMult(svd->AT,ui1,vi));
288: PetscCall(BVRestoreColumn(V,i,&vi));
289: PetscCall(BVRestoreColumn(U,i-1,&ui1));
290: PetscCall(BVOrthonormalizeColumn(V,i,PETSC_FALSE,&b,NULL));
291: beta[i-1] = b;
293: PetscCall(BVGetColumn(V,i,&vi));
294: PetscCall(BVGetColumn(U,i,&ui));
295: PetscCall(MatMult(svd->A,vi,ui));
296: PetscCall(BVRestoreColumn(V,i,&vi));
297: PetscCall(BVGetColumn(U,i-1,&ui1));
298: PetscCall(VecAXPY(ui,-b,ui1));
299: PetscCall(BVRestoreColumn(U,i-1,&ui1));
300: PetscCall(BVRestoreColumn(U,i,&ui));
301: PetscCall(BVNormColumn(U,i,NORM_2,&a));
302: PetscCall(BVScaleColumn(U,i,1.0/a));
303: alpha[i] = a;
304: }
306: PetscCall(BVGetColumn(V,n,&vi));
307: PetscCall(BVGetColumn(U,n-1,&ui1));
308: PetscCall(MatMult(svd->AT,ui1,vi));
309: PetscCall(BVRestoreColumn(V,n,&vi));
310: PetscCall(BVRestoreColumn(U,n-1,&ui1));
311: PetscCall(BVOrthogonalizeColumn(V,n,NULL,&b,NULL));
312: beta[n-1] = b;
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: /*
317: Custom CGS orthogonalization, preprocess after first orthogonalization
318: */
319: static PetscErrorCode SVDOrthogonalizeCGS(BV V,PetscInt i,PetscScalar* h,PetscReal a,BVOrthogRefineType refine,PetscReal eta,PetscReal *norm)
320: {
321: PetscReal sum,onorm;
322: PetscScalar dot;
323: PetscInt j;
325: PetscFunctionBegin;
326: switch (refine) {
327: case BV_ORTHOG_REFINE_NEVER:
328: PetscCall(BVNormColumn(V,i,NORM_2,norm));
329: break;
330: case BV_ORTHOG_REFINE_ALWAYS:
331: PetscCall(BVSetActiveColumns(V,0,i));
332: PetscCall(BVDotColumn(V,i,h));
333: PetscCall(BVMultColumn(V,-1.0,1.0,i,h));
334: PetscCall(BVNormColumn(V,i,NORM_2,norm));
335: break;
336: case BV_ORTHOG_REFINE_IFNEEDED:
337: dot = h[i];
338: onorm = PetscSqrtReal(PetscRealPart(dot)) / a;
339: sum = 0.0;
340: for (j=0;j<i;j++) {
341: sum += PetscRealPart(h[j] * PetscConj(h[j]));
342: }
343: *norm = PetscRealPart(dot)/(a*a) - sum;
344: if (*norm>0.0) *norm = PetscSqrtReal(*norm);
345: else PetscCall(BVNormColumn(V,i,NORM_2,norm));
346: if (*norm < eta*onorm) {
347: PetscCall(BVSetActiveColumns(V,0,i));
348: PetscCall(BVDotColumn(V,i,h));
349: PetscCall(BVMultColumn(V,-1.0,1.0,i,h));
350: PetscCall(BVNormColumn(V,i,NORM_2,norm));
351: }
352: break;
353: }
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: static PetscErrorCode SVDOneSideTRLanczosCGS(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,BV U,PetscInt nconv,PetscInt l,PetscInt n,PetscScalar* work)
358: {
359: PetscReal a,b,eta;
360: PetscInt i,j,k=nconv+l;
361: Vec ui,ui1,vi;
362: BVOrthogRefineType refine;
364: PetscFunctionBegin;
365: PetscCall(BVGetColumn(V,k,&vi));
366: PetscCall(BVGetColumn(U,k,&ui));
367: PetscCall(MatMult(svd->A,vi,ui));
368: PetscCall(BVRestoreColumn(V,k,&vi));
369: PetscCall(BVRestoreColumn(U,k,&ui));
370: if (l>0) {
371: PetscCall(BVSetActiveColumns(U,nconv,n));
372: for (i=0;i<l;i++) work[i]=beta[i+nconv];
373: PetscCall(BVMultColumn(U,-1.0,1.0,k,work));
374: }
375: PetscCall(BVGetOrthogonalization(V,NULL,&refine,&eta,NULL));
377: for (i=k+1;i<n;i++) {
378: PetscCall(BVGetColumn(V,i,&vi));
379: PetscCall(BVGetColumn(U,i-1,&ui1));
380: PetscCall(MatMult(svd->AT,ui1,vi));
381: PetscCall(BVRestoreColumn(V,i,&vi));
382: PetscCall(BVRestoreColumn(U,i-1,&ui1));
383: PetscCall(BVNormColumnBegin(U,i-1,NORM_2,&a));
384: if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
385: PetscCall(BVSetActiveColumns(V,0,i+1));
386: PetscCall(BVGetColumn(V,i,&vi));
387: PetscCall(BVDotVecBegin(V,vi,work));
388: } else {
389: PetscCall(BVSetActiveColumns(V,0,i));
390: PetscCall(BVDotColumnBegin(V,i,work));
391: }
392: PetscCall(BVNormColumnEnd(U,i-1,NORM_2,&a));
393: if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
394: PetscCall(BVDotVecEnd(V,vi,work));
395: PetscCall(BVRestoreColumn(V,i,&vi));
396: PetscCall(BVSetActiveColumns(V,0,i));
397: } else PetscCall(BVDotColumnEnd(V,i,work));
399: PetscCall(BVScaleColumn(U,i-1,1.0/a));
400: for (j=0;j<i;j++) work[j] = work[j] / a;
401: PetscCall(BVMultColumn(V,-1.0,1.0/a,i,work));
402: PetscCall(SVDOrthogonalizeCGS(V,i,work,a,refine,eta,&b));
403: PetscCall(BVScaleColumn(V,i,1.0/b));
404: PetscCheck(PetscAbsReal(b)>10*PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)svd),PETSC_ERR_PLIB,"Recurrence generated a zero vector; use a two-sided variant");
406: PetscCall(BVGetColumn(V,i,&vi));
407: PetscCall(BVGetColumn(U,i,&ui));
408: PetscCall(BVGetColumn(U,i-1,&ui1));
409: PetscCall(MatMult(svd->A,vi,ui));
410: PetscCall(VecAXPY(ui,-b,ui1));
411: PetscCall(BVRestoreColumn(V,i,&vi));
412: PetscCall(BVRestoreColumn(U,i,&ui));
413: PetscCall(BVRestoreColumn(U,i-1,&ui1));
415: alpha[i-1] = a;
416: beta[i-1] = b;
417: }
419: PetscCall(BVGetColumn(V,n,&vi));
420: PetscCall(BVGetColumn(U,n-1,&ui1));
421: PetscCall(MatMult(svd->AT,ui1,vi));
422: PetscCall(BVRestoreColumn(V,n,&vi));
423: PetscCall(BVRestoreColumn(U,n-1,&ui1));
425: PetscCall(BVNormColumnBegin(svd->U,n-1,NORM_2,&a));
426: if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
427: PetscCall(BVSetActiveColumns(V,0,n+1));
428: PetscCall(BVGetColumn(V,n,&vi));
429: PetscCall(BVDotVecBegin(V,vi,work));
430: } else {
431: PetscCall(BVSetActiveColumns(V,0,n));
432: PetscCall(BVDotColumnBegin(V,n,work));
433: }
434: PetscCall(BVNormColumnEnd(svd->U,n-1,NORM_2,&a));
435: if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
436: PetscCall(BVDotVecEnd(V,vi,work));
437: PetscCall(BVRestoreColumn(V,n,&vi));
438: } else PetscCall(BVDotColumnEnd(V,n,work));
440: PetscCall(BVScaleColumn(U,n-1,1.0/a));
441: for (j=0;j<n;j++) work[j] = work[j] / a;
442: PetscCall(BVMultColumn(V,-1.0,1.0/a,n,work));
443: PetscCall(SVDOrthogonalizeCGS(V,n,work,a,refine,eta,&b));
444: PetscCall(BVSetActiveColumns(V,nconv,n));
445: alpha[n-1] = a;
446: beta[n-1] = b;
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: static PetscErrorCode SVDSolve_TRLanczos(SVD svd)
451: {
452: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
453: PetscReal *alpha,*beta;
454: PetscScalar *swork=NULL,*w;
455: PetscInt i,k,l,nv,ld;
456: Mat U,V;
457: PetscBool breakdown=PETSC_FALSE;
458: BVOrthogType orthog;
460: PetscFunctionBegin;
461: PetscCall(PetscCitationsRegister(citation,&cited));
462: /* allocate working space */
463: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
464: PetscCall(BVGetOrthogonalization(svd->V,&orthog,NULL,NULL,NULL));
465: PetscCall(PetscMalloc1(ld,&w));
466: if (lanczos->oneside) PetscCall(PetscMalloc1(svd->ncv+1,&swork));
468: /* normalize start vector */
469: if (!svd->nini) {
470: PetscCall(BVSetRandomColumn(svd->V,0));
471: PetscCall(BVOrthonormalizeColumn(svd->V,0,PETSC_TRUE,NULL,NULL));
472: }
474: l = 0;
475: while (svd->reason == SVD_CONVERGED_ITERATING) {
476: svd->its++;
478: /* inner loop */
479: nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
480: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
481: beta = alpha + ld;
482: if (lanczos->oneside) {
483: if (orthog == BV_ORTHOG_MGS) PetscCall(SVDOneSideTRLanczosMGS(svd,alpha,beta,svd->V,svd->U,svd->nconv,l,nv,swork));
484: else PetscCall(SVDOneSideTRLanczosCGS(svd,alpha,beta,svd->V,svd->U,svd->nconv,l,nv,swork));
485: } else PetscCall(SVDTwoSideLanczos(svd,alpha,beta,svd->V,svd->U,svd->nconv+l,&nv,&breakdown));
486: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
487: PetscCall(BVScaleColumn(svd->V,nv,1.0/beta[nv-1]));
488: PetscCall(BVSetActiveColumns(svd->V,svd->nconv,nv));
489: PetscCall(BVSetActiveColumns(svd->U,svd->nconv,nv));
491: /* solve projected problem */
492: PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l));
493: PetscCall(DSSVDSetDimensions(svd->ds,nv));
494: PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
495: PetscCall(DSSolve(svd->ds,w,NULL));
496: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
497: PetscCall(DSUpdateExtraRow(svd->ds));
498: PetscCall(DSSynchronize(svd->ds,w,NULL));
499: for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
501: /* check convergence */
502: PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,1.0,&k));
503: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
505: /* update l */
506: if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
507: else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
508: if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
509: if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
511: if (svd->reason == SVD_CONVERGED_ITERATING) {
512: if (PetscUnlikely(breakdown || k==nv)) {
513: /* Start a new bidiagonalization */
514: PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
515: if (k<svd->nsv) {
516: PetscCall(BVSetRandomColumn(svd->V,k));
517: PetscCall(BVOrthonormalizeColumn(svd->V,k,PETSC_FALSE,NULL,&breakdown));
518: if (breakdown) {
519: svd->reason = SVD_DIVERGED_BREAKDOWN;
520: PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
521: }
522: }
523: } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
524: }
526: /* compute converged singular vectors and restart vectors */
527: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&V));
528: PetscCall(BVMultInPlace(svd->V,V,svd->nconv,k+l));
529: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&V));
530: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
531: PetscCall(BVMultInPlace(svd->U,U,svd->nconv,k+l));
532: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
534: /* copy the last vector to be the next initial vector */
535: if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(svd->V,nv,k+l));
537: svd->nconv = k;
538: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv));
539: }
541: /* orthonormalize U columns in one side method */
542: if (lanczos->oneside) {
543: for (i=0;i<svd->nconv;i++) PetscCall(BVOrthonormalizeColumn(svd->U,i,PETSC_FALSE,NULL,NULL));
544: }
546: /* free working space */
547: PetscCall(PetscFree(w));
548: if (swork) PetscCall(PetscFree(swork));
549: PetscCall(DSTruncate(svd->ds,svd->nconv,PETSC_TRUE));
550: PetscFunctionReturn(PETSC_SUCCESS);
551: }
553: static PetscErrorCode SVDLanczosHSVD(SVD svd,PetscReal *alpha,PetscReal *beta,PetscReal *omega,Mat A,Mat AT,BV V,BV U,PetscInt k,PetscInt *n,PetscBool *breakdown)
554: {
555: PetscInt i;
556: Vec u,v,ou=svd->workl[0];
557: PetscBool lindep=PETSC_FALSE;
558: PetscReal norm;
560: PetscFunctionBegin;
561: for (i=k;i<*n;i++) {
562: PetscCall(BVGetColumn(V,i,&v));
563: PetscCall(BVGetColumn(U,i,&u));
564: PetscCall(MatMult(A,v,u));
565: PetscCall(BVRestoreColumn(V,i,&v));
566: PetscCall(BVRestoreColumn(U,i,&u));
567: PetscCall(BVOrthonormalizeColumn(U,i,PETSC_FALSE,alpha+i,&lindep));
568: omega[i] = PetscSign(alpha[i]);
569: if (PetscUnlikely(lindep)) {
570: *n = i;
571: break;
572: }
574: PetscCall(BVGetColumn(V,i+1,&v));
575: PetscCall(BVGetColumn(U,i,&u));
576: PetscCall(VecPointwiseMult(ou,svd->omega,u));
577: PetscCall(MatMult(AT,ou,v));
578: PetscCall(BVRestoreColumn(V,i+1,&v));
579: PetscCall(BVRestoreColumn(U,i,&u));
580: PetscCall(BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,&norm,&lindep));
581: beta[i] = omega[i]*norm;
582: if (PetscUnlikely(lindep)) {
583: *n = i+1;
584: break;
585: }
586: }
588: if (breakdown) *breakdown = lindep;
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: static PetscErrorCode SVDSolve_TRLanczos_HSVD(SVD svd)
593: {
594: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
595: PetscReal *alpha,*beta,*omega;
596: PetscScalar *w;
597: PetscInt i,k,l,nv,ld,nini;
598: Mat UU,VV,D,A,AT;
599: BV U,V;
600: PetscBool breakdown=PETSC_FALSE;
601: BVOrthogType orthog;
602: Vec vomega;
604: PetscFunctionBegin;
605: /* undo the effect of swapping in this function */
606: if (svd->swapped) {
607: A = svd->AT;
608: AT = svd->A;
609: U = svd->V;
610: V = svd->U;
611: nini = svd->ninil;
612: } else {
613: A = svd->A;
614: AT = svd->AT;
615: U = svd->U;
616: V = svd->V;
617: nini = svd->nini;
618: }
619: /* allocate working space */
620: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
621: PetscCall(BVGetOrthogonalization(V,&orthog,NULL,NULL,NULL));
622: PetscCall(PetscMalloc1(ld,&w));
623: PetscCheck(!lanczos->oneside,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Oneside orthogonalization not supported for HSVD");
625: /* normalize start vector */
626: if (!nini) {
627: PetscCall(BVSetRandomColumn(V,0));
628: PetscCall(BVOrthonormalizeColumn(V,0,PETSC_TRUE,NULL,NULL));
629: }
631: l = 0;
632: while (svd->reason == SVD_CONVERGED_ITERATING) {
633: svd->its++;
635: /* inner loop */
636: nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
637: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
638: beta = alpha + ld;
639: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&omega));
640: PetscCall(SVDLanczosHSVD(svd,alpha,beta,omega,A,AT,V,U,svd->nconv+l,&nv,&breakdown));
641: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
642: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&omega));
643: PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
644: PetscCall(BVSetActiveColumns(U,svd->nconv,nv));
646: /* solve projected problem */
647: PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l));
648: PetscCall(DSHSVDSetDimensions(svd->ds,nv));
649: PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
650: PetscCall(DSSolve(svd->ds,w,NULL));
651: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
652: PetscCall(DSUpdateExtraRow(svd->ds));
653: PetscCall(DSSynchronize(svd->ds,w,NULL));
654: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&omega));
655: for (i=svd->nconv;i<nv;i++) {
656: svd->sigma[i] = PetscRealPart(w[i]);
657: svd->sign[i] = omega[i];
658: }
659: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&omega));
661: /* check convergence */
662: PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,1.0,&k));
663: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
665: /* update l */
666: if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
667: else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
668: if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
669: if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
671: if (svd->reason == SVD_CONVERGED_ITERATING) {
672: if (PetscUnlikely(breakdown || k==nv)) {
673: /* Start a new bidiagonalization */
674: PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
675: if (k<svd->nsv) {
676: PetscCall(BVSetRandomColumn(V,k));
677: PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,NULL,&breakdown));
678: if (breakdown) {
679: svd->reason = SVD_DIVERGED_BREAKDOWN;
680: PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
681: }
682: }
683: } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
684: }
686: /* compute converged singular vectors and restart vectors */
687: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&VV));
688: PetscCall(BVMultInPlace(V,VV,svd->nconv,k+l));
689: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&VV));
690: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&UU));
691: PetscCall(BVMultInPlace(U,UU,svd->nconv,k+l));
692: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&UU));
694: /* copy the last vector of V to be the next initial vector
695: and change signature matrix of U */
696: if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) {
697: PetscCall(BVCopyColumn(V,nv,k+l));
698: PetscCall(BVSetActiveColumns(U,0,k+l));
699: PetscCall(DSGetMatAndColumn(svd->ds,DS_MAT_D,0,&D,&vomega));
700: PetscCall(BVSetSignature(U,vomega));
701: PetscCall(DSRestoreMatAndColumn(svd->ds,DS_MAT_D,0,&D,&vomega));
702: }
704: svd->nconv = k;
705: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv));
706: }
708: /* free working space */
709: PetscCall(PetscFree(w));
710: PetscCall(DSTruncate(svd->ds,svd->nconv,PETSC_TRUE));
711: PetscFunctionReturn(PETSC_SUCCESS);
712: }
714: /* Given n computed generalized singular values in sigmain, backtransform them
715: in sigmaout by undoing scaling and reciprocating if swapped=true. Also updates vectors V
716: if given. If sigmaout=NULL then the result overwrites sigmain. */
717: static PetscErrorCode SVDLanczosBackTransform(SVD svd,PetscInt n,PetscReal *sigmain,PetscReal *sigmaout,BV V)
718: {
719: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
720: PetscInt i;
721: PetscReal c,s,r,f,scalef;
723: PetscFunctionBegin;
724: scalef = svd->swapped? 1.0/lanczos->scalef: lanczos->scalef;
725: for (i=0;i<n;i++) {
726: if (V && scalef != 1.0) {
727: s = 1.0/PetscSqrtReal(1.0+sigmain[i]*sigmain[i]);
728: c = sigmain[i]*s;
729: r = s/scalef;
730: f = 1.0/PetscSqrtReal(c*c+r*r);
731: PetscCall(BVScaleColumn(V,i,f));
732: }
733: if (sigmaout) {
734: if (svd->swapped) sigmaout[i] = 1.0/(sigmain[i]*scalef);
735: else sigmaout[i] = sigmain[i]*scalef;
736: } else {
737: sigmain[i] *= scalef;
738: if (svd->swapped) sigmain[i] = 1.0/sigmain[i];
739: }
740: }
741: PetscFunctionReturn(PETSC_SUCCESS);
742: }
744: static PetscErrorCode SVDLanczosGSingle(SVD svd,PetscReal *alpha,PetscReal *beta,Mat Z,BV V,BV U,KSP ksp,PetscInt k,PetscInt *n,PetscBool *breakdown)
745: {
746: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
747: PetscInt i,j,m;
748: const PetscScalar *carr;
749: PetscScalar *arr;
750: Vec u,v,ut=svd->workl[0],x=svd->workr[0],v1,u1,u2;
751: PetscBool lindep=PETSC_FALSE;
753: PetscFunctionBegin;
754: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&v1));
755: PetscCall(BVGetColumn(V,k,&v));
756: PetscCall(BVGetColumn(U,k,&u));
758: /* Form ut=[u;0] */
759: PetscCall(VecZeroEntries(ut));
760: PetscCall(VecGetLocalSize(u,&m));
761: PetscCall(VecGetArrayRead(u,&carr));
762: PetscCall(VecGetArray(ut,&arr));
763: for (j=0; j<m; j++) arr[j] = carr[j];
764: PetscCall(VecRestoreArrayRead(u,&carr));
765: PetscCall(VecRestoreArray(ut,&arr));
767: /* Solve least squares problem */
768: PetscCall(KSPSolve(ksp,ut,x));
770: PetscCall(MatMult(Z,x,v));
772: PetscCall(BVRestoreColumn(U,k,&u));
773: PetscCall(BVRestoreColumn(V,k,&v));
774: PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,alpha+k,&lindep));
775: if (PetscUnlikely(lindep)) {
776: *n = k;
777: if (breakdown) *breakdown = lindep;
778: PetscFunctionReturn(PETSC_SUCCESS);
779: }
781: for (i=k+1; i<*n; i++) {
783: /* Compute vector i of BV U */
784: PetscCall(BVGetColumn(V,i-1,&v));
785: PetscCall(VecGetArray(v,&arr));
786: PetscCall(VecPlaceArray(v1,arr));
787: PetscCall(VecRestoreArray(v,&arr));
788: PetscCall(BVRestoreColumn(V,i-1,&v));
789: PetscCall(BVInsertVec(U,i,v1));
790: PetscCall(VecResetArray(v1));
791: PetscCall(BVOrthonormalizeColumn(U,i,PETSC_FALSE,beta+i-1,&lindep));
792: if (PetscUnlikely(lindep)) {
793: *n = i;
794: break;
795: }
797: /* Compute vector i of BV V */
799: PetscCall(BVGetColumn(V,i,&v));
800: PetscCall(BVGetColumn(U,i,&u));
802: /* Form ut=[u;0] */
803: PetscCall(VecGetArrayRead(u,&carr));
804: PetscCall(VecGetArray(ut,&arr));
805: for (j=0; j<m; j++) arr[j] = carr[j];
806: PetscCall(VecRestoreArrayRead(u,&carr));
807: PetscCall(VecRestoreArray(ut,&arr));
809: /* Solve least squares problem */
810: PetscCall(KSPSolve(ksp,ut,x));
812: PetscCall(MatMult(Z,x,v));
814: PetscCall(BVRestoreColumn(U,i,&u));
815: PetscCall(BVRestoreColumn(V,i,&v));
816: if (!lanczos->oneside || i==k+1) PetscCall(BVOrthonormalizeColumn(V,i,PETSC_FALSE,alpha+i,&lindep));
817: else { /* cheap computation of V[i], if restart (i==k+1) do a full reorthogonalization */
818: PetscCall(BVGetColumn(V,i,&u2));
819: PetscCall(BVGetColumn(V,i-1,&u1));
820: PetscCall(VecAXPY(u2,-beta[i-1],u1));
821: PetscCall(BVRestoreColumn(V,i-1,&u1));
822: PetscCall(VecNorm(u2,NORM_2,&alpha[i]));
823: if (alpha[i]==0.0) lindep = PETSC_TRUE;
824: else PetscCall(VecScale(u2,1.0/alpha[i]));
825: PetscCall(BVRestoreColumn(V,i,&u2));
826: }
827: if (PetscUnlikely(lindep)) {
828: *n = i;
829: break;
830: }
831: }
833: /* Compute vector n of BV U */
834: if (!lindep) {
835: PetscCall(BVGetColumn(V,*n-1,&v));
836: PetscCall(VecGetArray(v,&arr));
837: PetscCall(VecPlaceArray(v1,arr));
838: PetscCall(VecRestoreArray(v,&arr));
839: PetscCall(BVRestoreColumn(V,*n-1,&v));
840: PetscCall(BVInsertVec(U,*n,v1));
841: PetscCall(VecResetArray(v1));
842: PetscCall(BVOrthonormalizeColumn(U,*n,PETSC_FALSE,beta+*n-1,&lindep));
843: }
844: if (breakdown) *breakdown = lindep;
845: PetscCall(VecDestroy(&v1));
846: PetscFunctionReturn(PETSC_SUCCESS);
847: }
849: /* solve generalized problem with single bidiagonalization of Q_A */
850: static PetscErrorCode SVDSolve_TRLanczosGSingle(SVD svd,BV U1,BV V)
851: {
852: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
853: PetscReal *alpha,*beta,normr,scaleth,sigma0,*sigma;
854: PetscScalar *w;
855: PetscInt i,k,l,nv,ld;
856: Mat U,VV;
857: PetscBool breakdown=PETSC_FALSE;
859: PetscFunctionBegin;
860: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
861: PetscCall(PetscMalloc2(ld,&w,ld,&sigma));
862: normr = (svd->conv==SVD_CONV_ABS)? PetscMax(svd->nrma,svd->nrmb*lanczos->scalef): 1.0;
863: /* Convert scale threshold th=c/s to the corresponding c */
864: scaleth = (lanczos->scaleth!=0)? lanczos->scaleth/PetscSqrtReal(lanczos->scaleth*lanczos->scaleth+1): 0.0;
866: /* normalize start vector */
867: if (!svd->ninil) {
868: PetscCall(BVSetRandomColumn(U1,0));
869: PetscCall(BVOrthonormalizeColumn(U1,0,PETSC_TRUE,NULL,NULL));
870: }
872: l = 0;
873: while (svd->reason == SVD_CONVERGED_ITERATING) {
874: svd->its++;
876: /* inner loop */
877: nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
878: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
879: beta = alpha + ld;
880: PetscCall(SVDLanczosGSingle(svd,alpha,beta,lanczos->Z,V,U1,lanczos->ksp,svd->nconv+l,&nv,&breakdown));
881: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
882: PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
883: PetscCall(BVSetActiveColumns(U1,svd->nconv,nv));
885: /* solve projected problem */
886: PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l));
887: PetscCall(DSSVDSetDimensions(svd->ds,nv));
888: PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
889: PetscCall(DSSolve(svd->ds,w,NULL));
890: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
891: PetscCall(DSUpdateExtraRow(svd->ds));
892: PetscCall(DSSynchronize(svd->ds,w,NULL));
893: for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
895: /* check convergence */
896: PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,normr,&k));
897: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
899: sigma0 = svd->which==SVD_LARGEST? svd->sigma[0] : 1.0/svd->sigma[0];
900: if (scaleth!=0 && k==0 && sigma0>scaleth) {
902: /* Scale and start from scratch */
903: lanczos->scalef *= svd->sigma[0]/PetscSqrtReal(1-svd->sigma[0]*svd->sigma[0]);
904: PetscCall(PetscInfo(svd,"Scaling by factor %g and starting from scratch\n",(double)lanczos->scalef));
905: PetscCall(MatZUpdateScale(svd));
906: if (svd->conv==SVD_CONV_ABS) normr = PetscMax(svd->nrma,svd->nrmb*lanczos->scalef);
907: l = 0;
909: } else {
911: /* update l */
912: if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
913: else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
914: if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
915: if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
917: if (svd->reason == SVD_CONVERGED_ITERATING) {
918: if (PetscUnlikely(breakdown || k==nv)) {
919: /* Start a new bidiagonalization */
920: PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
921: if (k<svd->nsv) {
922: PetscCall(BVSetRandomColumn(U1,k));
923: PetscCall(BVOrthonormalizeColumn(U1,k,PETSC_FALSE,NULL,&breakdown));
924: if (breakdown) {
925: svd->reason = SVD_DIVERGED_BREAKDOWN;
926: PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
927: }
928: }
929: } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
930: }
932: /* compute converged singular vectors and restart vectors */
933: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
934: PetscCall(BVMultInPlace(V,U,svd->nconv,k+l));
935: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
936: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&VV));
937: PetscCall(BVMultInPlace(U1,VV,svd->nconv,k+l));
938: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&VV));
940: /* copy the last vector to be the next initial vector */
941: if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(U1,nv,k+l));
942: }
944: svd->nconv = k;
945: PetscCall(SVDLanczosBackTransform(svd,nv,svd->sigma,sigma,NULL));
946: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,sigma,svd->errest,nv));
947: }
949: PetscCall(PetscFree2(w,sigma));
950: PetscFunctionReturn(PETSC_SUCCESS);
951: }
953: /* Move generalized left singular vectors (0..nconv) from U1 and U2 to its final destination svd->U (single variant) */
954: static inline PetscErrorCode SVDLeftSingularVectors_Single(SVD svd,BV U1,BV U2)
955: {
956: PetscInt i,k,m,p;
957: Vec u,u1,u2;
958: PetscScalar *ua,*u2a;
959: const PetscScalar *u1a;
960: PetscReal s;
962: PetscFunctionBegin;
963: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
964: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
965: for (i=0;i<svd->nconv;i++) {
966: PetscCall(BVGetColumn(U1,i,&u1));
967: PetscCall(BVGetColumn(U2,i,&u2));
968: PetscCall(BVGetColumn(svd->U,i,&u));
969: PetscCall(VecGetArrayRead(u1,&u1a));
970: PetscCall(VecGetArray(u,&ua));
971: PetscCall(VecGetArray(u2,&u2a));
972: /* Copy column from U1 to upper part of u */
973: for (k=0;k<m;k++) ua[k] = u1a[k];
974: /* Copy column from lower part of U to U2. Orthogonalize column in U2 and copy back to U */
975: for (k=0;k<p;k++) u2a[k] = ua[m+k];
976: PetscCall(VecRestoreArray(u2,&u2a));
977: PetscCall(BVRestoreColumn(U2,i,&u2));
978: PetscCall(BVOrthonormalizeColumn(U2,i,PETSC_FALSE,&s,NULL));
979: PetscCall(BVGetColumn(U2,i,&u2));
980: PetscCall(VecGetArray(u2,&u2a));
981: for (k=0;k<p;k++) ua[m+k] = u2a[k];
982: /* Update singular value */
983: svd->sigma[i] /= s;
984: PetscCall(VecRestoreArrayRead(u1,&u1a));
985: PetscCall(VecRestoreArray(u,&ua));
986: PetscCall(VecRestoreArray(u2,&u2a));
987: PetscCall(BVRestoreColumn(U1,i,&u1));
988: PetscCall(BVRestoreColumn(U2,i,&u2));
989: PetscCall(BVRestoreColumn(svd->U,i,&u));
990: }
991: PetscFunctionReturn(PETSC_SUCCESS);
992: }
994: static PetscErrorCode SVDLanczosGUpper(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)
995: {
996: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
997: PetscInt i,j,m,p;
998: const PetscScalar *carr;
999: PetscScalar *arr,*u2arr;
1000: Vec u,v,ut=svd->workl[0],x=svd->workr[0],v1,u1,u2;
1001: PetscBool lindep=PETSC_FALSE,lindep1=PETSC_FALSE,lindep2=PETSC_FALSE;
1003: PetscFunctionBegin;
1004: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&v1));
1005: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1006: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1008: for (i=k; i<*n; i++) {
1009: /* Compute vector i of BV U1 */
1010: PetscCall(BVGetColumn(V,i,&v));
1011: PetscCall(VecGetArrayRead(v,&carr));
1012: PetscCall(VecPlaceArray(v1,carr));
1013: PetscCall(BVInsertVec(U1,i,v1));
1014: PetscCall(VecResetArray(v1));
1015: if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(U1,i,PETSC_FALSE,alpha+i,&lindep1));
1016: else { /* cheap computation of U1[i], if restart (i==k) do a full reorthogonalization */
1017: PetscCall(BVGetColumn(U1,i,&u2));
1018: if (i>0) {
1019: PetscCall(BVGetColumn(U1,i-1,&u1));
1020: PetscCall(VecAXPY(u2,-beta[i-1],u1));
1021: PetscCall(BVRestoreColumn(U1,i-1,&u1));
1022: }
1023: PetscCall(VecNorm(u2,NORM_2,&alpha[i]));
1024: if (alpha[i]==0.0) lindep = PETSC_TRUE;
1025: else PetscCall(VecScale(u2,1.0/alpha[i]));
1026: PetscCall(BVRestoreColumn(U1,i,&u2));
1027: }
1029: /* Compute vector i of BV U2 */
1030: PetscCall(BVGetColumn(U2,i,&u2));
1031: PetscCall(VecGetArray(u2,&u2arr));
1032: if (i%2) {
1033: for (j=0; j<p; j++) u2arr[j] = -carr[m+j];
1034: } else {
1035: for (j=0; j<p; j++) u2arr[j] = carr[m+j];
1036: }
1037: PetscCall(VecRestoreArray(u2,&u2arr));
1038: PetscCall(VecRestoreArrayRead(v,&carr));
1039: PetscCall(BVRestoreColumn(V,i,&v));
1040: if (lanczos->oneside && i>k) { /* cheap computation of U2[i], if restart (i==k) do a full reorthogonalization */
1041: if (i>0) {
1042: PetscCall(BVGetColumn(U2,i-1,&u1));
1043: PetscCall(VecAXPY(u2,(i%2)?betah[i-1]:-betah[i-1],u1));
1044: PetscCall(BVRestoreColumn(U2,i-1,&u1));
1045: }
1046: PetscCall(VecNorm(u2,NORM_2,&alphah[i]));
1047: if (alphah[i]==0.0) lindep = PETSC_TRUE;
1048: else PetscCall(VecScale(u2,1.0/alphah[i]));
1049: }
1050: PetscCall(BVRestoreColumn(U2,i,&u2));
1051: if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(U2,i,PETSC_FALSE,alphah+i,&lindep2));
1052: if (i%2) alphah[i] = -alphah[i];
1053: if (PetscUnlikely(lindep1 || lindep2)) {
1054: lindep = PETSC_TRUE;
1055: *n = i;
1056: break;
1057: }
1059: /* Compute vector i+1 of BV V */
1060: PetscCall(BVGetColumn(V,i+1,&v));
1061: /* Form ut=[u;0] */
1062: PetscCall(BVGetColumn(U1,i,&u));
1063: PetscCall(VecZeroEntries(ut));
1064: PetscCall(VecGetArrayRead(u,&carr));
1065: PetscCall(VecGetArray(ut,&arr));
1066: for (j=0; j<m; j++) arr[j] = carr[j];
1067: PetscCall(VecRestoreArrayRead(u,&carr));
1068: PetscCall(VecRestoreArray(ut,&arr));
1069: /* Solve least squares problem */
1070: PetscCall(KSPSolve(ksp,ut,x));
1071: PetscCall(MatMult(Z,x,v));
1072: PetscCall(BVRestoreColumn(U1,i,&u));
1073: PetscCall(BVRestoreColumn(V,i+1,&v));
1074: PetscCall(BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,beta+i,&lindep));
1075: betah[i] = -alpha[i]*beta[i]/alphah[i];
1076: if (PetscUnlikely(lindep)) {
1077: *n = i;
1078: break;
1079: }
1080: }
1081: if (breakdown) *breakdown = lindep;
1082: PetscCall(VecDestroy(&v1));
1083: PetscFunctionReturn(PETSC_SUCCESS);
1084: }
1086: /* generate random initial vector in column k for joint upper-upper bidiagonalization */
1087: static inline PetscErrorCode SVDInitialVectorGUpper(SVD svd,BV V,BV U1,PetscInt k,PetscBool *breakdown)
1088: {
1089: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1090: Vec u,v,ut=svd->workl[0],x=svd->workr[0];
1091: PetscInt m,j;
1092: PetscScalar *arr;
1093: const PetscScalar *carr;
1095: PetscFunctionBegin;
1096: /* Form ut=[u;0] where u is the k-th column of U1 */
1097: PetscCall(VecZeroEntries(ut));
1098: PetscCall(BVGetColumn(U1,k,&u));
1099: PetscCall(VecGetLocalSize(u,&m));
1100: PetscCall(VecGetArrayRead(u,&carr));
1101: PetscCall(VecGetArray(ut,&arr));
1102: for (j=0; j<m; j++) arr[j] = carr[j];
1103: PetscCall(VecRestoreArrayRead(u,&carr));
1104: PetscCall(VecRestoreArray(ut,&arr));
1105: PetscCall(BVRestoreColumn(U1,k,&u));
1106: /* Solve least squares problem Z*x=ut for x. Then set v=Z*x */
1107: PetscCall(KSPSolve(lanczos->ksp,ut,x));
1108: PetscCall(BVGetColumn(V,k,&v));
1109: PetscCall(MatMult(lanczos->Z,x,v));
1110: PetscCall(BVRestoreColumn(V,k,&v));
1111: if (breakdown) PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,NULL,breakdown));
1112: else PetscCall(BVOrthonormalizeColumn(V,k,PETSC_TRUE,NULL,NULL));
1113: PetscFunctionReturn(PETSC_SUCCESS);
1114: }
1116: /* solve generalized problem with joint upper-upper bidiagonalization */
1117: static PetscErrorCode SVDSolve_TRLanczosGUpper(SVD svd,BV U1,BV U2,BV V)
1118: {
1119: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1120: PetscReal *alpha,*beta,*alphah,*betah,normr,sigma0,*sigma;
1121: PetscScalar *w;
1122: PetscInt i,k,l,nv,ld;
1123: Mat U,Vmat,X;
1124: PetscBool breakdown=PETSC_FALSE;
1126: PetscFunctionBegin;
1127: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
1128: PetscCall(PetscMalloc2(ld,&w,ld,&sigma));
1129: normr = (svd->conv==SVD_CONV_ABS)? PetscMax(svd->nrma,svd->nrmb*lanczos->scalef): 1.0;
1131: /* normalize start vector */
1132: if (!svd->ninil) PetscCall(BVSetRandomColumn(U1,0));
1133: PetscCall(SVDInitialVectorGUpper(svd,V,U1,0,NULL));
1135: l = 0;
1136: while (svd->reason == SVD_CONVERGED_ITERATING) {
1137: svd->its++;
1139: /* inner loop */
1140: nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
1141: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
1142: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&alphah));
1143: beta = alpha + ld;
1144: betah = alpha + 2*ld;
1145: PetscCall(SVDLanczosGUpper(svd,alpha,beta,alphah,betah,lanczos->Z,U1,U2,V,lanczos->ksp,svd->nconv+l,&nv,&breakdown));
1146: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
1147: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&alphah));
1148: PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
1149: PetscCall(BVSetActiveColumns(U1,svd->nconv,nv));
1150: PetscCall(BVSetActiveColumns(U2,svd->nconv,nv));
1152: /* solve projected problem */
1153: PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l));
1154: PetscCall(DSGSVDSetDimensions(svd->ds,nv,nv));
1155: PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
1156: PetscCall(DSSolve(svd->ds,w,NULL));
1157: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
1158: PetscCall(DSUpdateExtraRow(svd->ds));
1159: PetscCall(DSSynchronize(svd->ds,w,NULL));
1160: for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
1162: /* check convergence */
1163: PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,normr,&k));
1164: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
1166: sigma0 = svd->which==SVD_LARGEST? svd->sigma[0] : 1.0/svd->sigma[0];
1167: if (lanczos->scaleth!=0 && k==0 && sigma0>lanczos->scaleth) {
1169: /* Scale and start from scratch */
1170: lanczos->scalef *= svd->sigma[0];
1171: PetscCall(PetscInfo(svd,"Scaling by factor %g and starting from scratch\n",(double)lanczos->scalef));
1172: PetscCall(MatZUpdateScale(svd));
1173: if (svd->conv==SVD_CONV_ABS) normr = PetscMax(svd->nrma,svd->nrmb*lanczos->scalef);
1174: l = 0;
1175: if (!svd->ninil) PetscCall(BVSetRandomColumn(U1,0));
1176: PetscCall(SVDInitialVectorGUpper(svd,V,U1,0,NULL));
1178: } else {
1180: /* update l */
1181: if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
1182: else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
1183: if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
1184: if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
1186: if (svd->reason == SVD_CONVERGED_ITERATING) {
1187: if (PetscUnlikely(breakdown || k==nv)) {
1188: /* Start a new bidiagonalization */
1189: PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
1190: if (k<svd->nsv) {
1191: PetscCall(BVSetRandomColumn(U1,k));
1192: PetscCall(SVDInitialVectorGUpper(svd,V,U1,k,&breakdown));
1193: if (breakdown) {
1194: svd->reason = SVD_DIVERGED_BREAKDOWN;
1195: PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
1196: }
1197: }
1198: } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
1199: }
1200: /* compute converged singular vectors and restart vectors */
1201: PetscCall(DSGetMat(svd->ds,DS_MAT_X,&X));
1202: PetscCall(BVMultInPlace(V,X,svd->nconv,k+l));
1203: PetscCall(DSRestoreMat(svd->ds,DS_MAT_X,&X));
1204: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
1205: PetscCall(BVMultInPlace(U1,U,svd->nconv,k+l));
1206: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
1207: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&Vmat));
1208: PetscCall(BVMultInPlace(U2,Vmat,svd->nconv,k+l));
1209: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&Vmat));
1211: /* copy the last vector to be the next initial vector */
1212: if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(V,nv,k+l));
1213: }
1215: svd->nconv = k;
1216: PetscCall(SVDLanczosBackTransform(svd,nv,svd->sigma,sigma,NULL));
1217: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,sigma,svd->errest,nv));
1218: }
1220: PetscCall(PetscFree2(w,sigma));
1221: PetscFunctionReturn(PETSC_SUCCESS);
1222: }
1224: /* Move generalized left singular vectors (0..nconv) from U1 and U2 to its final destination svd->U (upper and lower variants) */
1225: static inline PetscErrorCode SVDLeftSingularVectors(SVD svd,BV U1,BV U2)
1226: {
1227: PetscInt i,k,m,p;
1228: Vec u,u1,u2;
1229: PetscScalar *ua;
1230: const PetscScalar *u1a,*u2a;
1232: PetscFunctionBegin;
1233: PetscCall(BVGetSizes(U1,&m,NULL,NULL));
1234: PetscCall(BVGetSizes(U2,&p,NULL,NULL));
1235: for (i=0;i<svd->nconv;i++) {
1236: PetscCall(BVGetColumn(U1,i,&u1));
1237: PetscCall(BVGetColumn(U2,i,&u2));
1238: PetscCall(BVGetColumn(svd->U,i,&u));
1239: PetscCall(VecGetArrayRead(u1,&u1a));
1240: PetscCall(VecGetArrayRead(u2,&u2a));
1241: PetscCall(VecGetArray(u,&ua));
1242: /* Copy column from u1 to upper part of u */
1243: for (k=0;k<m;k++) ua[k] = u1a[k];
1244: /* Copy column from u2 to lower part of u */
1245: for (k=0;k<p;k++) ua[m+k] = u2a[k];
1246: PetscCall(VecRestoreArrayRead(u1,&u1a));
1247: PetscCall(VecRestoreArrayRead(u2,&u2a));
1248: PetscCall(VecRestoreArray(u,&ua));
1249: PetscCall(BVRestoreColumn(U1,i,&u1));
1250: PetscCall(BVRestoreColumn(U2,i,&u2));
1251: PetscCall(BVRestoreColumn(svd->U,i,&u));
1252: }
1253: PetscFunctionReturn(PETSC_SUCCESS);
1254: }
1256: static PetscErrorCode SVDLanczosGLower(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)
1257: {
1258: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1259: PetscInt i,j,m,p;
1260: const PetscScalar *carr;
1261: PetscScalar *arr,*u2arr;
1262: Vec u,v,ut=svd->workl[0],x=svd->workr[0],v1,u1,u2;
1263: PetscBool lindep=PETSC_FALSE;
1265: PetscFunctionBegin;
1266: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&v1));
1267: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1268: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1270: for (i=k; i<*n; i++) {
1271: /* Compute vector i of BV U2 */
1272: PetscCall(BVGetColumn(V,i,&v));
1273: PetscCall(VecGetArrayRead(v,&carr));
1274: PetscCall(BVGetColumn(U2,i,&u2));
1275: PetscCall(VecGetArray(u2,&u2arr));
1276: if (i%2) {
1277: for (j=0; j<p; j++) u2arr[j] = -carr[m+j];
1278: } else {
1279: for (j=0; j<p; j++) u2arr[j] = carr[m+j];
1280: }
1281: PetscCall(VecRestoreArray(u2,&u2arr));
1282: if (lanczos->oneside && i>k) { /* cheap computation of U2[i], if restart (i==k) do a full reorthogonalization */
1283: if (i>0) {
1284: PetscCall(BVGetColumn(U2,i-1,&u1));
1285: PetscCall(VecAXPY(u2,(i%2)?betah[i-1]:-betah[i-1],u1));
1286: PetscCall(BVRestoreColumn(U2,i-1,&u1));
1287: }
1288: PetscCall(VecNorm(u2,NORM_2,&alphah[i]));
1289: if (alphah[i]==0.0) lindep = PETSC_TRUE;
1290: else PetscCall(VecScale(u2,1.0/alphah[i]));
1291: }
1292: PetscCall(BVRestoreColumn(U2,i,&u2));
1293: if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(U2,i,PETSC_FALSE,alphah+i,&lindep));
1294: if (i%2) alphah[i] = -alphah[i];
1295: if (PetscUnlikely(lindep)) {
1296: PetscCall(BVRestoreColumn(V,i,&v));
1297: *n = i;
1298: break;
1299: }
1301: /* Compute vector i+1 of BV U1 */
1302: PetscCall(VecPlaceArray(v1,carr));
1303: PetscCall(BVInsertVec(U1,i+1,v1));
1304: PetscCall(VecResetArray(v1));
1305: PetscCall(BVOrthonormalizeColumn(U1,i+1,PETSC_FALSE,beta+i,&lindep));
1306: PetscCall(VecRestoreArrayRead(v,&carr));
1307: PetscCall(BVRestoreColumn(V,i,&v));
1308: if (PetscUnlikely(lindep)) {
1309: *n = i+1;
1310: break;
1311: }
1313: /* Compute vector i+1 of BV V */
1314: PetscCall(BVGetColumn(V,i+1,&v));
1315: /* Form ut=[u;0] where u is column i+1 of BV U1 */
1316: PetscCall(BVGetColumn(U1,i+1,&u));
1317: PetscCall(VecZeroEntries(ut));
1318: PetscCall(VecGetArrayRead(u,&carr));
1319: PetscCall(VecGetArray(ut,&arr));
1320: for (j=0; j<m; j++) arr[j] = carr[j];
1321: PetscCall(VecRestoreArrayRead(u,&carr));
1322: PetscCall(VecRestoreArray(ut,&arr));
1323: /* Solve least squares problem */
1324: PetscCall(KSPSolve(ksp,ut,x));
1325: PetscCall(MatMult(Z,x,v));
1326: PetscCall(BVRestoreColumn(U1,i+1,&u));
1327: PetscCall(BVRestoreColumn(V,i+1,&v));
1328: if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,alpha+i+1,&lindep));
1329: else { /* cheap computation of V[i+1], if restart (i==k) do a full reorthogonalization */
1330: PetscCall(BVGetColumn(V,i+1,&u2));
1331: PetscCall(BVGetColumn(V,i,&u1));
1332: PetscCall(VecAXPY(u2,-beta[i],u1));
1333: PetscCall(BVRestoreColumn(V,i,&u1));
1334: PetscCall(VecNorm(u2,NORM_2,&alpha[i+1]));
1335: if (alpha[i+1]==0.0) lindep = PETSC_TRUE;
1336: else PetscCall(VecScale(u2,1.0/alpha[i+1]));
1337: PetscCall(BVRestoreColumn(V,i+1,&u2));
1338: }
1339: betah[i] = -alpha[i+1]*beta[i]/alphah[i];
1340: if (PetscUnlikely(lindep)) {
1341: *n = i+1;
1342: break;
1343: }
1344: }
1345: if (breakdown) *breakdown = lindep;
1346: PetscCall(VecDestroy(&v1));
1347: PetscFunctionReturn(PETSC_SUCCESS);
1348: }
1350: /* generate random initial vector in column k for joint lower-upper bidiagonalization */
1351: static inline PetscErrorCode SVDInitialVectorGLower(SVD svd,BV V,BV U1,BV U2,PetscInt k,PetscBool *breakdown)
1352: {
1353: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1354: const PetscScalar *carr;
1355: PetscScalar *arr;
1356: PetscReal *alpha;
1357: PetscInt j,m,p;
1358: Vec u,uh,v,ut=svd->workl[0],x=svd->workr[0];
1360: PetscFunctionBegin;
1361: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1362: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1363: /* Form ut=[0;uh], where uh is the k-th column of U2 */
1364: PetscCall(BVGetColumn(U2,k,&uh));
1365: PetscCall(VecZeroEntries(ut));
1366: PetscCall(VecGetArrayRead(uh,&carr));
1367: PetscCall(VecGetArray(ut,&arr));
1368: for (j=0; j<p; j++) arr[m+j] = carr[j];
1369: PetscCall(VecRestoreArrayRead(uh,&carr));
1370: PetscCall(VecRestoreArray(ut,&arr));
1371: PetscCall(BVRestoreColumn(U2,k,&uh));
1372: /* Solve least squares problem Z*x=ut for x. Then set ut=Z*x */
1373: PetscCall(KSPSolve(lanczos->ksp,ut,x));
1374: PetscCall(MatMult(lanczos->Z,x,ut));
1375: /* Form u, column k of BV U1, as the upper part of ut and orthonormalize */
1376: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&u));
1377: PetscCall(VecGetArrayRead(ut,&carr));
1378: PetscCall(VecPlaceArray(u,carr));
1379: PetscCall(BVInsertVec(U1,k,u));
1380: PetscCall(VecResetArray(u));
1381: PetscCall(VecRestoreArrayRead(ut,&carr));
1382: PetscCall(VecDestroy(&u));
1383: if (breakdown) PetscCall(BVOrthonormalizeColumn(U1,k,PETSC_FALSE,NULL,breakdown));
1384: else PetscCall(BVOrthonormalizeColumn(U1,k,PETSC_TRUE,NULL,NULL));
1386: if (!breakdown || !*breakdown) {
1387: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1388: /* Compute k-th vector of BV V */
1389: PetscCall(BVGetColumn(V,k,&v));
1390: /* Form ut=[u;0] where u is the 1st column of U1 */
1391: PetscCall(BVGetColumn(U1,k,&u));
1392: PetscCall(VecZeroEntries(ut));
1393: PetscCall(VecGetArrayRead(u,&carr));
1394: PetscCall(VecGetArray(ut,&arr));
1395: for (j=0; j<m; j++) arr[j] = carr[j];
1396: PetscCall(VecRestoreArrayRead(u,&carr));
1397: PetscCall(VecRestoreArray(ut,&arr));
1398: /* Solve least squares problem */
1399: PetscCall(KSPSolve(lanczos->ksp,ut,x));
1400: PetscCall(MatMult(lanczos->Z,x,v));
1401: PetscCall(BVRestoreColumn(U1,k,&u));
1402: PetscCall(BVRestoreColumn(V,k,&v));
1403: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
1404: if (breakdown) PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,alpha+k,breakdown));
1405: else PetscCall(BVOrthonormalizeColumn(V,k,PETSC_TRUE,alpha+k,NULL));
1406: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
1407: }
1408: PetscFunctionReturn(PETSC_SUCCESS);
1409: }
1411: /* solve generalized problem with joint lower-upper bidiagonalization */
1412: static PetscErrorCode SVDSolve_TRLanczosGLower(SVD svd,BV U1,BV U2,BV V)
1413: {
1414: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1415: PetscReal *alpha,*beta,*alphah,*betah,normr,scalef,*sigma,sigma0;
1416: PetscScalar *w;
1417: PetscInt i,k,l,nv,ld;
1418: Mat U,Vmat,X;
1419: PetscBool breakdown=PETSC_FALSE,inverted;
1421: PetscFunctionBegin;
1422: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
1423: PetscCall(PetscMalloc2(ld,&w,ld,&sigma));
1424: inverted = ((svd->which==SVD_LARGEST && svd->swapped) || (svd->which==SVD_SMALLEST && !svd->swapped))? PETSC_TRUE: PETSC_FALSE;
1425: scalef = svd->swapped? 1.0/lanczos->scalef : lanczos->scalef;
1426: normr = (svd->conv==SVD_CONV_ABS)? PetscMax(svd->nrma,svd->nrmb*scalef): 1.0;
1428: /* normalize start vector */
1429: if (!svd->ninil) PetscCall(BVSetRandomColumn(U2,0));
1430: PetscCall(SVDInitialVectorGLower(svd,V,U1,U2,0,NULL));
1432: l = 0;
1433: while (svd->reason == SVD_CONVERGED_ITERATING) {
1434: svd->its++;
1436: /* inner loop */
1437: nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
1438: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
1439: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&alphah));
1440: beta = alpha + ld;
1441: betah = alpha + 2*ld;
1442: PetscCall(SVDLanczosGLower(svd,alpha,beta,alphah,betah,lanczos->Z,U1,U2,V,lanczos->ksp,svd->nconv+l,&nv,&breakdown));
1443: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
1444: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&alphah));
1445: PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
1446: PetscCall(BVSetActiveColumns(U1,svd->nconv,nv+1));
1447: PetscCall(BVSetActiveColumns(U2,svd->nconv,nv));
1449: /* solve projected problem */
1450: PetscCall(DSSetDimensions(svd->ds,nv+1,svd->nconv,svd->nconv+l));
1451: PetscCall(DSGSVDSetDimensions(svd->ds,nv,nv));
1452: PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
1453: PetscCall(DSSolve(svd->ds,w,NULL));
1454: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
1455: PetscCall(DSUpdateExtraRow(svd->ds));
1456: PetscCall(DSSynchronize(svd->ds,w,NULL));
1457: for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
1459: /* check convergence */
1460: PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,normr,&k));
1461: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
1463: sigma0 = inverted? 1.0/svd->sigma[0] : svd->sigma[0];
1464: if (lanczos->scaleth!=0 && k==0 && sigma0>lanczos->scaleth) {
1466: /* Scale and start from scratch */
1467: lanczos->scalef *= svd->swapped? 1.0/svd->sigma[0] : svd->sigma[0];
1468: PetscCall(PetscInfo(svd,"Scaling by factor %g and starting from scratch\n",(double)lanczos->scalef));
1469: PetscCall(MatZUpdateScale(svd));
1470: scalef = svd->swapped? 1.0/lanczos->scalef : lanczos->scalef;
1471: if (svd->conv==SVD_CONV_ABS) normr = PetscMax(svd->nrma,svd->nrmb*scalef);
1472: l = 0;
1473: if (!svd->ninil) PetscCall(BVSetRandomColumn(U2,0));
1474: PetscCall(SVDInitialVectorGLower(svd,V,U1,U2,0,NULL));
1476: } else {
1478: /* update l */
1479: if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
1480: else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
1481: if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
1482: if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
1484: if (svd->reason == SVD_CONVERGED_ITERATING) {
1485: if (PetscUnlikely(breakdown || k==nv)) {
1486: /* Start a new bidiagonalization */
1487: PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
1488: if (k<svd->nsv) {
1489: PetscCall(BVSetRandomColumn(U2,k));
1490: PetscCall(SVDInitialVectorGLower(svd,V,U1,U2,k,&breakdown));
1491: if (breakdown) {
1492: svd->reason = SVD_DIVERGED_BREAKDOWN;
1493: PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
1494: }
1495: }
1496: } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
1497: }
1499: /* compute converged singular vectors and restart vectors */
1500: PetscCall(DSGetMat(svd->ds,DS_MAT_X,&X));
1501: PetscCall(BVMultInPlace(V,X,svd->nconv,k+l));
1502: PetscCall(DSRestoreMat(svd->ds,DS_MAT_X,&X));
1503: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
1504: PetscCall(BVMultInPlace(U1,U,svd->nconv,k+l+1));
1505: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
1506: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&Vmat));
1507: PetscCall(BVMultInPlace(U2,Vmat,svd->nconv,k+l));
1508: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&Vmat));
1510: /* copy the last vector to be the next initial vector */
1511: if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(V,nv,k+l));
1512: }
1514: svd->nconv = k;
1515: PetscCall(SVDLanczosBackTransform(svd,nv,svd->sigma,sigma,NULL));
1516: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,sigma,svd->errest,nv));
1517: }
1519: PetscCall(PetscFree2(w,sigma));
1520: PetscFunctionReturn(PETSC_SUCCESS);
1521: }
1523: static PetscErrorCode SVDSolve_TRLanczos_GSVD(SVD svd)
1524: {
1525: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1526: PetscInt k,m,p;
1527: PetscBool convchg=PETSC_FALSE;
1528: BV U1,U2,UU;
1529: BVType type;
1530: VecType vtype;
1531: Mat U,V;
1532: SlepcSC sc;
1534: PetscFunctionBegin;
1535: PetscCall(PetscCitationsRegister(citationg,&citedg));
1537: if (svd->swapped) {
1538: PetscCall(DSGetSlepcSC(svd->ds,&sc));
1539: if (svd->which==SVD_LARGEST) sc->comparison = SlepcCompareSmallestReal;
1540: else sc->comparison = SlepcCompareLargestReal;
1541: }
1542: if (svd->converged==SVDConvergedNorm) { /* override temporarily since computed residual is already relative to the norms */
1543: svd->converged = SVDConvergedAbsolute;
1544: convchg = PETSC_TRUE;
1545: }
1546: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1547: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1549: /* Create BV for U1 */
1550: PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&U1));
1551: PetscCall(BVGetType(svd->U,&type));
1552: PetscCall(BVSetType(U1,type));
1553: PetscCall(BVGetSizes(svd->U,NULL,NULL,&k));
1554: PetscCall(BVSetSizes(U1,m,PETSC_DECIDE,k));
1555: PetscCall(BVGetVecType(svd->U,&vtype));
1556: PetscCall(BVSetVecType(U1,vtype));
1558: /* Create BV for U2 */
1559: PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&U2));
1560: PetscCall(BVSetType(U2,type));
1561: PetscCall(BVSetSizes(U2,p,PETSC_DECIDE,k));
1562: PetscCall(BVSetVecType(U2,vtype));
1564: /* Copy initial vectors from svd->U to U1 and U2 */
1565: if (svd->ninil) {
1566: Vec u, uh, nest, aux[2];
1567: PetscCall(BVGetColumn(U1,0,&u));
1568: PetscCall(BVGetColumn(U2,0,&uh));
1569: aux[0] = u;
1570: aux[1] = uh;
1571: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,aux,&nest));
1572: PetscCall(BVCopyVec(svd->U,0,nest));
1573: PetscCall(BVRestoreColumn(U1,0,&u));
1574: PetscCall(BVRestoreColumn(U2,0,&uh));
1575: PetscCall(VecDestroy(&nest));
1576: }
1578: switch (lanczos->bidiag) {
1579: case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1580: PetscCall(SVDSolve_TRLanczosGSingle(svd,U1,svd->U));
1581: break;
1582: case SVD_TRLANCZOS_GBIDIAG_UPPER:
1583: PetscCall(SVDSolve_TRLanczosGUpper(svd,U1,U2,svd->U));
1584: break;
1585: case SVD_TRLANCZOS_GBIDIAG_LOWER:
1586: PetscCall(SVDSolve_TRLanczosGLower(svd,U1,U2,svd->U));
1587: break;
1588: }
1590: /* Compute converged right singular vectors */
1591: PetscCall(BVSetActiveColumns(svd->U,0,svd->nconv));
1592: PetscCall(BVSetActiveColumns(svd->V,0,svd->nconv));
1593: PetscCall(BVGetMat(svd->U,&U));
1594: PetscCall(BVGetMat(svd->V,&V));
1595: PetscCall(KSPMatSolve(lanczos->ksp,U,V));
1596: PetscCall(BVRestoreMat(svd->U,&U));
1597: PetscCall(BVRestoreMat(svd->V,&V));
1599: /* Finish computing left singular vectors and move them to its place */
1600: if (svd->swapped) SlepcSwap(U1,U2,UU);
1601: switch (lanczos->bidiag) {
1602: case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1603: PetscCall(SVDLeftSingularVectors_Single(svd,U1,U2));
1604: break;
1605: case SVD_TRLANCZOS_GBIDIAG_UPPER:
1606: case SVD_TRLANCZOS_GBIDIAG_LOWER:
1607: PetscCall(SVDLeftSingularVectors(svd,U1,U2));
1608: break;
1609: }
1611: /* undo scaling and compute the reciprocals of sigma if matrices were swapped */
1612: PetscCall(SVDLanczosBackTransform(svd,svd->nconv,svd->sigma,NULL,svd->V));
1614: PetscCall(BVDestroy(&U1));
1615: PetscCall(BVDestroy(&U2));
1616: PetscCall(DSTruncate(svd->ds,svd->nconv,PETSC_TRUE));
1617: if (convchg) svd->converged = SVDConvergedNorm;
1618: PetscFunctionReturn(PETSC_SUCCESS);
1619: }
1621: static PetscErrorCode SVDSetFromOptions_TRLanczos(SVD svd,PetscOptionItems *PetscOptionsObject)
1622: {
1623: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1624: PetscBool flg,val,lock;
1625: PetscReal keep,scale;
1626: SVDTRLanczosGBidiag bidiag;
1628: PetscFunctionBegin;
1629: PetscOptionsHeadBegin(PetscOptionsObject,"SVD TRLanczos Options");
1631: PetscCall(PetscOptionsBool("-svd_trlanczos_oneside","Use one-side reorthogonalization","SVDTRLanczosSetOneSide",lanczos->oneside,&val,&flg));
1632: if (flg) PetscCall(SVDTRLanczosSetOneSide(svd,val));
1634: PetscCall(PetscOptionsReal("-svd_trlanczos_restart","Proportion of vectors kept after restart","SVDTRLanczosSetRestart",0.5,&keep,&flg));
1635: if (flg) PetscCall(SVDTRLanczosSetRestart(svd,keep));
1637: PetscCall(PetscOptionsBool("-svd_trlanczos_locking","Choose between locking and non-locking variants","SVDTRLanczosSetLocking",PETSC_TRUE,&lock,&flg));
1638: if (flg) PetscCall(SVDTRLanczosSetLocking(svd,lock));
1640: PetscCall(PetscOptionsEnum("-svd_trlanczos_gbidiag","Bidiagonalization choice for Generalized Problem","SVDTRLanczosSetGBidiag",SVDTRLanczosGBidiags,(PetscEnum)lanczos->bidiag,(PetscEnum*)&bidiag,&flg));
1641: if (flg) PetscCall(SVDTRLanczosSetGBidiag(svd,bidiag));
1643: PetscCall(PetscOptionsBool("-svd_trlanczos_explicitmatrix","Build explicit matrix for KSP solver","SVDTRLanczosSetExplicitMatrix",lanczos->explicitmatrix,&val,&flg));
1644: if (flg) PetscCall(SVDTRLanczosSetExplicitMatrix(svd,val));
1646: PetscCall(SVDTRLanczosGetScale(svd,&scale));
1647: PetscCall(PetscOptionsReal("-svd_trlanczos_scale","Scale parameter for matrix B","SVDTRLanczosSetScale",scale,&scale,&flg));
1648: if (flg) PetscCall(SVDTRLanczosSetScale(svd,scale));
1650: PetscOptionsHeadEnd();
1652: if (svd->OPb) {
1653: if (!lanczos->ksp) PetscCall(SVDTRLanczosGetKSP(svd,&lanczos->ksp));
1654: PetscCall(KSPSetFromOptions(lanczos->ksp));
1655: }
1656: PetscFunctionReturn(PETSC_SUCCESS);
1657: }
1659: static PetscErrorCode SVDTRLanczosSetOneSide_TRLanczos(SVD svd,PetscBool oneside)
1660: {
1661: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1663: PetscFunctionBegin;
1664: if (lanczos->oneside != oneside) {
1665: lanczos->oneside = oneside;
1666: svd->state = SVD_STATE_INITIAL;
1667: }
1668: PetscFunctionReturn(PETSC_SUCCESS);
1669: }
1671: /*@
1672: SVDTRLanczosSetOneSide - Indicate if the variant of the Lanczos method
1673: to be used is one-sided or two-sided.
1675: Logically Collective
1677: Input Parameters:
1678: + svd - singular value solver
1679: - oneside - boolean flag indicating if the method is one-sided or not
1681: Options Database Key:
1682: . -svd_trlanczos_oneside <boolean> - Indicates the boolean flag
1684: Notes:
1685: By default, a two-sided variant is selected, which is sometimes slightly
1686: more robust. However, the one-sided variant is faster because it avoids
1687: the orthogonalization associated to left singular vectors.
1689: One-sided orthogonalization is also available for the GSVD, in which case
1690: two orthogonalizations out of three are avoided.
1692: Level: advanced
1694: .seealso: SVDLanczosSetOneSide()
1695: @*/
1696: PetscErrorCode SVDTRLanczosSetOneSide(SVD svd,PetscBool oneside)
1697: {
1698: PetscFunctionBegin;
1701: PetscTryMethod(svd,"SVDTRLanczosSetOneSide_C",(SVD,PetscBool),(svd,oneside));
1702: PetscFunctionReturn(PETSC_SUCCESS);
1703: }
1705: static PetscErrorCode SVDTRLanczosGetOneSide_TRLanczos(SVD svd,PetscBool *oneside)
1706: {
1707: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1709: PetscFunctionBegin;
1710: *oneside = lanczos->oneside;
1711: PetscFunctionReturn(PETSC_SUCCESS);
1712: }
1714: /*@
1715: SVDTRLanczosGetOneSide - Gets if the variant of the Lanczos method
1716: to be used is one-sided or two-sided.
1718: Not Collective
1720: Input Parameters:
1721: . svd - singular value solver
1723: Output Parameters:
1724: . oneside - boolean flag indicating if the method is one-sided or not
1726: Level: advanced
1728: .seealso: SVDTRLanczosSetOneSide()
1729: @*/
1730: PetscErrorCode SVDTRLanczosGetOneSide(SVD svd,PetscBool *oneside)
1731: {
1732: PetscFunctionBegin;
1734: PetscAssertPointer(oneside,2);
1735: PetscUseMethod(svd,"SVDTRLanczosGetOneSide_C",(SVD,PetscBool*),(svd,oneside));
1736: PetscFunctionReturn(PETSC_SUCCESS);
1737: }
1739: static PetscErrorCode SVDTRLanczosSetGBidiag_TRLanczos(SVD svd,SVDTRLanczosGBidiag bidiag)
1740: {
1741: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1743: PetscFunctionBegin;
1744: switch (bidiag) {
1745: case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1746: case SVD_TRLANCZOS_GBIDIAG_UPPER:
1747: case SVD_TRLANCZOS_GBIDIAG_LOWER:
1748: if (lanczos->bidiag != bidiag) {
1749: lanczos->bidiag = bidiag;
1750: svd->state = SVD_STATE_INITIAL;
1751: }
1752: break;
1753: default:
1754: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid bidiagonalization choice");
1755: }
1756: PetscFunctionReturn(PETSC_SUCCESS);
1757: }
1759: /*@
1760: SVDTRLanczosSetGBidiag - Sets the bidiagonalization choice to use in
1761: the GSVD TRLanczos solver.
1763: Logically Collective
1765: Input Parameters:
1766: + svd - the singular value solver
1767: - bidiag - the bidiagonalization choice
1769: Options Database Key:
1770: . -svd_trlanczos_gbidiag - Sets the bidiagonalization choice (either 's' or 'juu'
1771: or 'jlu')
1773: Level: advanced
1775: .seealso: SVDTRLanczosGetGBidiag(), SVDTRLanczosGBidiag
1776: @*/
1777: PetscErrorCode SVDTRLanczosSetGBidiag(SVD svd,SVDTRLanczosGBidiag bidiag)
1778: {
1779: PetscFunctionBegin;
1782: PetscTryMethod(svd,"SVDTRLanczosSetGBidiag_C",(SVD,SVDTRLanczosGBidiag),(svd,bidiag));
1783: PetscFunctionReturn(PETSC_SUCCESS);
1784: }
1786: static PetscErrorCode SVDTRLanczosGetGBidiag_TRLanczos(SVD svd,SVDTRLanczosGBidiag *bidiag)
1787: {
1788: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1790: PetscFunctionBegin;
1791: *bidiag = lanczos->bidiag;
1792: PetscFunctionReturn(PETSC_SUCCESS);
1793: }
1795: /*@
1796: SVDTRLanczosGetGBidiag - Gets the bidiagonalization choice used in the GSVD
1797: TRLanczos solver.
1799: Not Collective
1801: Input Parameter:
1802: . svd - the singular value solver
1804: Output Parameter:
1805: . bidiag - the bidiagonalization choice
1807: Level: advanced
1809: .seealso: SVDTRLanczosSetGBidiag(), SVDTRLanczosGBidiag
1810: @*/
1811: PetscErrorCode SVDTRLanczosGetGBidiag(SVD svd,SVDTRLanczosGBidiag *bidiag)
1812: {
1813: PetscFunctionBegin;
1815: PetscAssertPointer(bidiag,2);
1816: PetscUseMethod(svd,"SVDTRLanczosGetGBidiag_C",(SVD,SVDTRLanczosGBidiag*),(svd,bidiag));
1817: PetscFunctionReturn(PETSC_SUCCESS);
1818: }
1820: static PetscErrorCode SVDTRLanczosSetKSP_TRLanczos(SVD svd,KSP ksp)
1821: {
1822: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1824: PetscFunctionBegin;
1825: PetscCall(PetscObjectReference((PetscObject)ksp));
1826: PetscCall(KSPDestroy(&ctx->ksp));
1827: ctx->ksp = ksp;
1828: svd->state = SVD_STATE_INITIAL;
1829: PetscFunctionReturn(PETSC_SUCCESS);
1830: }
1832: /*@
1833: SVDTRLanczosSetKSP - Associate a linear solver object (KSP) to the SVD solver.
1835: Collective
1837: Input Parameters:
1838: + svd - SVD solver
1839: - ksp - the linear solver object
1841: Note:
1842: Only used for the GSVD problem.
1844: Level: advanced
1846: .seealso: SVDTRLanczosGetKSP()
1847: @*/
1848: PetscErrorCode SVDTRLanczosSetKSP(SVD svd,KSP ksp)
1849: {
1850: PetscFunctionBegin;
1853: PetscCheckSameComm(svd,1,ksp,2);
1854: PetscTryMethod(svd,"SVDTRLanczosSetKSP_C",(SVD,KSP),(svd,ksp));
1855: PetscFunctionReturn(PETSC_SUCCESS);
1856: }
1858: static PetscErrorCode SVDTRLanczosGetKSP_TRLanczos(SVD svd,KSP *ksp)
1859: {
1860: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1861: PC pc;
1863: PetscFunctionBegin;
1864: if (!ctx->ksp) {
1865: /* Create linear solver */
1866: PetscCall(KSPCreate(PetscObjectComm((PetscObject)svd),&ctx->ksp));
1867: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)svd,1));
1868: PetscCall(KSPSetOptionsPrefix(ctx->ksp,((PetscObject)svd)->prefix));
1869: PetscCall(KSPAppendOptionsPrefix(ctx->ksp,"svd_trlanczos_"));
1870: PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)svd)->options));
1871: PetscCall(KSPSetType(ctx->ksp,KSPLSQR));
1872: PetscCall(KSPGetPC(ctx->ksp,&pc));
1873: PetscCall(PCSetType(pc,PCNONE));
1874: PetscCall(KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE));
1875: PetscCall(KSPSetTolerances(ctx->ksp,SlepcDefaultTol(svd->tol)/10.0,PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
1876: }
1877: *ksp = ctx->ksp;
1878: PetscFunctionReturn(PETSC_SUCCESS);
1879: }
1881: /*@
1882: SVDTRLanczosGetKSP - Retrieve the linear solver object (KSP) associated with
1883: the SVD solver.
1885: Collective
1887: Input Parameter:
1888: . svd - SVD solver
1890: Output Parameter:
1891: . ksp - the linear solver object
1893: Level: advanced
1895: .seealso: SVDTRLanczosSetKSP()
1896: @*/
1897: PetscErrorCode SVDTRLanczosGetKSP(SVD svd,KSP *ksp)
1898: {
1899: PetscFunctionBegin;
1901: PetscAssertPointer(ksp,2);
1902: PetscUseMethod(svd,"SVDTRLanczosGetKSP_C",(SVD,KSP*),(svd,ksp));
1903: PetscFunctionReturn(PETSC_SUCCESS);
1904: }
1906: static PetscErrorCode SVDTRLanczosSetRestart_TRLanczos(SVD svd,PetscReal keep)
1907: {
1908: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1910: PetscFunctionBegin;
1911: if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) ctx->keep = 0.5;
1912: else {
1913: PetscCheck(keep>=0.1 && keep<=0.9,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument %g must be in the range [0.1,0.9]",(double)keep);
1914: ctx->keep = keep;
1915: }
1916: PetscFunctionReturn(PETSC_SUCCESS);
1917: }
1919: /*@
1920: SVDTRLanczosSetRestart - Sets the restart parameter for the thick-restart
1921: Lanczos method, in particular the proportion of basis vectors that must be
1922: kept after restart.
1924: Logically Collective
1926: Input Parameters:
1927: + svd - the singular value solver
1928: - keep - the number of vectors to be kept at restart
1930: Options Database Key:
1931: . -svd_trlanczos_restart - Sets the restart parameter
1933: Notes:
1934: Allowed values are in the range [0.1,0.9]. The default is 0.5.
1936: Level: advanced
1938: .seealso: SVDTRLanczosGetRestart()
1939: @*/
1940: PetscErrorCode SVDTRLanczosSetRestart(SVD svd,PetscReal keep)
1941: {
1942: PetscFunctionBegin;
1945: PetscTryMethod(svd,"SVDTRLanczosSetRestart_C",(SVD,PetscReal),(svd,keep));
1946: PetscFunctionReturn(PETSC_SUCCESS);
1947: }
1949: static PetscErrorCode SVDTRLanczosGetRestart_TRLanczos(SVD svd,PetscReal *keep)
1950: {
1951: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1953: PetscFunctionBegin;
1954: *keep = ctx->keep;
1955: PetscFunctionReturn(PETSC_SUCCESS);
1956: }
1958: /*@
1959: SVDTRLanczosGetRestart - Gets the restart parameter used in the thick-restart
1960: Lanczos method.
1962: Not Collective
1964: Input Parameter:
1965: . svd - the singular value solver
1967: Output Parameter:
1968: . keep - the restart parameter
1970: Level: advanced
1972: .seealso: SVDTRLanczosSetRestart()
1973: @*/
1974: PetscErrorCode SVDTRLanczosGetRestart(SVD svd,PetscReal *keep)
1975: {
1976: PetscFunctionBegin;
1978: PetscAssertPointer(keep,2);
1979: PetscUseMethod(svd,"SVDTRLanczosGetRestart_C",(SVD,PetscReal*),(svd,keep));
1980: PetscFunctionReturn(PETSC_SUCCESS);
1981: }
1983: static PetscErrorCode SVDTRLanczosSetLocking_TRLanczos(SVD svd,PetscBool lock)
1984: {
1985: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1987: PetscFunctionBegin;
1988: ctx->lock = lock;
1989: PetscFunctionReturn(PETSC_SUCCESS);
1990: }
1992: /*@
1993: SVDTRLanczosSetLocking - Choose between locking and non-locking variants of
1994: the thick-restart Lanczos method.
1996: Logically Collective
1998: Input Parameters:
1999: + svd - the singular value solver
2000: - lock - true if the locking variant must be selected
2002: Options Database Key:
2003: . -svd_trlanczos_locking - Sets the locking flag
2005: Notes:
2006: The default is to lock converged singular triplets when the method restarts.
2007: This behaviour can be changed so that all directions are kept in the
2008: working subspace even if already converged to working accuracy (the
2009: non-locking variant).
2011: Level: advanced
2013: .seealso: SVDTRLanczosGetLocking()
2014: @*/
2015: PetscErrorCode SVDTRLanczosSetLocking(SVD svd,PetscBool lock)
2016: {
2017: PetscFunctionBegin;
2020: PetscTryMethod(svd,"SVDTRLanczosSetLocking_C",(SVD,PetscBool),(svd,lock));
2021: PetscFunctionReturn(PETSC_SUCCESS);
2022: }
2024: static PetscErrorCode SVDTRLanczosGetLocking_TRLanczos(SVD svd,PetscBool *lock)
2025: {
2026: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
2028: PetscFunctionBegin;
2029: *lock = ctx->lock;
2030: PetscFunctionReturn(PETSC_SUCCESS);
2031: }
2033: /*@
2034: SVDTRLanczosGetLocking - Gets the locking flag used in the thick-restart
2035: Lanczos method.
2037: Not Collective
2039: Input Parameter:
2040: . svd - the singular value solver
2042: Output Parameter:
2043: . lock - the locking flag
2045: Level: advanced
2047: .seealso: SVDTRLanczosSetLocking()
2048: @*/
2049: PetscErrorCode SVDTRLanczosGetLocking(SVD svd,PetscBool *lock)
2050: {
2051: PetscFunctionBegin;
2053: PetscAssertPointer(lock,2);
2054: PetscUseMethod(svd,"SVDTRLanczosGetLocking_C",(SVD,PetscBool*),(svd,lock));
2055: PetscFunctionReturn(PETSC_SUCCESS);
2056: }
2058: static PetscErrorCode SVDTRLanczosSetExplicitMatrix_TRLanczos(SVD svd,PetscBool explicitmat)
2059: {
2060: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2062: PetscFunctionBegin;
2063: if (lanczos->explicitmatrix != explicitmat) {
2064: lanczos->explicitmatrix = explicitmat;
2065: svd->state = SVD_STATE_INITIAL;
2066: }
2067: PetscFunctionReturn(PETSC_SUCCESS);
2068: }
2070: /*@
2071: SVDTRLanczosSetExplicitMatrix - Indicate if the matrix Z=[A;B] must
2072: be built explicitly.
2074: Logically Collective
2076: Input Parameters:
2077: + svd - singular value solver
2078: - explicitmat - Boolean flag indicating if Z=[A;B] is built explicitly
2080: Options Database Key:
2081: . -svd_trlanczos_explicitmatrix <boolean> - Indicates the boolean flag
2083: Notes:
2084: This option is relevant for the GSVD case only.
2085: Z is the coefficient matrix of the KSP solver used internally.
2087: Level: advanced
2089: .seealso: SVDTRLanczosGetExplicitMatrix()
2090: @*/
2091: PetscErrorCode SVDTRLanczosSetExplicitMatrix(SVD svd,PetscBool explicitmat)
2092: {
2093: PetscFunctionBegin;
2096: PetscTryMethod(svd,"SVDTRLanczosSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
2097: PetscFunctionReturn(PETSC_SUCCESS);
2098: }
2100: static PetscErrorCode SVDTRLanczosGetExplicitMatrix_TRLanczos(SVD svd,PetscBool *explicitmat)
2101: {
2102: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2104: PetscFunctionBegin;
2105: *explicitmat = lanczos->explicitmatrix;
2106: PetscFunctionReturn(PETSC_SUCCESS);
2107: }
2109: /*@
2110: SVDTRLanczosGetExplicitMatrix - Returns the flag indicating if Z=[A;B] is built explicitly.
2112: Not Collective
2114: Input Parameter:
2115: . svd - singular value solver
2117: Output Parameter:
2118: . explicitmat - the mode flag
2120: Level: advanced
2122: .seealso: SVDTRLanczosSetExplicitMatrix()
2123: @*/
2124: PetscErrorCode SVDTRLanczosGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
2125: {
2126: PetscFunctionBegin;
2128: PetscAssertPointer(explicitmat,2);
2129: PetscUseMethod(svd,"SVDTRLanczosGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
2130: PetscFunctionReturn(PETSC_SUCCESS);
2131: }
2133: static PetscErrorCode SVDTRLanczosSetScale_TRLanczos(SVD svd,PetscReal scale)
2134: {
2135: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
2137: PetscFunctionBegin;
2138: if (scale<0) {
2139: ctx->scalef = 1.0;
2140: ctx->scaleth = -scale;
2141: } else {
2142: ctx->scalef = scale;
2143: ctx->scaleth = 0.0;
2144: }
2145: PetscFunctionReturn(PETSC_SUCCESS);
2146: }
2148: /*@
2149: SVDTRLanczosSetScale - Sets the scale parameter for the GSVD.
2151: Logically Collective
2153: Input Parameters:
2154: + svd - singular value solver
2155: - scale - scale parameter
2157: Options Database Key:
2158: . -svd_trlanczos_scale <real> - scale factor/threshold
2160: Notes:
2161: This parameter is relevant for the GSVD case only. If the parameter is
2162: positive, it indicates the scale factor for B in matrix Z=[A;B]. If
2163: negative, its absolute value is the threshold for automatic scaling.
2164: In automatic scaling, whenever the largest approximate generalized singular
2165: value (or the inverse of the smallest value, if SVD_SMALLEST is used)
2166: exceeds the threshold, the computation is restarted with matrix B
2167: scaled by that value.
2169: Level: advanced
2171: .seealso: SVDTRLanczosGetScale()
2172: @*/
2173: PetscErrorCode SVDTRLanczosSetScale(SVD svd,PetscReal scale)
2174: {
2175: PetscFunctionBegin;
2178: PetscTryMethod(svd,"SVDTRLanczosSetScale_C",(SVD,PetscReal),(svd,scale));
2179: PetscFunctionReturn(PETSC_SUCCESS);
2180: }
2182: static PetscErrorCode SVDTRLanczosGetScale_TRLanczos(SVD svd,PetscReal *scale)
2183: {
2184: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
2186: PetscFunctionBegin;
2187: if (ctx->scaleth==0) *scale = ctx->scalef;
2188: else *scale = -ctx->scaleth;
2189: PetscFunctionReturn(PETSC_SUCCESS);
2190: }
2192: /*@
2193: SVDTRLanczosGetScale - Gets the scale parameter for the GSVD.
2195: Not Collective
2197: Input Parameter:
2198: . svd - the singular value solver
2200: Output Parameter:
2201: . scale - the scale parameter
2203: Notes:
2204: This parameter is relevant for the GSVD case only. If the parameter is
2205: positive, it indicates the scale factor for B in matrix Z=[A;B]. If
2206: negative, its absolute value is the threshold for automatic scaling.
2208: Level: advanced
2210: .seealso: SVDTRLanczosSetScale()
2211: @*/
2212: PetscErrorCode SVDTRLanczosGetScale(SVD svd,PetscReal *scale)
2213: {
2214: PetscFunctionBegin;
2216: PetscAssertPointer(scale,2);
2217: PetscUseMethod(svd,"SVDTRLanczosGetScale_C",(SVD,PetscReal*),(svd,scale));
2218: PetscFunctionReturn(PETSC_SUCCESS);
2219: }
2221: static PetscErrorCode SVDReset_TRLanczos(SVD svd)
2222: {
2223: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2225: PetscFunctionBegin;
2226: if (svd->isgeneralized || (!svd->problem_type && svd->OPb)) {
2227: PetscCall(KSPReset(lanczos->ksp));
2228: PetscCall(MatDestroy(&lanczos->Z));
2229: }
2230: PetscFunctionReturn(PETSC_SUCCESS);
2231: }
2233: static PetscErrorCode SVDDestroy_TRLanczos(SVD svd)
2234: {
2235: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2237: PetscFunctionBegin;
2238: if (svd->isgeneralized || (!svd->problem_type && svd->OPb)) PetscCall(KSPDestroy(&lanczos->ksp));
2239: PetscCall(PetscFree(svd->data));
2240: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",NULL));
2241: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",NULL));
2242: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetGBidiag_C",NULL));
2243: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetGBidiag_C",NULL));
2244: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetKSP_C",NULL));
2245: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetKSP_C",NULL));
2246: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetRestart_C",NULL));
2247: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetRestart_C",NULL));
2248: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetLocking_C",NULL));
2249: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetLocking_C",NULL));
2250: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetExplicitMatrix_C",NULL));
2251: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetExplicitMatrix_C",NULL));
2252: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetScale_C",NULL));
2253: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetScale_C",NULL));
2254: PetscFunctionReturn(PETSC_SUCCESS);
2255: }
2257: static PetscErrorCode SVDView_TRLanczos(SVD svd,PetscViewer viewer)
2258: {
2259: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2260: PetscBool isascii;
2262: PetscFunctionBegin;
2263: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
2264: if (isascii) {
2265: PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*lanczos->keep)));
2266: PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",lanczos->lock?"":"non-"));
2267: if (svd->isgeneralized) {
2268: const char *bidiag="";
2270: switch (lanczos->bidiag) {
2271: case SVD_TRLANCZOS_GBIDIAG_SINGLE: bidiag = "single"; break;
2272: case SVD_TRLANCZOS_GBIDIAG_UPPER: bidiag = "joint upper-upper"; break;
2273: case SVD_TRLANCZOS_GBIDIAG_LOWER: bidiag = "joint lower-upper"; break;
2274: }
2275: PetscCall(PetscViewerASCIIPrintf(viewer," bidiagonalization choice: %s\n",bidiag));
2276: PetscCall(PetscViewerASCIIPrintf(viewer," %s matrix\n",lanczos->explicitmatrix?"explicit":"implicit"));
2277: if (lanczos->scaleth==0) PetscCall(PetscViewerASCIIPrintf(viewer," scale factor for matrix B: %g\n",(double)lanczos->scalef));
2278: else PetscCall(PetscViewerASCIIPrintf(viewer," automatic scaling for matrix B with threshold: %g\n",(double)lanczos->scaleth));
2279: if (!lanczos->ksp) PetscCall(SVDTRLanczosGetKSP(svd,&lanczos->ksp));
2280: PetscCall(PetscViewerASCIIPushTab(viewer));
2281: PetscCall(KSPView(lanczos->ksp,viewer));
2282: PetscCall(PetscViewerASCIIPopTab(viewer));
2283: } else PetscCall(PetscViewerASCIIPrintf(viewer," %s-sided reorthogonalization\n",lanczos->oneside? "one": "two"));
2284: }
2285: PetscFunctionReturn(PETSC_SUCCESS);
2286: }
2288: static PetscErrorCode SVDSetDSType_TRLanczos(SVD svd)
2289: {
2290: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2291: DSType dstype;
2293: PetscFunctionBegin;
2294: dstype = svd->ishyperbolic? DSHSVD: DSSVD;
2295: if (svd->OPb && (lanczos->bidiag==SVD_TRLANCZOS_GBIDIAG_UPPER || lanczos->bidiag==SVD_TRLANCZOS_GBIDIAG_LOWER)) dstype = DSGSVD;
2296: PetscCall(DSSetType(svd->ds,dstype));
2297: PetscFunctionReturn(PETSC_SUCCESS);
2298: }
2300: SLEPC_EXTERN PetscErrorCode SVDCreate_TRLanczos(SVD svd)
2301: {
2302: SVD_TRLANCZOS *ctx;
2304: PetscFunctionBegin;
2305: PetscCall(PetscNew(&ctx));
2306: svd->data = (void*)ctx;
2308: ctx->lock = PETSC_TRUE;
2309: ctx->bidiag = SVD_TRLANCZOS_GBIDIAG_LOWER;
2310: ctx->scalef = 1.0;
2311: ctx->scaleth = 0.0;
2313: svd->ops->setup = SVDSetUp_TRLanczos;
2314: svd->ops->solve = SVDSolve_TRLanczos;
2315: svd->ops->solveg = SVDSolve_TRLanczos_GSVD;
2316: svd->ops->solveh = SVDSolve_TRLanczos_HSVD;
2317: svd->ops->destroy = SVDDestroy_TRLanczos;
2318: svd->ops->reset = SVDReset_TRLanczos;
2319: svd->ops->setfromoptions = SVDSetFromOptions_TRLanczos;
2320: svd->ops->view = SVDView_TRLanczos;
2321: svd->ops->setdstype = SVDSetDSType_TRLanczos;
2322: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",SVDTRLanczosSetOneSide_TRLanczos));
2323: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",SVDTRLanczosGetOneSide_TRLanczos));
2324: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetGBidiag_C",SVDTRLanczosSetGBidiag_TRLanczos));
2325: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetGBidiag_C",SVDTRLanczosGetGBidiag_TRLanczos));
2326: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetKSP_C",SVDTRLanczosSetKSP_TRLanczos));
2327: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetKSP_C",SVDTRLanczosGetKSP_TRLanczos));
2328: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetRestart_C",SVDTRLanczosSetRestart_TRLanczos));
2329: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetRestart_C",SVDTRLanczosGetRestart_TRLanczos));
2330: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetLocking_C",SVDTRLanczosSetLocking_TRLanczos));
2331: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetLocking_C",SVDTRLanczosGetLocking_TRLanczos));
2332: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetExplicitMatrix_C",SVDTRLanczosSetExplicitMatrix_TRLanczos));
2333: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetExplicitMatrix_C",SVDTRLanczosGetExplicitMatrix_TRLanczos));
2334: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetScale_C",SVDTRLanczosSetScale_TRLanczos));
2335: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetScale_C",SVDTRLanczosGetScale_TRLanczos));
2336: PetscFunctionReturn(PETSC_SUCCESS);
2337: }