Actual source code: trlanczos.c
slepc-3.21.0 2024-03-30
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: #define SWAP(a,b,t) do {t=a;a=b;b=t;} while (0)
205: static PetscErrorCode SVDSetUp_TRLanczos(SVD svd)
206: {
207: PetscInt M,N,P,m,n,p;
208: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
209: MatZData *zdata;
210: Mat aux;
212: PetscFunctionBegin;
213: PetscCall(MatGetSize(svd->A,&M,&N));
214: PetscCall(SVDSetDimensions_Default(svd));
215: 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");
216: PetscCheck(lanczos->lock || svd->mpd>=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
217: if (svd->max_it==PETSC_DEFAULT) svd->max_it = PetscMax(N/svd->ncv,100);
218: if (!lanczos->keep) lanczos->keep = 0.5;
219: svd->leftbasis = PETSC_TRUE;
220: PetscCall(SVDAllocateSolution(svd,1));
221: if (svd->isgeneralized) {
222: PetscCall(MatGetSize(svd->B,&P,NULL));
223: if (lanczos->bidiag == SVD_TRLANCZOS_GBIDIAG_LOWER && ((svd->which==SVD_LARGEST && P<=N) || (svd->which==SVD_SMALLEST && M>N && P<=N))) {
224: SWAP(svd->A,svd->B,aux);
225: SWAP(svd->AT,svd->BT,aux);
226: svd->swapped = PETSC_TRUE;
227: } else svd->swapped = PETSC_FALSE;
229: PetscCall(SVDSetWorkVecs(svd,1,1));
231: if (svd->conv==SVD_CONV_ABS) { /* Residual norms are multiplied by matrix norms */
232: if (!svd->nrma) PetscCall(MatNorm(svd->A,NORM_INFINITY,&svd->nrma));
233: if (!svd->nrmb) PetscCall(MatNorm(svd->B,NORM_INFINITY,&svd->nrmb));
234: }
236: /* Create the matrix Z=[A;B] */
237: PetscCall(MatGetLocalSize(svd->A,&m,&n));
238: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
239: if (!lanczos->explicitmatrix) {
240: PetscCall(MatDestroy(&lanczos->Z));
241: PetscCall(MatZCreateContext(svd,&zdata));
242: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)svd),m+p,n,PETSC_DECIDE,PETSC_DECIDE,zdata,&lanczos->Z));
243: PetscCall(MatShellSetOperation(lanczos->Z,MATOP_MULT,(void(*)(void))MatMult_Z));
244: #if defined(PETSC_USE_COMPLEX)
245: PetscCall(MatShellSetOperation(lanczos->Z,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultTranspose_Z));
246: #else
247: PetscCall(MatShellSetOperation(lanczos->Z,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Z));
248: #endif
249: PetscCall(MatShellSetOperation(lanczos->Z,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_Z));
250: PetscCall(MatShellSetOperation(lanczos->Z,MATOP_DESTROY,(void(*)(void))MatDestroy_Z));
251: }
252: /* Explicit matrix is created here, when updating the scale */
253: PetscCall(MatZUpdateScale(svd));
255: } else if (svd->ishyperbolic) {
256: PetscCall(BV_SetMatrixDiagonal(svd->swapped?svd->V:svd->U,svd->omega,svd->OP));
257: PetscCall(SVDSetWorkVecs(svd,1,0));
258: }
259: PetscCall(DSSetCompact(svd->ds,PETSC_TRUE));
260: PetscCall(DSSetExtraRow(svd->ds,PETSC_TRUE));
261: PetscCall(DSAllocate(svd->ds,svd->ncv+1));
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: static PetscErrorCode SVDOneSideTRLanczosMGS(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,BV U,PetscInt nconv,PetscInt l,PetscInt n,PetscScalar* work)
266: {
267: PetscReal a,b;
268: PetscInt i,k=nconv+l;
269: Vec ui,ui1,vi;
271: PetscFunctionBegin;
272: PetscCall(BVGetColumn(V,k,&vi));
273: PetscCall(BVGetColumn(U,k,&ui));
274: PetscCall(MatMult(svd->A,vi,ui));
275: PetscCall(BVRestoreColumn(V,k,&vi));
276: PetscCall(BVRestoreColumn(U,k,&ui));
277: if (l>0) {
278: PetscCall(BVSetActiveColumns(U,nconv,n));
279: for (i=0;i<l;i++) work[i]=beta[i+nconv];
280: PetscCall(BVMultColumn(U,-1.0,1.0,k,work));
281: }
282: PetscCall(BVNormColumn(U,k,NORM_2,&a));
283: PetscCall(BVScaleColumn(U,k,1.0/a));
284: alpha[k] = a;
286: for (i=k+1;i<n;i++) {
287: PetscCall(BVGetColumn(V,i,&vi));
288: PetscCall(BVGetColumn(U,i-1,&ui1));
289: PetscCall(MatMult(svd->AT,ui1,vi));
290: PetscCall(BVRestoreColumn(V,i,&vi));
291: PetscCall(BVRestoreColumn(U,i-1,&ui1));
292: PetscCall(BVOrthonormalizeColumn(V,i,PETSC_FALSE,&b,NULL));
293: beta[i-1] = b;
295: PetscCall(BVGetColumn(V,i,&vi));
296: PetscCall(BVGetColumn(U,i,&ui));
297: PetscCall(MatMult(svd->A,vi,ui));
298: PetscCall(BVRestoreColumn(V,i,&vi));
299: PetscCall(BVGetColumn(U,i-1,&ui1));
300: PetscCall(VecAXPY(ui,-b,ui1));
301: PetscCall(BVRestoreColumn(U,i-1,&ui1));
302: PetscCall(BVRestoreColumn(U,i,&ui));
303: PetscCall(BVNormColumn(U,i,NORM_2,&a));
304: PetscCall(BVScaleColumn(U,i,1.0/a));
305: alpha[i] = a;
306: }
308: PetscCall(BVGetColumn(V,n,&vi));
309: PetscCall(BVGetColumn(U,n-1,&ui1));
310: PetscCall(MatMult(svd->AT,ui1,vi));
311: PetscCall(BVRestoreColumn(V,n,&vi));
312: PetscCall(BVRestoreColumn(U,n-1,&ui1));
313: PetscCall(BVOrthogonalizeColumn(V,n,NULL,&b,NULL));
314: beta[n-1] = b;
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: /*
319: Custom CGS orthogonalization, preprocess after first orthogonalization
320: */
321: static PetscErrorCode SVDOrthogonalizeCGS(BV V,PetscInt i,PetscScalar* h,PetscReal a,BVOrthogRefineType refine,PetscReal eta,PetscReal *norm)
322: {
323: PetscReal sum,onorm;
324: PetscScalar dot;
325: PetscInt j;
327: PetscFunctionBegin;
328: switch (refine) {
329: case BV_ORTHOG_REFINE_NEVER:
330: PetscCall(BVNormColumn(V,i,NORM_2,norm));
331: break;
332: case BV_ORTHOG_REFINE_ALWAYS:
333: PetscCall(BVSetActiveColumns(V,0,i));
334: PetscCall(BVDotColumn(V,i,h));
335: PetscCall(BVMultColumn(V,-1.0,1.0,i,h));
336: PetscCall(BVNormColumn(V,i,NORM_2,norm));
337: break;
338: case BV_ORTHOG_REFINE_IFNEEDED:
339: dot = h[i];
340: onorm = PetscSqrtReal(PetscRealPart(dot)) / a;
341: sum = 0.0;
342: for (j=0;j<i;j++) {
343: sum += PetscRealPart(h[j] * PetscConj(h[j]));
344: }
345: *norm = PetscRealPart(dot)/(a*a) - sum;
346: if (*norm>0.0) *norm = PetscSqrtReal(*norm);
347: else PetscCall(BVNormColumn(V,i,NORM_2,norm));
348: if (*norm < eta*onorm) {
349: PetscCall(BVSetActiveColumns(V,0,i));
350: PetscCall(BVDotColumn(V,i,h));
351: PetscCall(BVMultColumn(V,-1.0,1.0,i,h));
352: PetscCall(BVNormColumn(V,i,NORM_2,norm));
353: }
354: break;
355: }
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
359: static PetscErrorCode SVDOneSideTRLanczosCGS(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,BV U,PetscInt nconv,PetscInt l,PetscInt n,PetscScalar* work)
360: {
361: PetscReal a,b,eta;
362: PetscInt i,j,k=nconv+l;
363: Vec ui,ui1,vi;
364: BVOrthogRefineType refine;
366: PetscFunctionBegin;
367: PetscCall(BVGetColumn(V,k,&vi));
368: PetscCall(BVGetColumn(U,k,&ui));
369: PetscCall(MatMult(svd->A,vi,ui));
370: PetscCall(BVRestoreColumn(V,k,&vi));
371: PetscCall(BVRestoreColumn(U,k,&ui));
372: if (l>0) {
373: PetscCall(BVSetActiveColumns(U,nconv,n));
374: for (i=0;i<l;i++) work[i]=beta[i+nconv];
375: PetscCall(BVMultColumn(U,-1.0,1.0,k,work));
376: }
377: PetscCall(BVGetOrthogonalization(V,NULL,&refine,&eta,NULL));
379: for (i=k+1;i<n;i++) {
380: PetscCall(BVGetColumn(V,i,&vi));
381: PetscCall(BVGetColumn(U,i-1,&ui1));
382: PetscCall(MatMult(svd->AT,ui1,vi));
383: PetscCall(BVRestoreColumn(V,i,&vi));
384: PetscCall(BVRestoreColumn(U,i-1,&ui1));
385: PetscCall(BVNormColumnBegin(U,i-1,NORM_2,&a));
386: if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
387: PetscCall(BVSetActiveColumns(V,0,i+1));
388: PetscCall(BVGetColumn(V,i,&vi));
389: PetscCall(BVDotVecBegin(V,vi,work));
390: } else {
391: PetscCall(BVSetActiveColumns(V,0,i));
392: PetscCall(BVDotColumnBegin(V,i,work));
393: }
394: PetscCall(BVNormColumnEnd(U,i-1,NORM_2,&a));
395: if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
396: PetscCall(BVDotVecEnd(V,vi,work));
397: PetscCall(BVRestoreColumn(V,i,&vi));
398: PetscCall(BVSetActiveColumns(V,0,i));
399: } else PetscCall(BVDotColumnEnd(V,i,work));
401: PetscCall(BVScaleColumn(U,i-1,1.0/a));
402: for (j=0;j<i;j++) work[j] = work[j] / a;
403: PetscCall(BVMultColumn(V,-1.0,1.0/a,i,work));
404: PetscCall(SVDOrthogonalizeCGS(V,i,work,a,refine,eta,&b));
405: PetscCall(BVScaleColumn(V,i,1.0/b));
406: PetscCheck(PetscAbsReal(b)>10*PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)svd),PETSC_ERR_PLIB,"Recurrence generated a zero vector; use a two-sided variant");
408: PetscCall(BVGetColumn(V,i,&vi));
409: PetscCall(BVGetColumn(U,i,&ui));
410: PetscCall(BVGetColumn(U,i-1,&ui1));
411: PetscCall(MatMult(svd->A,vi,ui));
412: PetscCall(VecAXPY(ui,-b,ui1));
413: PetscCall(BVRestoreColumn(V,i,&vi));
414: PetscCall(BVRestoreColumn(U,i,&ui));
415: PetscCall(BVRestoreColumn(U,i-1,&ui1));
417: alpha[i-1] = a;
418: beta[i-1] = b;
419: }
421: PetscCall(BVGetColumn(V,n,&vi));
422: PetscCall(BVGetColumn(U,n-1,&ui1));
423: PetscCall(MatMult(svd->AT,ui1,vi));
424: PetscCall(BVRestoreColumn(V,n,&vi));
425: PetscCall(BVRestoreColumn(U,n-1,&ui1));
427: PetscCall(BVNormColumnBegin(svd->U,n-1,NORM_2,&a));
428: if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
429: PetscCall(BVSetActiveColumns(V,0,n+1));
430: PetscCall(BVGetColumn(V,n,&vi));
431: PetscCall(BVDotVecBegin(V,vi,work));
432: } else {
433: PetscCall(BVSetActiveColumns(V,0,n));
434: PetscCall(BVDotColumnBegin(V,n,work));
435: }
436: PetscCall(BVNormColumnEnd(svd->U,n-1,NORM_2,&a));
437: if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
438: PetscCall(BVDotVecEnd(V,vi,work));
439: PetscCall(BVRestoreColumn(V,n,&vi));
440: } else PetscCall(BVDotColumnEnd(V,n,work));
442: PetscCall(BVScaleColumn(U,n-1,1.0/a));
443: for (j=0;j<n;j++) work[j] = work[j] / a;
444: PetscCall(BVMultColumn(V,-1.0,1.0/a,n,work));
445: PetscCall(SVDOrthogonalizeCGS(V,n,work,a,refine,eta,&b));
446: PetscCall(BVSetActiveColumns(V,nconv,n));
447: alpha[n-1] = a;
448: beta[n-1] = b;
449: PetscFunctionReturn(PETSC_SUCCESS);
450: }
452: static PetscErrorCode SVDSolve_TRLanczos(SVD svd)
453: {
454: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
455: PetscReal *alpha,*beta;
456: PetscScalar *swork=NULL,*w;
457: PetscInt i,k,l,nv,ld;
458: Mat U,V;
459: PetscBool breakdown=PETSC_FALSE;
460: BVOrthogType orthog;
462: PetscFunctionBegin;
463: PetscCall(PetscCitationsRegister(citation,&cited));
464: /* allocate working space */
465: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
466: PetscCall(BVGetOrthogonalization(svd->V,&orthog,NULL,NULL,NULL));
467: PetscCall(PetscMalloc1(ld,&w));
468: if (lanczos->oneside) PetscCall(PetscMalloc1(svd->ncv+1,&swork));
470: /* normalize start vector */
471: if (!svd->nini) {
472: PetscCall(BVSetRandomColumn(svd->V,0));
473: PetscCall(BVOrthonormalizeColumn(svd->V,0,PETSC_TRUE,NULL,NULL));
474: }
476: l = 0;
477: while (svd->reason == SVD_CONVERGED_ITERATING) {
478: svd->its++;
480: /* inner loop */
481: nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
482: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
483: beta = alpha + ld;
484: if (lanczos->oneside) {
485: if (orthog == BV_ORTHOG_MGS) PetscCall(SVDOneSideTRLanczosMGS(svd,alpha,beta,svd->V,svd->U,svd->nconv,l,nv,swork));
486: else PetscCall(SVDOneSideTRLanczosCGS(svd,alpha,beta,svd->V,svd->U,svd->nconv,l,nv,swork));
487: } else PetscCall(SVDTwoSideLanczos(svd,alpha,beta,svd->V,svd->U,svd->nconv+l,&nv,&breakdown));
488: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
489: PetscCall(BVScaleColumn(svd->V,nv,1.0/beta[nv-1]));
490: PetscCall(BVSetActiveColumns(svd->V,svd->nconv,nv));
491: PetscCall(BVSetActiveColumns(svd->U,svd->nconv,nv));
493: /* solve projected problem */
494: PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l));
495: PetscCall(DSSVDSetDimensions(svd->ds,nv));
496: PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
497: PetscCall(DSSolve(svd->ds,w,NULL));
498: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
499: PetscCall(DSUpdateExtraRow(svd->ds));
500: PetscCall(DSSynchronize(svd->ds,w,NULL));
501: for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
503: /* check convergence */
504: PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,1.0,&k));
505: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
507: /* update l */
508: if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
509: else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
510: if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
511: if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
513: if (svd->reason == SVD_CONVERGED_ITERATING) {
514: if (PetscUnlikely(breakdown || k==nv)) {
515: /* Start a new bidiagonalization */
516: PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
517: if (k<svd->nsv) {
518: PetscCall(BVSetRandomColumn(svd->V,k));
519: PetscCall(BVOrthonormalizeColumn(svd->V,k,PETSC_FALSE,NULL,&breakdown));
520: if (breakdown) {
521: svd->reason = SVD_DIVERGED_BREAKDOWN;
522: PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
523: }
524: }
525: } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
526: }
528: /* compute converged singular vectors and restart vectors */
529: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&V));
530: PetscCall(BVMultInPlace(svd->V,V,svd->nconv,k+l));
531: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&V));
532: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
533: PetscCall(BVMultInPlace(svd->U,U,svd->nconv,k+l));
534: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
536: /* copy the last vector to be the next initial vector */
537: if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(svd->V,nv,k+l));
539: svd->nconv = k;
540: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv));
541: }
543: /* orthonormalize U columns in one side method */
544: if (lanczos->oneside) {
545: for (i=0;i<svd->nconv;i++) PetscCall(BVOrthonormalizeColumn(svd->U,i,PETSC_FALSE,NULL,NULL));
546: }
548: /* free working space */
549: PetscCall(PetscFree(w));
550: if (swork) PetscCall(PetscFree(swork));
551: PetscCall(DSTruncate(svd->ds,svd->nconv,PETSC_TRUE));
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
555: 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)
556: {
557: PetscInt i;
558: Vec u,v,ou=svd->workl[0];
559: PetscBool lindep=PETSC_FALSE;
560: PetscReal norm;
562: PetscFunctionBegin;
563: for (i=k;i<*n;i++) {
564: PetscCall(BVGetColumn(V,i,&v));
565: PetscCall(BVGetColumn(U,i,&u));
566: PetscCall(MatMult(A,v,u));
567: PetscCall(BVRestoreColumn(V,i,&v));
568: PetscCall(BVRestoreColumn(U,i,&u));
569: PetscCall(BVOrthonormalizeColumn(U,i,PETSC_FALSE,alpha+i,&lindep));
570: omega[i] = PetscSign(alpha[i]);
571: if (PetscUnlikely(lindep)) {
572: *n = i;
573: break;
574: }
576: PetscCall(BVGetColumn(V,i+1,&v));
577: PetscCall(BVGetColumn(U,i,&u));
578: PetscCall(VecPointwiseMult(ou,svd->omega,u));
579: PetscCall(MatMult(AT,ou,v));
580: PetscCall(BVRestoreColumn(V,i+1,&v));
581: PetscCall(BVRestoreColumn(U,i,&u));
582: PetscCall(BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,&norm,&lindep));
583: beta[i] = omega[i]*norm;
584: if (PetscUnlikely(lindep)) {
585: *n = i+1;
586: break;
587: }
588: }
590: if (breakdown) *breakdown = lindep;
591: PetscFunctionReturn(PETSC_SUCCESS);
592: }
594: static PetscErrorCode SVDSolve_TRLanczos_HSVD(SVD svd)
595: {
596: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
597: PetscReal *alpha,*beta,*omega;
598: PetscScalar *w;
599: PetscInt i,k,l,nv,ld,nini;
600: Mat UU,VV,D,A,AT;
601: BV U,V;
602: PetscBool breakdown=PETSC_FALSE;
603: BVOrthogType orthog;
604: Vec vomega;
606: PetscFunctionBegin;
607: /* undo the effect of swapping in this function */
608: if (svd->swapped) {
609: A = svd->AT;
610: AT = svd->A;
611: U = svd->V;
612: V = svd->U;
613: nini = svd->ninil;
614: } else {
615: A = svd->A;
616: AT = svd->AT;
617: U = svd->U;
618: V = svd->V;
619: nini = svd->nini;
620: }
621: /* allocate working space */
622: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
623: PetscCall(BVGetOrthogonalization(V,&orthog,NULL,NULL,NULL));
624: PetscCall(PetscMalloc1(ld,&w));
625: PetscCheck(!lanczos->oneside,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Oneside orthogonalization not supported for HSVD");
627: /* normalize start vector */
628: if (!nini) {
629: PetscCall(BVSetRandomColumn(V,0));
630: PetscCall(BVOrthonormalizeColumn(V,0,PETSC_TRUE,NULL,NULL));
631: }
633: l = 0;
634: while (svd->reason == SVD_CONVERGED_ITERATING) {
635: svd->its++;
637: /* inner loop */
638: nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
639: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
640: beta = alpha + ld;
641: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&omega));
642: PetscCall(SVDLanczosHSVD(svd,alpha,beta,omega,A,AT,V,U,svd->nconv+l,&nv,&breakdown));
643: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
644: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&omega));
645: PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
646: PetscCall(BVSetActiveColumns(U,svd->nconv,nv));
648: /* solve projected problem */
649: PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l));
650: PetscCall(DSHSVDSetDimensions(svd->ds,nv));
651: PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
652: PetscCall(DSSolve(svd->ds,w,NULL));
653: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
654: PetscCall(DSUpdateExtraRow(svd->ds));
655: PetscCall(DSSynchronize(svd->ds,w,NULL));
656: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&omega));
657: for (i=svd->nconv;i<nv;i++) {
658: svd->sigma[i] = PetscRealPart(w[i]);
659: svd->sign[i] = omega[i];
660: }
661: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&omega));
663: /* check convergence */
664: PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,1.0,&k));
665: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
667: /* update l */
668: if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
669: else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
670: if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
671: if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
673: if (svd->reason == SVD_CONVERGED_ITERATING) {
674: if (PetscUnlikely(breakdown || k==nv)) {
675: /* Start a new bidiagonalization */
676: PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
677: if (k<svd->nsv) {
678: PetscCall(BVSetRandomColumn(V,k));
679: PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,NULL,&breakdown));
680: if (breakdown) {
681: svd->reason = SVD_DIVERGED_BREAKDOWN;
682: PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
683: }
684: }
685: } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
686: }
688: /* compute converged singular vectors and restart vectors */
689: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&VV));
690: PetscCall(BVMultInPlace(V,VV,svd->nconv,k+l));
691: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&VV));
692: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&UU));
693: PetscCall(BVMultInPlace(U,UU,svd->nconv,k+l));
694: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&UU));
696: /* copy the last vector of V to be the next initial vector
697: and change signature matrix of U */
698: if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) {
699: PetscCall(BVCopyColumn(V,nv,k+l));
700: PetscCall(BVSetActiveColumns(U,0,k+l));
701: PetscCall(DSGetMatAndColumn(svd->ds,DS_MAT_D,0,&D,&vomega));
702: PetscCall(BVSetSignature(U,vomega));
703: PetscCall(DSRestoreMatAndColumn(svd->ds,DS_MAT_D,0,&D,&vomega));
704: }
706: svd->nconv = k;
707: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv));
708: }
710: /* free working space */
711: PetscCall(PetscFree(w));
712: PetscCall(DSTruncate(svd->ds,svd->nconv,PETSC_TRUE));
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }
716: /* Given n computed generalized singular values in sigmain, backtransform them
717: in sigmaout by undoing scaling and reciprocating if swapped=true. Also updates vectors V
718: if given. If sigmaout=NULL then the result overwrites sigmain. */
719: static PetscErrorCode SVDLanczosBackTransform(SVD svd,PetscInt n,PetscReal *sigmain,PetscReal *sigmaout,BV V)
720: {
721: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
722: PetscInt i;
723: PetscReal c,s,r,f,scalef;
725: PetscFunctionBegin;
726: scalef = svd->swapped? 1.0/lanczos->scalef: lanczos->scalef;
727: for (i=0;i<n;i++) {
728: if (V && scalef != 1.0) {
729: s = 1.0/PetscSqrtReal(1.0+sigmain[i]*sigmain[i]);
730: c = sigmain[i]*s;
731: r = s/scalef;
732: f = 1.0/PetscSqrtReal(c*c+r*r);
733: PetscCall(BVScaleColumn(V,i,f));
734: }
735: if (sigmaout) {
736: if (svd->swapped) sigmaout[i] = 1.0/(sigmain[i]*scalef);
737: else sigmaout[i] = sigmain[i]*scalef;
738: } else {
739: sigmain[i] *= scalef;
740: if (svd->swapped) sigmain[i] = 1.0/sigmain[i];
741: }
742: }
743: PetscFunctionReturn(PETSC_SUCCESS);
744: }
746: static PetscErrorCode SVDLanczosGSingle(SVD svd,PetscReal *alpha,PetscReal *beta,Mat Z,BV V,BV U,KSP ksp,PetscInt k,PetscInt *n,PetscBool *breakdown)
747: {
748: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
749: PetscInt i,j,m;
750: const PetscScalar *carr;
751: PetscScalar *arr;
752: Vec u,v,ut=svd->workl[0],x=svd->workr[0],v1,u1,u2;
753: PetscBool lindep=PETSC_FALSE;
755: PetscFunctionBegin;
756: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&v1));
757: PetscCall(BVGetColumn(V,k,&v));
758: PetscCall(BVGetColumn(U,k,&u));
760: /* Form ut=[u;0] */
761: PetscCall(VecZeroEntries(ut));
762: PetscCall(VecGetLocalSize(u,&m));
763: PetscCall(VecGetArrayRead(u,&carr));
764: PetscCall(VecGetArray(ut,&arr));
765: for (j=0; j<m; j++) arr[j] = carr[j];
766: PetscCall(VecRestoreArrayRead(u,&carr));
767: PetscCall(VecRestoreArray(ut,&arr));
769: /* Solve least squares problem */
770: PetscCall(KSPSolve(ksp,ut,x));
772: PetscCall(MatMult(Z,x,v));
774: PetscCall(BVRestoreColumn(U,k,&u));
775: PetscCall(BVRestoreColumn(V,k,&v));
776: PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,alpha+k,&lindep));
777: if (PetscUnlikely(lindep)) {
778: *n = k;
779: if (breakdown) *breakdown = lindep;
780: PetscFunctionReturn(PETSC_SUCCESS);
781: }
783: for (i=k+1; i<*n; i++) {
785: /* Compute vector i of BV U */
786: PetscCall(BVGetColumn(V,i-1,&v));
787: PetscCall(VecGetArray(v,&arr));
788: PetscCall(VecPlaceArray(v1,arr));
789: PetscCall(VecRestoreArray(v,&arr));
790: PetscCall(BVRestoreColumn(V,i-1,&v));
791: PetscCall(BVInsertVec(U,i,v1));
792: PetscCall(VecResetArray(v1));
793: PetscCall(BVOrthonormalizeColumn(U,i,PETSC_FALSE,beta+i-1,&lindep));
794: if (PetscUnlikely(lindep)) {
795: *n = i;
796: break;
797: }
799: /* Compute vector i of BV V */
801: PetscCall(BVGetColumn(V,i,&v));
802: PetscCall(BVGetColumn(U,i,&u));
804: /* Form ut=[u;0] */
805: PetscCall(VecGetArrayRead(u,&carr));
806: PetscCall(VecGetArray(ut,&arr));
807: for (j=0; j<m; j++) arr[j] = carr[j];
808: PetscCall(VecRestoreArrayRead(u,&carr));
809: PetscCall(VecRestoreArray(ut,&arr));
811: /* Solve least squares problem */
812: PetscCall(KSPSolve(ksp,ut,x));
814: PetscCall(MatMult(Z,x,v));
816: PetscCall(BVRestoreColumn(U,i,&u));
817: PetscCall(BVRestoreColumn(V,i,&v));
818: if (!lanczos->oneside || i==k+1) PetscCall(BVOrthonormalizeColumn(V,i,PETSC_FALSE,alpha+i,&lindep));
819: else { /* cheap computation of V[i], if restart (i==k+1) do a full reorthogonalization */
820: PetscCall(BVGetColumn(V,i,&u2));
821: PetscCall(BVGetColumn(V,i-1,&u1));
822: PetscCall(VecAXPY(u2,-beta[i-1],u1));
823: PetscCall(BVRestoreColumn(V,i-1,&u1));
824: PetscCall(VecNorm(u2,NORM_2,&alpha[i]));
825: if (alpha[i]==0.0) lindep = PETSC_TRUE;
826: else PetscCall(VecScale(u2,1.0/alpha[i]));
827: PetscCall(BVRestoreColumn(V,i,&u2));
828: }
829: if (PetscUnlikely(lindep)) {
830: *n = i;
831: break;
832: }
833: }
835: /* Compute vector n of BV U */
836: if (!lindep) {
837: PetscCall(BVGetColumn(V,*n-1,&v));
838: PetscCall(VecGetArray(v,&arr));
839: PetscCall(VecPlaceArray(v1,arr));
840: PetscCall(VecRestoreArray(v,&arr));
841: PetscCall(BVRestoreColumn(V,*n-1,&v));
842: PetscCall(BVInsertVec(U,*n,v1));
843: PetscCall(VecResetArray(v1));
844: PetscCall(BVOrthonormalizeColumn(U,*n,PETSC_FALSE,beta+*n-1,&lindep));
845: }
846: if (breakdown) *breakdown = lindep;
847: PetscCall(VecDestroy(&v1));
848: PetscFunctionReturn(PETSC_SUCCESS);
849: }
851: /* solve generalized problem with single bidiagonalization of Q_A */
852: static PetscErrorCode SVDSolve_TRLanczosGSingle(SVD svd,BV U1,BV V)
853: {
854: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
855: PetscReal *alpha,*beta,normr,scaleth,sigma0,*sigma;
856: PetscScalar *w;
857: PetscInt i,k,l,nv,ld;
858: Mat U,VV;
859: PetscBool breakdown=PETSC_FALSE;
861: PetscFunctionBegin;
862: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
863: PetscCall(PetscMalloc2(ld,&w,ld,&sigma));
864: normr = (svd->conv==SVD_CONV_ABS)? PetscMax(svd->nrma,svd->nrmb*lanczos->scalef): 1.0;
865: /* Convert scale threshold th=c/s to the corresponding c */
866: scaleth = (lanczos->scaleth!=0)? lanczos->scaleth/PetscSqrtReal(lanczos->scaleth*lanczos->scaleth+1): 0.0;
868: /* normalize start vector */
869: if (!svd->ninil) {
870: PetscCall(BVSetRandomColumn(U1,0));
871: PetscCall(BVOrthonormalizeColumn(U1,0,PETSC_TRUE,NULL,NULL));
872: }
874: l = 0;
875: while (svd->reason == SVD_CONVERGED_ITERATING) {
876: svd->its++;
878: /* inner loop */
879: nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
880: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
881: beta = alpha + ld;
882: PetscCall(SVDLanczosGSingle(svd,alpha,beta,lanczos->Z,V,U1,lanczos->ksp,svd->nconv+l,&nv,&breakdown));
883: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
884: PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
885: PetscCall(BVSetActiveColumns(U1,svd->nconv,nv));
887: /* solve projected problem */
888: PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l));
889: PetscCall(DSSVDSetDimensions(svd->ds,nv));
890: PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
891: PetscCall(DSSolve(svd->ds,w,NULL));
892: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
893: PetscCall(DSUpdateExtraRow(svd->ds));
894: PetscCall(DSSynchronize(svd->ds,w,NULL));
895: for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
897: /* check convergence */
898: PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,normr,&k));
899: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
901: sigma0 = svd->which==SVD_LARGEST? svd->sigma[0] : 1.0/svd->sigma[0];
902: if (scaleth!=0 && k==0 && sigma0>scaleth) {
904: /* Scale and start from scratch */
905: lanczos->scalef *= svd->sigma[0]/PetscSqrtReal(1-svd->sigma[0]*svd->sigma[0]);
906: PetscCall(PetscInfo(svd,"Scaling by factor %g and starting from scratch\n",(double)lanczos->scalef));
907: PetscCall(MatZUpdateScale(svd));
908: if (svd->conv==SVD_CONV_ABS) normr = PetscMax(svd->nrma,svd->nrmb*lanczos->scalef);
909: l = 0;
911: } else {
913: /* update l */
914: if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
915: else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
916: if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
917: if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
919: if (svd->reason == SVD_CONVERGED_ITERATING) {
920: if (PetscUnlikely(breakdown || k==nv)) {
921: /* Start a new bidiagonalization */
922: PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
923: if (k<svd->nsv) {
924: PetscCall(BVSetRandomColumn(U1,k));
925: PetscCall(BVOrthonormalizeColumn(U1,k,PETSC_FALSE,NULL,&breakdown));
926: if (breakdown) {
927: svd->reason = SVD_DIVERGED_BREAKDOWN;
928: PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
929: }
930: }
931: } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
932: }
934: /* compute converged singular vectors and restart vectors */
935: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
936: PetscCall(BVMultInPlace(V,U,svd->nconv,k+l));
937: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
938: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&VV));
939: PetscCall(BVMultInPlace(U1,VV,svd->nconv,k+l));
940: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&VV));
942: /* copy the last vector to be the next initial vector */
943: if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(U1,nv,k+l));
944: }
946: svd->nconv = k;
947: PetscCall(SVDLanczosBackTransform(svd,nv,svd->sigma,sigma,NULL));
948: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,sigma,svd->errest,nv));
949: }
951: PetscCall(PetscFree2(w,sigma));
952: PetscFunctionReturn(PETSC_SUCCESS);
953: }
955: /* Move generalized left singular vectors (0..nconv) from U1 and U2 to its final destination svd->U (single variant) */
956: static inline PetscErrorCode SVDLeftSingularVectors_Single(SVD svd,BV U1,BV U2)
957: {
958: PetscInt i,k,m,p;
959: Vec u,u1,u2;
960: PetscScalar *ua,*u2a;
961: const PetscScalar *u1a;
962: PetscReal s;
964: PetscFunctionBegin;
965: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
966: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
967: for (i=0;i<svd->nconv;i++) {
968: PetscCall(BVGetColumn(U1,i,&u1));
969: PetscCall(BVGetColumn(U2,i,&u2));
970: PetscCall(BVGetColumn(svd->U,i,&u));
971: PetscCall(VecGetArrayRead(u1,&u1a));
972: PetscCall(VecGetArray(u,&ua));
973: PetscCall(VecGetArray(u2,&u2a));
974: /* Copy column from U1 to upper part of u */
975: for (k=0;k<m;k++) ua[k] = u1a[k];
976: /* Copy column from lower part of U to U2. Orthogonalize column in U2 and copy back to U */
977: for (k=0;k<p;k++) u2a[k] = ua[m+k];
978: PetscCall(VecRestoreArray(u2,&u2a));
979: PetscCall(BVRestoreColumn(U2,i,&u2));
980: PetscCall(BVOrthonormalizeColumn(U2,i,PETSC_FALSE,&s,NULL));
981: PetscCall(BVGetColumn(U2,i,&u2));
982: PetscCall(VecGetArray(u2,&u2a));
983: for (k=0;k<p;k++) ua[m+k] = u2a[k];
984: /* Update singular value */
985: svd->sigma[i] /= s;
986: PetscCall(VecRestoreArrayRead(u1,&u1a));
987: PetscCall(VecRestoreArray(u,&ua));
988: PetscCall(VecRestoreArray(u2,&u2a));
989: PetscCall(BVRestoreColumn(U1,i,&u1));
990: PetscCall(BVRestoreColumn(U2,i,&u2));
991: PetscCall(BVRestoreColumn(svd->U,i,&u));
992: }
993: PetscFunctionReturn(PETSC_SUCCESS);
994: }
996: 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)
997: {
998: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
999: PetscInt i,j,m,p;
1000: const PetscScalar *carr;
1001: PetscScalar *arr,*u2arr;
1002: Vec u,v,ut=svd->workl[0],x=svd->workr[0],v1,u1,u2;
1003: PetscBool lindep=PETSC_FALSE,lindep1=PETSC_FALSE,lindep2=PETSC_FALSE;
1005: PetscFunctionBegin;
1006: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&v1));
1007: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1008: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1010: for (i=k; i<*n; i++) {
1011: /* Compute vector i of BV U1 */
1012: PetscCall(BVGetColumn(V,i,&v));
1013: PetscCall(VecGetArrayRead(v,&carr));
1014: PetscCall(VecPlaceArray(v1,carr));
1015: PetscCall(BVInsertVec(U1,i,v1));
1016: PetscCall(VecResetArray(v1));
1017: if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(U1,i,PETSC_FALSE,alpha+i,&lindep1));
1018: else { /* cheap computation of U1[i], if restart (i==k) do a full reorthogonalization */
1019: PetscCall(BVGetColumn(U1,i,&u2));
1020: if (i>0) {
1021: PetscCall(BVGetColumn(U1,i-1,&u1));
1022: PetscCall(VecAXPY(u2,-beta[i-1],u1));
1023: PetscCall(BVRestoreColumn(U1,i-1,&u1));
1024: }
1025: PetscCall(VecNorm(u2,NORM_2,&alpha[i]));
1026: if (alpha[i]==0.0) lindep = PETSC_TRUE;
1027: else PetscCall(VecScale(u2,1.0/alpha[i]));
1028: PetscCall(BVRestoreColumn(U1,i,&u2));
1029: }
1031: /* Compute vector i of BV U2 */
1032: PetscCall(BVGetColumn(U2,i,&u2));
1033: PetscCall(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: PetscCall(VecRestoreArray(u2,&u2arr));
1040: PetscCall(VecRestoreArrayRead(v,&carr));
1041: PetscCall(BVRestoreColumn(V,i,&v));
1042: if (lanczos->oneside && i>k) { /* cheap computation of U2[i], if restart (i==k) do a full reorthogonalization */
1043: if (i>0) {
1044: PetscCall(BVGetColumn(U2,i-1,&u1));
1045: PetscCall(VecAXPY(u2,(i%2)?betah[i-1]:-betah[i-1],u1));
1046: PetscCall(BVRestoreColumn(U2,i-1,&u1));
1047: }
1048: PetscCall(VecNorm(u2,NORM_2,&alphah[i]));
1049: if (alphah[i]==0.0) lindep = PETSC_TRUE;
1050: else PetscCall(VecScale(u2,1.0/alphah[i]));
1051: }
1052: PetscCall(BVRestoreColumn(U2,i,&u2));
1053: if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(U2,i,PETSC_FALSE,alphah+i,&lindep2));
1054: if (i%2) alphah[i] = -alphah[i];
1055: if (PetscUnlikely(lindep1 || lindep2)) {
1056: lindep = PETSC_TRUE;
1057: *n = i;
1058: break;
1059: }
1061: /* Compute vector i+1 of BV V */
1062: PetscCall(BVGetColumn(V,i+1,&v));
1063: /* Form ut=[u;0] */
1064: PetscCall(BVGetColumn(U1,i,&u));
1065: PetscCall(VecZeroEntries(ut));
1066: PetscCall(VecGetArrayRead(u,&carr));
1067: PetscCall(VecGetArray(ut,&arr));
1068: for (j=0; j<m; j++) arr[j] = carr[j];
1069: PetscCall(VecRestoreArrayRead(u,&carr));
1070: PetscCall(VecRestoreArray(ut,&arr));
1071: /* Solve least squares problem */
1072: PetscCall(KSPSolve(ksp,ut,x));
1073: PetscCall(MatMult(Z,x,v));
1074: PetscCall(BVRestoreColumn(U1,i,&u));
1075: PetscCall(BVRestoreColumn(V,i+1,&v));
1076: PetscCall(BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,beta+i,&lindep));
1077: betah[i] = -alpha[i]*beta[i]/alphah[i];
1078: if (PetscUnlikely(lindep)) {
1079: *n = i;
1080: break;
1081: }
1082: }
1083: if (breakdown) *breakdown = lindep;
1084: PetscCall(VecDestroy(&v1));
1085: PetscFunctionReturn(PETSC_SUCCESS);
1086: }
1088: /* generate random initial vector in column k for joint upper-upper bidiagonalization */
1089: static inline PetscErrorCode SVDInitialVectorGUpper(SVD svd,BV V,BV U1,PetscInt k,PetscBool *breakdown)
1090: {
1091: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1092: Vec u,v,ut=svd->workl[0],x=svd->workr[0];
1093: PetscInt m,j;
1094: PetscScalar *arr;
1095: const PetscScalar *carr;
1097: PetscFunctionBegin;
1098: /* Form ut=[u;0] where u is the k-th column of U1 */
1099: PetscCall(VecZeroEntries(ut));
1100: PetscCall(BVGetColumn(U1,k,&u));
1101: PetscCall(VecGetLocalSize(u,&m));
1102: PetscCall(VecGetArrayRead(u,&carr));
1103: PetscCall(VecGetArray(ut,&arr));
1104: for (j=0; j<m; j++) arr[j] = carr[j];
1105: PetscCall(VecRestoreArrayRead(u,&carr));
1106: PetscCall(VecRestoreArray(ut,&arr));
1107: PetscCall(BVRestoreColumn(U1,k,&u));
1108: /* Solve least squares problem Z*x=ut for x. Then set v=Z*x */
1109: PetscCall(KSPSolve(lanczos->ksp,ut,x));
1110: PetscCall(BVGetColumn(V,k,&v));
1111: PetscCall(MatMult(lanczos->Z,x,v));
1112: PetscCall(BVRestoreColumn(V,k,&v));
1113: if (breakdown) PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,NULL,breakdown));
1114: else PetscCall(BVOrthonormalizeColumn(V,k,PETSC_TRUE,NULL,NULL));
1115: PetscFunctionReturn(PETSC_SUCCESS);
1116: }
1118: /* solve generalized problem with joint upper-upper bidiagonalization */
1119: static PetscErrorCode SVDSolve_TRLanczosGUpper(SVD svd,BV U1,BV U2,BV V)
1120: {
1121: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1122: PetscReal *alpha,*beta,*alphah,*betah,normr,sigma0,*sigma;
1123: PetscScalar *w;
1124: PetscInt i,k,l,nv,ld;
1125: Mat U,Vmat,X;
1126: PetscBool breakdown=PETSC_FALSE;
1128: PetscFunctionBegin;
1129: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
1130: PetscCall(PetscMalloc2(ld,&w,ld,&sigma));
1131: normr = (svd->conv==SVD_CONV_ABS)? PetscMax(svd->nrma,svd->nrmb*lanczos->scalef): 1.0;
1133: /* normalize start vector */
1134: if (!svd->ninil) PetscCall(BVSetRandomColumn(U1,0));
1135: PetscCall(SVDInitialVectorGUpper(svd,V,U1,0,NULL));
1137: l = 0;
1138: while (svd->reason == SVD_CONVERGED_ITERATING) {
1139: svd->its++;
1141: /* inner loop */
1142: nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
1143: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
1144: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&alphah));
1145: beta = alpha + ld;
1146: betah = alpha + 2*ld;
1147: PetscCall(SVDLanczosGUpper(svd,alpha,beta,alphah,betah,lanczos->Z,U1,U2,V,lanczos->ksp,svd->nconv+l,&nv,&breakdown));
1148: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
1149: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&alphah));
1150: PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
1151: PetscCall(BVSetActiveColumns(U1,svd->nconv,nv));
1152: PetscCall(BVSetActiveColumns(U2,svd->nconv,nv));
1154: /* solve projected problem */
1155: PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l));
1156: PetscCall(DSGSVDSetDimensions(svd->ds,nv,nv));
1157: PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
1158: PetscCall(DSSolve(svd->ds,w,NULL));
1159: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
1160: PetscCall(DSUpdateExtraRow(svd->ds));
1161: PetscCall(DSSynchronize(svd->ds,w,NULL));
1162: for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
1164: /* check convergence */
1165: PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,normr,&k));
1166: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
1168: sigma0 = svd->which==SVD_LARGEST? svd->sigma[0] : 1.0/svd->sigma[0];
1169: if (lanczos->scaleth!=0 && k==0 && sigma0>lanczos->scaleth) {
1171: /* Scale and start from scratch */
1172: lanczos->scalef *= svd->sigma[0];
1173: PetscCall(PetscInfo(svd,"Scaling by factor %g and starting from scratch\n",(double)lanczos->scalef));
1174: PetscCall(MatZUpdateScale(svd));
1175: if (svd->conv==SVD_CONV_ABS) normr = PetscMax(svd->nrma,svd->nrmb*lanczos->scalef);
1176: l = 0;
1177: if (!svd->ninil) PetscCall(BVSetRandomColumn(U1,0));
1178: PetscCall(SVDInitialVectorGUpper(svd,V,U1,0,NULL));
1180: } else {
1182: /* update l */
1183: if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
1184: else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
1185: if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
1186: if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
1188: if (svd->reason == SVD_CONVERGED_ITERATING) {
1189: if (PetscUnlikely(breakdown || k==nv)) {
1190: /* Start a new bidiagonalization */
1191: PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
1192: if (k<svd->nsv) {
1193: PetscCall(BVSetRandomColumn(U1,k));
1194: PetscCall(SVDInitialVectorGUpper(svd,V,U1,k,&breakdown));
1195: if (breakdown) {
1196: svd->reason = SVD_DIVERGED_BREAKDOWN;
1197: PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
1198: }
1199: }
1200: } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
1201: }
1202: /* compute converged singular vectors and restart vectors */
1203: PetscCall(DSGetMat(svd->ds,DS_MAT_X,&X));
1204: PetscCall(BVMultInPlace(V,X,svd->nconv,k+l));
1205: PetscCall(DSRestoreMat(svd->ds,DS_MAT_X,&X));
1206: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
1207: PetscCall(BVMultInPlace(U1,U,svd->nconv,k+l));
1208: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
1209: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&Vmat));
1210: PetscCall(BVMultInPlace(U2,Vmat,svd->nconv,k+l));
1211: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&Vmat));
1213: /* copy the last vector to be the next initial vector */
1214: if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(V,nv,k+l));
1215: }
1217: svd->nconv = k;
1218: PetscCall(SVDLanczosBackTransform(svd,nv,svd->sigma,sigma,NULL));
1219: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,sigma,svd->errest,nv));
1220: }
1222: PetscCall(PetscFree2(w,sigma));
1223: PetscFunctionReturn(PETSC_SUCCESS);
1224: }
1226: /* Move generalized left singular vectors (0..nconv) from U1 and U2 to its final destination svd->U (upper and lower variants) */
1227: static inline PetscErrorCode SVDLeftSingularVectors(SVD svd,BV U1,BV U2)
1228: {
1229: PetscInt i,k,m,p;
1230: Vec u,u1,u2;
1231: PetscScalar *ua;
1232: const PetscScalar *u1a,*u2a;
1234: PetscFunctionBegin;
1235: PetscCall(BVGetSizes(U1,&m,NULL,NULL));
1236: PetscCall(BVGetSizes(U2,&p,NULL,NULL));
1237: for (i=0;i<svd->nconv;i++) {
1238: PetscCall(BVGetColumn(U1,i,&u1));
1239: PetscCall(BVGetColumn(U2,i,&u2));
1240: PetscCall(BVGetColumn(svd->U,i,&u));
1241: PetscCall(VecGetArrayRead(u1,&u1a));
1242: PetscCall(VecGetArrayRead(u2,&u2a));
1243: PetscCall(VecGetArray(u,&ua));
1244: /* Copy column from u1 to upper part of u */
1245: for (k=0;k<m;k++) ua[k] = u1a[k];
1246: /* Copy column from u2 to lower part of u */
1247: for (k=0;k<p;k++) ua[m+k] = u2a[k];
1248: PetscCall(VecRestoreArrayRead(u1,&u1a));
1249: PetscCall(VecRestoreArrayRead(u2,&u2a));
1250: PetscCall(VecRestoreArray(u,&ua));
1251: PetscCall(BVRestoreColumn(U1,i,&u1));
1252: PetscCall(BVRestoreColumn(U2,i,&u2));
1253: PetscCall(BVRestoreColumn(svd->U,i,&u));
1254: }
1255: PetscFunctionReturn(PETSC_SUCCESS);
1256: }
1258: 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)
1259: {
1260: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1261: PetscInt i,j,m,p;
1262: const PetscScalar *carr;
1263: PetscScalar *arr,*u2arr;
1264: Vec u,v,ut=svd->workl[0],x=svd->workr[0],v1,u1,u2;
1265: PetscBool lindep=PETSC_FALSE;
1267: PetscFunctionBegin;
1268: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&v1));
1269: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1270: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1272: for (i=k; i<*n; i++) {
1273: /* Compute vector i of BV U2 */
1274: PetscCall(BVGetColumn(V,i,&v));
1275: PetscCall(VecGetArrayRead(v,&carr));
1276: PetscCall(BVGetColumn(U2,i,&u2));
1277: PetscCall(VecGetArray(u2,&u2arr));
1278: if (i%2) {
1279: for (j=0; j<p; j++) u2arr[j] = -carr[m+j];
1280: } else {
1281: for (j=0; j<p; j++) u2arr[j] = carr[m+j];
1282: }
1283: PetscCall(VecRestoreArray(u2,&u2arr));
1284: if (lanczos->oneside && i>k) { /* cheap computation of U2[i], if restart (i==k) do a full reorthogonalization */
1285: if (i>0) {
1286: PetscCall(BVGetColumn(U2,i-1,&u1));
1287: PetscCall(VecAXPY(u2,(i%2)?betah[i-1]:-betah[i-1],u1));
1288: PetscCall(BVRestoreColumn(U2,i-1,&u1));
1289: }
1290: PetscCall(VecNorm(u2,NORM_2,&alphah[i]));
1291: if (alphah[i]==0.0) lindep = PETSC_TRUE;
1292: else PetscCall(VecScale(u2,1.0/alphah[i]));
1293: }
1294: PetscCall(BVRestoreColumn(U2,i,&u2));
1295: if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(U2,i,PETSC_FALSE,alphah+i,&lindep));
1296: if (i%2) alphah[i] = -alphah[i];
1297: if (PetscUnlikely(lindep)) {
1298: PetscCall(BVRestoreColumn(V,i,&v));
1299: *n = i;
1300: break;
1301: }
1303: /* Compute vector i+1 of BV U1 */
1304: PetscCall(VecPlaceArray(v1,carr));
1305: PetscCall(BVInsertVec(U1,i+1,v1));
1306: PetscCall(VecResetArray(v1));
1307: PetscCall(BVOrthonormalizeColumn(U1,i+1,PETSC_FALSE,beta+i,&lindep));
1308: PetscCall(VecRestoreArrayRead(v,&carr));
1309: PetscCall(BVRestoreColumn(V,i,&v));
1310: if (PetscUnlikely(lindep)) {
1311: *n = i+1;
1312: break;
1313: }
1315: /* Compute vector i+1 of BV V */
1316: PetscCall(BVGetColumn(V,i+1,&v));
1317: /* Form ut=[u;0] where u is column i+1 of BV U1 */
1318: PetscCall(BVGetColumn(U1,i+1,&u));
1319: PetscCall(VecZeroEntries(ut));
1320: PetscCall(VecGetArrayRead(u,&carr));
1321: PetscCall(VecGetArray(ut,&arr));
1322: for (j=0; j<m; j++) arr[j] = carr[j];
1323: PetscCall(VecRestoreArrayRead(u,&carr));
1324: PetscCall(VecRestoreArray(ut,&arr));
1325: /* Solve least squares problem */
1326: PetscCall(KSPSolve(ksp,ut,x));
1327: PetscCall(MatMult(Z,x,v));
1328: PetscCall(BVRestoreColumn(U1,i+1,&u));
1329: PetscCall(BVRestoreColumn(V,i+1,&v));
1330: if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,alpha+i+1,&lindep));
1331: else { /* cheap computation of V[i+1], if restart (i==k) do a full reorthogonalization */
1332: PetscCall(BVGetColumn(V,i+1,&u2));
1333: PetscCall(BVGetColumn(V,i,&u1));
1334: PetscCall(VecAXPY(u2,-beta[i],u1));
1335: PetscCall(BVRestoreColumn(V,i,&u1));
1336: PetscCall(VecNorm(u2,NORM_2,&alpha[i+1]));
1337: if (alpha[i+1]==0.0) lindep = PETSC_TRUE;
1338: else PetscCall(VecScale(u2,1.0/alpha[i+1]));
1339: PetscCall(BVRestoreColumn(V,i+1,&u2));
1340: }
1341: betah[i] = -alpha[i+1]*beta[i]/alphah[i];
1342: if (PetscUnlikely(lindep)) {
1343: *n = i+1;
1344: break;
1345: }
1346: }
1347: if (breakdown) *breakdown = lindep;
1348: PetscCall(VecDestroy(&v1));
1349: PetscFunctionReturn(PETSC_SUCCESS);
1350: }
1352: /* generate random initial vector in column k for joint lower-upper bidiagonalization */
1353: static inline PetscErrorCode SVDInitialVectorGLower(SVD svd,BV V,BV U1,BV U2,PetscInt k,PetscBool *breakdown)
1354: {
1355: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1356: const PetscScalar *carr;
1357: PetscScalar *arr;
1358: PetscReal *alpha;
1359: PetscInt j,m,p;
1360: Vec u,uh,v,ut=svd->workl[0],x=svd->workr[0];
1362: PetscFunctionBegin;
1363: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1364: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1365: /* Form ut=[0;uh], where uh is the k-th column of U2 */
1366: PetscCall(BVGetColumn(U2,k,&uh));
1367: PetscCall(VecZeroEntries(ut));
1368: PetscCall(VecGetArrayRead(uh,&carr));
1369: PetscCall(VecGetArray(ut,&arr));
1370: for (j=0; j<p; j++) arr[m+j] = carr[j];
1371: PetscCall(VecRestoreArrayRead(uh,&carr));
1372: PetscCall(VecRestoreArray(ut,&arr));
1373: PetscCall(BVRestoreColumn(U2,k,&uh));
1374: /* Solve least squares problem Z*x=ut for x. Then set ut=Z*x */
1375: PetscCall(KSPSolve(lanczos->ksp,ut,x));
1376: PetscCall(MatMult(lanczos->Z,x,ut));
1377: /* Form u, column k of BV U1, as the upper part of ut and orthonormalize */
1378: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&u));
1379: PetscCall(VecGetArrayRead(ut,&carr));
1380: PetscCall(VecPlaceArray(u,carr));
1381: PetscCall(BVInsertVec(U1,k,u));
1382: PetscCall(VecResetArray(u));
1383: PetscCall(VecRestoreArrayRead(ut,&carr));
1384: PetscCall(VecDestroy(&u));
1385: if (breakdown) PetscCall(BVOrthonormalizeColumn(U1,k,PETSC_FALSE,NULL,breakdown));
1386: else PetscCall(BVOrthonormalizeColumn(U1,k,PETSC_TRUE,NULL,NULL));
1388: if (!breakdown || !*breakdown) {
1389: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1390: /* Compute k-th vector of BV V */
1391: PetscCall(BVGetColumn(V,k,&v));
1392: /* Form ut=[u;0] where u is the 1st column of U1 */
1393: PetscCall(BVGetColumn(U1,k,&u));
1394: PetscCall(VecZeroEntries(ut));
1395: PetscCall(VecGetArrayRead(u,&carr));
1396: PetscCall(VecGetArray(ut,&arr));
1397: for (j=0; j<m; j++) arr[j] = carr[j];
1398: PetscCall(VecRestoreArrayRead(u,&carr));
1399: PetscCall(VecRestoreArray(ut,&arr));
1400: /* Solve least squares problem */
1401: PetscCall(KSPSolve(lanczos->ksp,ut,x));
1402: PetscCall(MatMult(lanczos->Z,x,v));
1403: PetscCall(BVRestoreColumn(U1,k,&u));
1404: PetscCall(BVRestoreColumn(V,k,&v));
1405: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
1406: if (breakdown) PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,alpha+k,breakdown));
1407: else PetscCall(BVOrthonormalizeColumn(V,k,PETSC_TRUE,alpha+k,NULL));
1408: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
1409: }
1410: PetscFunctionReturn(PETSC_SUCCESS);
1411: }
1413: /* solve generalized problem with joint lower-upper bidiagonalization */
1414: static PetscErrorCode SVDSolve_TRLanczosGLower(SVD svd,BV U1,BV U2,BV V)
1415: {
1416: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1417: PetscReal *alpha,*beta,*alphah,*betah,normr,scalef,*sigma,sigma0;
1418: PetscScalar *w;
1419: PetscInt i,k,l,nv,ld;
1420: Mat U,Vmat,X;
1421: PetscBool breakdown=PETSC_FALSE,inverted;
1423: PetscFunctionBegin;
1424: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
1425: PetscCall(PetscMalloc2(ld,&w,ld,&sigma));
1426: inverted = ((svd->which==SVD_LARGEST && svd->swapped) || (svd->which==SVD_SMALLEST && !svd->swapped))? PETSC_TRUE: PETSC_FALSE;
1427: scalef = svd->swapped? 1.0/lanczos->scalef : lanczos->scalef;
1428: normr = (svd->conv==SVD_CONV_ABS)? PetscMax(svd->nrma,svd->nrmb*scalef): 1.0;
1430: /* normalize start vector */
1431: if (!svd->ninil) PetscCall(BVSetRandomColumn(U2,0));
1432: PetscCall(SVDInitialVectorGLower(svd,V,U1,U2,0,NULL));
1434: l = 0;
1435: while (svd->reason == SVD_CONVERGED_ITERATING) {
1436: svd->its++;
1438: /* inner loop */
1439: nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
1440: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
1441: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&alphah));
1442: beta = alpha + ld;
1443: betah = alpha + 2*ld;
1444: PetscCall(SVDLanczosGLower(svd,alpha,beta,alphah,betah,lanczos->Z,U1,U2,V,lanczos->ksp,svd->nconv+l,&nv,&breakdown));
1445: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
1446: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&alphah));
1447: PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
1448: PetscCall(BVSetActiveColumns(U1,svd->nconv,nv+1));
1449: PetscCall(BVSetActiveColumns(U2,svd->nconv,nv));
1451: /* solve projected problem */
1452: PetscCall(DSSetDimensions(svd->ds,nv+1,svd->nconv,svd->nconv+l));
1453: PetscCall(DSGSVDSetDimensions(svd->ds,nv,nv));
1454: PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
1455: PetscCall(DSSolve(svd->ds,w,NULL));
1456: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
1457: PetscCall(DSUpdateExtraRow(svd->ds));
1458: PetscCall(DSSynchronize(svd->ds,w,NULL));
1459: for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
1461: /* check convergence */
1462: PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,normr,&k));
1463: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
1465: sigma0 = inverted? 1.0/svd->sigma[0] : svd->sigma[0];
1466: if (lanczos->scaleth!=0 && k==0 && sigma0>lanczos->scaleth) {
1468: /* Scale and start from scratch */
1469: lanczos->scalef *= svd->swapped? 1.0/svd->sigma[0] : svd->sigma[0];
1470: PetscCall(PetscInfo(svd,"Scaling by factor %g and starting from scratch\n",(double)lanczos->scalef));
1471: PetscCall(MatZUpdateScale(svd));
1472: scalef = svd->swapped? 1.0/lanczos->scalef : lanczos->scalef;
1473: if (svd->conv==SVD_CONV_ABS) normr = PetscMax(svd->nrma,svd->nrmb*scalef);
1474: l = 0;
1475: if (!svd->ninil) PetscCall(BVSetRandomColumn(U2,0));
1476: PetscCall(SVDInitialVectorGLower(svd,V,U1,U2,0,NULL));
1478: } else {
1480: /* update l */
1481: if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
1482: else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
1483: if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
1484: if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
1486: if (svd->reason == SVD_CONVERGED_ITERATING) {
1487: if (PetscUnlikely(breakdown || k==nv)) {
1488: /* Start a new bidiagonalization */
1489: PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
1490: if (k<svd->nsv) {
1491: PetscCall(BVSetRandomColumn(U2,k));
1492: PetscCall(SVDInitialVectorGLower(svd,V,U1,U2,k,&breakdown));
1493: if (breakdown) {
1494: svd->reason = SVD_DIVERGED_BREAKDOWN;
1495: PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
1496: }
1497: }
1498: } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
1499: }
1501: /* compute converged singular vectors and restart vectors */
1502: PetscCall(DSGetMat(svd->ds,DS_MAT_X,&X));
1503: PetscCall(BVMultInPlace(V,X,svd->nconv,k+l));
1504: PetscCall(DSRestoreMat(svd->ds,DS_MAT_X,&X));
1505: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
1506: PetscCall(BVMultInPlace(U1,U,svd->nconv,k+l+1));
1507: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
1508: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&Vmat));
1509: PetscCall(BVMultInPlace(U2,Vmat,svd->nconv,k+l));
1510: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&Vmat));
1512: /* copy the last vector to be the next initial vector */
1513: if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(V,nv,k+l));
1514: }
1516: svd->nconv = k;
1517: PetscCall(SVDLanczosBackTransform(svd,nv,svd->sigma,sigma,NULL));
1518: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,sigma,svd->errest,nv));
1519: }
1521: PetscCall(PetscFree2(w,sigma));
1522: PetscFunctionReturn(PETSC_SUCCESS);
1523: }
1525: static PetscErrorCode SVDSolve_TRLanczos_GSVD(SVD svd)
1526: {
1527: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1528: PetscInt k,m,p;
1529: PetscBool convchg=PETSC_FALSE;
1530: BV U1,U2,UU;
1531: BVType type;
1532: VecType vtype;
1533: Mat U,V;
1534: SlepcSC sc;
1536: PetscFunctionBegin;
1537: PetscCall(PetscCitationsRegister(citationg,&citedg));
1539: if (svd->swapped) {
1540: PetscCall(DSGetSlepcSC(svd->ds,&sc));
1541: if (svd->which==SVD_LARGEST) sc->comparison = SlepcCompareSmallestReal;
1542: else sc->comparison = SlepcCompareLargestReal;
1543: }
1544: if (svd->converged==SVDConvergedNorm) { /* override temporarily since computed residual is already relative to the norms */
1545: svd->converged = SVDConvergedAbsolute;
1546: convchg = PETSC_TRUE;
1547: }
1548: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1549: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1551: /* Create BV for U1 */
1552: PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&U1));
1553: PetscCall(BVGetType(svd->U,&type));
1554: PetscCall(BVSetType(U1,type));
1555: PetscCall(BVGetSizes(svd->U,NULL,NULL,&k));
1556: PetscCall(BVSetSizes(U1,m,PETSC_DECIDE,k));
1557: PetscCall(BVGetVecType(svd->U,&vtype));
1558: PetscCall(BVSetVecType(U1,vtype));
1560: /* Create BV for U2 */
1561: PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&U2));
1562: PetscCall(BVSetType(U2,type));
1563: PetscCall(BVSetSizes(U2,p,PETSC_DECIDE,k));
1564: PetscCall(BVSetVecType(U2,vtype));
1566: /* Copy initial vectors from svd->U to U1 and U2 */
1567: if (svd->ninil) {
1568: Vec u, uh, nest, aux[2];
1569: PetscCall(BVGetColumn(U1,0,&u));
1570: PetscCall(BVGetColumn(U2,0,&uh));
1571: aux[0] = u;
1572: aux[1] = uh;
1573: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,aux,&nest));
1574: PetscCall(BVCopyVec(svd->U,0,nest));
1575: PetscCall(BVRestoreColumn(U1,0,&u));
1576: PetscCall(BVRestoreColumn(U2,0,&uh));
1577: PetscCall(VecDestroy(&nest));
1578: }
1580: switch (lanczos->bidiag) {
1581: case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1582: PetscCall(SVDSolve_TRLanczosGSingle(svd,U1,svd->U));
1583: break;
1584: case SVD_TRLANCZOS_GBIDIAG_UPPER:
1585: PetscCall(SVDSolve_TRLanczosGUpper(svd,U1,U2,svd->U));
1586: break;
1587: case SVD_TRLANCZOS_GBIDIAG_LOWER:
1588: PetscCall(SVDSolve_TRLanczosGLower(svd,U1,U2,svd->U));
1589: break;
1590: }
1592: /* Compute converged right singular vectors */
1593: PetscCall(BVSetActiveColumns(svd->U,0,svd->nconv));
1594: PetscCall(BVSetActiveColumns(svd->V,0,svd->nconv));
1595: PetscCall(BVGetMat(svd->U,&U));
1596: PetscCall(BVGetMat(svd->V,&V));
1597: PetscCall(KSPMatSolve(lanczos->ksp,U,V));
1598: PetscCall(BVRestoreMat(svd->U,&U));
1599: PetscCall(BVRestoreMat(svd->V,&V));
1601: /* Finish computing left singular vectors and move them to its place */
1602: if (svd->swapped) SWAP(U1,U2,UU);
1603: switch (lanczos->bidiag) {
1604: case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1605: PetscCall(SVDLeftSingularVectors_Single(svd,U1,U2));
1606: break;
1607: case SVD_TRLANCZOS_GBIDIAG_UPPER:
1608: case SVD_TRLANCZOS_GBIDIAG_LOWER:
1609: PetscCall(SVDLeftSingularVectors(svd,U1,U2));
1610: break;
1611: }
1613: /* undo scaling and compute the reciprocals of sigma if matrices were swapped */
1614: PetscCall(SVDLanczosBackTransform(svd,svd->nconv,svd->sigma,NULL,svd->V));
1616: PetscCall(BVDestroy(&U1));
1617: PetscCall(BVDestroy(&U2));
1618: PetscCall(DSTruncate(svd->ds,svd->nconv,PETSC_TRUE));
1619: if (convchg) svd->converged = SVDConvergedNorm;
1620: PetscFunctionReturn(PETSC_SUCCESS);
1621: }
1623: static PetscErrorCode SVDSetFromOptions_TRLanczos(SVD svd,PetscOptionItems *PetscOptionsObject)
1624: {
1625: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1626: PetscBool flg,val,lock;
1627: PetscReal keep,scale;
1628: SVDTRLanczosGBidiag bidiag;
1630: PetscFunctionBegin;
1631: PetscOptionsHeadBegin(PetscOptionsObject,"SVD TRLanczos Options");
1633: PetscCall(PetscOptionsBool("-svd_trlanczos_oneside","Use one-side reorthogonalization","SVDTRLanczosSetOneSide",lanczos->oneside,&val,&flg));
1634: if (flg) PetscCall(SVDTRLanczosSetOneSide(svd,val));
1636: PetscCall(PetscOptionsReal("-svd_trlanczos_restart","Proportion of vectors kept after restart","SVDTRLanczosSetRestart",0.5,&keep,&flg));
1637: if (flg) PetscCall(SVDTRLanczosSetRestart(svd,keep));
1639: PetscCall(PetscOptionsBool("-svd_trlanczos_locking","Choose between locking and non-locking variants","SVDTRLanczosSetLocking",PETSC_TRUE,&lock,&flg));
1640: if (flg) PetscCall(SVDTRLanczosSetLocking(svd,lock));
1642: PetscCall(PetscOptionsEnum("-svd_trlanczos_gbidiag","Bidiagonalization choice for Generalized Problem","SVDTRLanczosSetGBidiag",SVDTRLanczosGBidiags,(PetscEnum)lanczos->bidiag,(PetscEnum*)&bidiag,&flg));
1643: if (flg) PetscCall(SVDTRLanczosSetGBidiag(svd,bidiag));
1645: PetscCall(PetscOptionsBool("-svd_trlanczos_explicitmatrix","Build explicit matrix for KSP solver","SVDTRLanczosSetExplicitMatrix",lanczos->explicitmatrix,&val,&flg));
1646: if (flg) PetscCall(SVDTRLanczosSetExplicitMatrix(svd,val));
1648: PetscCall(SVDTRLanczosGetScale(svd,&scale));
1649: PetscCall(PetscOptionsReal("-svd_trlanczos_scale","Scale parameter for matrix B","SVDTRLanczosSetScale",scale,&scale,&flg));
1650: if (flg) PetscCall(SVDTRLanczosSetScale(svd,scale));
1652: PetscOptionsHeadEnd();
1654: if (svd->OPb) {
1655: if (!lanczos->ksp) PetscCall(SVDTRLanczosGetKSP(svd,&lanczos->ksp));
1656: PetscCall(KSPSetFromOptions(lanczos->ksp));
1657: }
1658: PetscFunctionReturn(PETSC_SUCCESS);
1659: }
1661: static PetscErrorCode SVDTRLanczosSetOneSide_TRLanczos(SVD svd,PetscBool oneside)
1662: {
1663: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1665: PetscFunctionBegin;
1666: if (lanczos->oneside != oneside) {
1667: lanczos->oneside = oneside;
1668: svd->state = SVD_STATE_INITIAL;
1669: }
1670: PetscFunctionReturn(PETSC_SUCCESS);
1671: }
1673: /*@
1674: SVDTRLanczosSetOneSide - Indicate if the variant of the Lanczos method
1675: to be used is one-sided or two-sided.
1677: Logically Collective
1679: Input Parameters:
1680: + svd - singular value solver
1681: - oneside - boolean flag indicating if the method is one-sided or not
1683: Options Database Key:
1684: . -svd_trlanczos_oneside <boolean> - Indicates the boolean flag
1686: Notes:
1687: By default, a two-sided variant is selected, which is sometimes slightly
1688: more robust. However, the one-sided variant is faster because it avoids
1689: the orthogonalization associated to left singular vectors.
1691: One-sided orthogonalization is also available for the GSVD, in which case
1692: two orthogonalizations out of three are avoided.
1694: Level: advanced
1696: .seealso: SVDLanczosSetOneSide()
1697: @*/
1698: PetscErrorCode SVDTRLanczosSetOneSide(SVD svd,PetscBool oneside)
1699: {
1700: PetscFunctionBegin;
1703: PetscTryMethod(svd,"SVDTRLanczosSetOneSide_C",(SVD,PetscBool),(svd,oneside));
1704: PetscFunctionReturn(PETSC_SUCCESS);
1705: }
1707: static PetscErrorCode SVDTRLanczosGetOneSide_TRLanczos(SVD svd,PetscBool *oneside)
1708: {
1709: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1711: PetscFunctionBegin;
1712: *oneside = lanczos->oneside;
1713: PetscFunctionReturn(PETSC_SUCCESS);
1714: }
1716: /*@
1717: SVDTRLanczosGetOneSide - Gets if the variant of the Lanczos method
1718: to be used is one-sided or two-sided.
1720: Not Collective
1722: Input Parameters:
1723: . svd - singular value solver
1725: Output Parameters:
1726: . oneside - boolean flag indicating if the method is one-sided or not
1728: Level: advanced
1730: .seealso: SVDTRLanczosSetOneSide()
1731: @*/
1732: PetscErrorCode SVDTRLanczosGetOneSide(SVD svd,PetscBool *oneside)
1733: {
1734: PetscFunctionBegin;
1736: PetscAssertPointer(oneside,2);
1737: PetscUseMethod(svd,"SVDTRLanczosGetOneSide_C",(SVD,PetscBool*),(svd,oneside));
1738: PetscFunctionReturn(PETSC_SUCCESS);
1739: }
1741: static PetscErrorCode SVDTRLanczosSetGBidiag_TRLanczos(SVD svd,SVDTRLanczosGBidiag bidiag)
1742: {
1743: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1745: PetscFunctionBegin;
1746: switch (bidiag) {
1747: case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1748: case SVD_TRLANCZOS_GBIDIAG_UPPER:
1749: case SVD_TRLANCZOS_GBIDIAG_LOWER:
1750: if (lanczos->bidiag != bidiag) {
1751: lanczos->bidiag = bidiag;
1752: svd->state = SVD_STATE_INITIAL;
1753: }
1754: break;
1755: default:
1756: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid bidiagonalization choice");
1757: }
1758: PetscFunctionReturn(PETSC_SUCCESS);
1759: }
1761: /*@
1762: SVDTRLanczosSetGBidiag - Sets the bidiagonalization choice to use in
1763: the GSVD TRLanczos solver.
1765: Logically Collective
1767: Input Parameters:
1768: + svd - the singular value solver
1769: - bidiag - the bidiagonalization choice
1771: Options Database Key:
1772: . -svd_trlanczos_gbidiag - Sets the bidiagonalization choice (either 's' or 'juu'
1773: or 'jlu')
1775: Level: advanced
1777: .seealso: SVDTRLanczosGetGBidiag(), SVDTRLanczosGBidiag
1778: @*/
1779: PetscErrorCode SVDTRLanczosSetGBidiag(SVD svd,SVDTRLanczosGBidiag bidiag)
1780: {
1781: PetscFunctionBegin;
1784: PetscTryMethod(svd,"SVDTRLanczosSetGBidiag_C",(SVD,SVDTRLanczosGBidiag),(svd,bidiag));
1785: PetscFunctionReturn(PETSC_SUCCESS);
1786: }
1788: static PetscErrorCode SVDTRLanczosGetGBidiag_TRLanczos(SVD svd,SVDTRLanczosGBidiag *bidiag)
1789: {
1790: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1792: PetscFunctionBegin;
1793: *bidiag = lanczos->bidiag;
1794: PetscFunctionReturn(PETSC_SUCCESS);
1795: }
1797: /*@
1798: SVDTRLanczosGetGBidiag - Gets the bidiagonalization choice used in the GSVD
1799: TRLanczos solver.
1801: Not Collective
1803: Input Parameter:
1804: . svd - the singular value solver
1806: Output Parameter:
1807: . bidiag - the bidiagonalization choice
1809: Level: advanced
1811: .seealso: SVDTRLanczosSetGBidiag(), SVDTRLanczosGBidiag
1812: @*/
1813: PetscErrorCode SVDTRLanczosGetGBidiag(SVD svd,SVDTRLanczosGBidiag *bidiag)
1814: {
1815: PetscFunctionBegin;
1817: PetscAssertPointer(bidiag,2);
1818: PetscUseMethod(svd,"SVDTRLanczosGetGBidiag_C",(SVD,SVDTRLanczosGBidiag*),(svd,bidiag));
1819: PetscFunctionReturn(PETSC_SUCCESS);
1820: }
1822: static PetscErrorCode SVDTRLanczosSetKSP_TRLanczos(SVD svd,KSP ksp)
1823: {
1824: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1826: PetscFunctionBegin;
1827: PetscCall(PetscObjectReference((PetscObject)ksp));
1828: PetscCall(KSPDestroy(&ctx->ksp));
1829: ctx->ksp = ksp;
1830: svd->state = SVD_STATE_INITIAL;
1831: PetscFunctionReturn(PETSC_SUCCESS);
1832: }
1834: /*@
1835: SVDTRLanczosSetKSP - Associate a linear solver object (KSP) to the SVD solver.
1837: Collective
1839: Input Parameters:
1840: + svd - SVD solver
1841: - ksp - the linear solver object
1843: Note:
1844: Only used for the GSVD problem.
1846: Level: advanced
1848: .seealso: SVDTRLanczosGetKSP()
1849: @*/
1850: PetscErrorCode SVDTRLanczosSetKSP(SVD svd,KSP ksp)
1851: {
1852: PetscFunctionBegin;
1855: PetscCheckSameComm(svd,1,ksp,2);
1856: PetscTryMethod(svd,"SVDTRLanczosSetKSP_C",(SVD,KSP),(svd,ksp));
1857: PetscFunctionReturn(PETSC_SUCCESS);
1858: }
1860: static PetscErrorCode SVDTRLanczosGetKSP_TRLanczos(SVD svd,KSP *ksp)
1861: {
1862: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1863: PC pc;
1865: PetscFunctionBegin;
1866: if (!ctx->ksp) {
1867: /* Create linear solver */
1868: PetscCall(KSPCreate(PetscObjectComm((PetscObject)svd),&ctx->ksp));
1869: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)svd,1));
1870: PetscCall(KSPSetOptionsPrefix(ctx->ksp,((PetscObject)svd)->prefix));
1871: PetscCall(KSPAppendOptionsPrefix(ctx->ksp,"svd_trlanczos_"));
1872: PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)svd)->options));
1873: PetscCall(KSPSetType(ctx->ksp,KSPLSQR));
1874: PetscCall(KSPGetPC(ctx->ksp,&pc));
1875: PetscCall(PCSetType(pc,PCNONE));
1876: PetscCall(KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE));
1877: PetscCall(KSPSetTolerances(ctx->ksp,SlepcDefaultTol(svd->tol)/10.0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
1878: }
1879: *ksp = ctx->ksp;
1880: PetscFunctionReturn(PETSC_SUCCESS);
1881: }
1883: /*@
1884: SVDTRLanczosGetKSP - Retrieve the linear solver object (KSP) associated with
1885: the SVD solver.
1887: Collective
1889: Input Parameter:
1890: . svd - SVD solver
1892: Output Parameter:
1893: . ksp - the linear solver object
1895: Level: advanced
1897: .seealso: SVDTRLanczosSetKSP()
1898: @*/
1899: PetscErrorCode SVDTRLanczosGetKSP(SVD svd,KSP *ksp)
1900: {
1901: PetscFunctionBegin;
1903: PetscAssertPointer(ksp,2);
1904: PetscUseMethod(svd,"SVDTRLanczosGetKSP_C",(SVD,KSP*),(svd,ksp));
1905: PetscFunctionReturn(PETSC_SUCCESS);
1906: }
1908: static PetscErrorCode SVDTRLanczosSetRestart_TRLanczos(SVD svd,PetscReal keep)
1909: {
1910: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1912: PetscFunctionBegin;
1913: if (keep==(PetscReal)PETSC_DEFAULT) ctx->keep = 0.5;
1914: else {
1915: 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);
1916: ctx->keep = keep;
1917: }
1918: PetscFunctionReturn(PETSC_SUCCESS);
1919: }
1921: /*@
1922: SVDTRLanczosSetRestart - Sets the restart parameter for the thick-restart
1923: Lanczos method, in particular the proportion of basis vectors that must be
1924: kept after restart.
1926: Logically Collective
1928: Input Parameters:
1929: + svd - the singular value solver
1930: - keep - the number of vectors to be kept at restart
1932: Options Database Key:
1933: . -svd_trlanczos_restart - Sets the restart parameter
1935: Notes:
1936: Allowed values are in the range [0.1,0.9]. The default is 0.5.
1938: Level: advanced
1940: .seealso: SVDTRLanczosGetRestart()
1941: @*/
1942: PetscErrorCode SVDTRLanczosSetRestart(SVD svd,PetscReal keep)
1943: {
1944: PetscFunctionBegin;
1947: PetscTryMethod(svd,"SVDTRLanczosSetRestart_C",(SVD,PetscReal),(svd,keep));
1948: PetscFunctionReturn(PETSC_SUCCESS);
1949: }
1951: static PetscErrorCode SVDTRLanczosGetRestart_TRLanczos(SVD svd,PetscReal *keep)
1952: {
1953: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1955: PetscFunctionBegin;
1956: *keep = ctx->keep;
1957: PetscFunctionReturn(PETSC_SUCCESS);
1958: }
1960: /*@
1961: SVDTRLanczosGetRestart - Gets the restart parameter used in the thick-restart
1962: Lanczos method.
1964: Not Collective
1966: Input Parameter:
1967: . svd - the singular value solver
1969: Output Parameter:
1970: . keep - the restart parameter
1972: Level: advanced
1974: .seealso: SVDTRLanczosSetRestart()
1975: @*/
1976: PetscErrorCode SVDTRLanczosGetRestart(SVD svd,PetscReal *keep)
1977: {
1978: PetscFunctionBegin;
1980: PetscAssertPointer(keep,2);
1981: PetscUseMethod(svd,"SVDTRLanczosGetRestart_C",(SVD,PetscReal*),(svd,keep));
1982: PetscFunctionReturn(PETSC_SUCCESS);
1983: }
1985: static PetscErrorCode SVDTRLanczosSetLocking_TRLanczos(SVD svd,PetscBool lock)
1986: {
1987: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1989: PetscFunctionBegin;
1990: ctx->lock = lock;
1991: PetscFunctionReturn(PETSC_SUCCESS);
1992: }
1994: /*@
1995: SVDTRLanczosSetLocking - Choose between locking and non-locking variants of
1996: the thick-restart Lanczos method.
1998: Logically Collective
2000: Input Parameters:
2001: + svd - the singular value solver
2002: - lock - true if the locking variant must be selected
2004: Options Database Key:
2005: . -svd_trlanczos_locking - Sets the locking flag
2007: Notes:
2008: The default is to lock converged singular triplets when the method restarts.
2009: This behaviour can be changed so that all directions are kept in the
2010: working subspace even if already converged to working accuracy (the
2011: non-locking variant).
2013: Level: advanced
2015: .seealso: SVDTRLanczosGetLocking()
2016: @*/
2017: PetscErrorCode SVDTRLanczosSetLocking(SVD svd,PetscBool lock)
2018: {
2019: PetscFunctionBegin;
2022: PetscTryMethod(svd,"SVDTRLanczosSetLocking_C",(SVD,PetscBool),(svd,lock));
2023: PetscFunctionReturn(PETSC_SUCCESS);
2024: }
2026: static PetscErrorCode SVDTRLanczosGetLocking_TRLanczos(SVD svd,PetscBool *lock)
2027: {
2028: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
2030: PetscFunctionBegin;
2031: *lock = ctx->lock;
2032: PetscFunctionReturn(PETSC_SUCCESS);
2033: }
2035: /*@
2036: SVDTRLanczosGetLocking - Gets the locking flag used in the thick-restart
2037: Lanczos method.
2039: Not Collective
2041: Input Parameter:
2042: . svd - the singular value solver
2044: Output Parameter:
2045: . lock - the locking flag
2047: Level: advanced
2049: .seealso: SVDTRLanczosSetLocking()
2050: @*/
2051: PetscErrorCode SVDTRLanczosGetLocking(SVD svd,PetscBool *lock)
2052: {
2053: PetscFunctionBegin;
2055: PetscAssertPointer(lock,2);
2056: PetscUseMethod(svd,"SVDTRLanczosGetLocking_C",(SVD,PetscBool*),(svd,lock));
2057: PetscFunctionReturn(PETSC_SUCCESS);
2058: }
2060: static PetscErrorCode SVDTRLanczosSetExplicitMatrix_TRLanczos(SVD svd,PetscBool explicitmat)
2061: {
2062: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2064: PetscFunctionBegin;
2065: if (lanczos->explicitmatrix != explicitmat) {
2066: lanczos->explicitmatrix = explicitmat;
2067: svd->state = SVD_STATE_INITIAL;
2068: }
2069: PetscFunctionReturn(PETSC_SUCCESS);
2070: }
2072: /*@
2073: SVDTRLanczosSetExplicitMatrix - Indicate if the matrix Z=[A;B] must
2074: be built explicitly.
2076: Logically Collective
2078: Input Parameters:
2079: + svd - singular value solver
2080: - explicitmat - Boolean flag indicating if Z=[A;B] is built explicitly
2082: Options Database Key:
2083: . -svd_trlanczos_explicitmatrix <boolean> - Indicates the boolean flag
2085: Notes:
2086: This option is relevant for the GSVD case only.
2087: Z is the coefficient matrix of the KSP solver used internally.
2089: Level: advanced
2091: .seealso: SVDTRLanczosGetExplicitMatrix()
2092: @*/
2093: PetscErrorCode SVDTRLanczosSetExplicitMatrix(SVD svd,PetscBool explicitmat)
2094: {
2095: PetscFunctionBegin;
2098: PetscTryMethod(svd,"SVDTRLanczosSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
2099: PetscFunctionReturn(PETSC_SUCCESS);
2100: }
2102: static PetscErrorCode SVDTRLanczosGetExplicitMatrix_TRLanczos(SVD svd,PetscBool *explicitmat)
2103: {
2104: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2106: PetscFunctionBegin;
2107: *explicitmat = lanczos->explicitmatrix;
2108: PetscFunctionReturn(PETSC_SUCCESS);
2109: }
2111: /*@
2112: SVDTRLanczosGetExplicitMatrix - Returns the flag indicating if Z=[A;B] is built explicitly.
2114: Not Collective
2116: Input Parameter:
2117: . svd - singular value solver
2119: Output Parameter:
2120: . explicitmat - the mode flag
2122: Level: advanced
2124: .seealso: SVDTRLanczosSetExplicitMatrix()
2125: @*/
2126: PetscErrorCode SVDTRLanczosGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
2127: {
2128: PetscFunctionBegin;
2130: PetscAssertPointer(explicitmat,2);
2131: PetscUseMethod(svd,"SVDTRLanczosGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
2132: PetscFunctionReturn(PETSC_SUCCESS);
2133: }
2135: static PetscErrorCode SVDTRLanczosSetScale_TRLanczos(SVD svd,PetscReal scale)
2136: {
2137: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
2139: PetscFunctionBegin;
2140: if (scale<0) {
2141: ctx->scalef = 1.0;
2142: ctx->scaleth = -scale;
2143: } else {
2144: ctx->scalef = scale;
2145: ctx->scaleth = 0.0;
2146: }
2147: PetscFunctionReturn(PETSC_SUCCESS);
2148: }
2150: /*@
2151: SVDTRLanczosSetScale - Sets the scale parameter for the GSVD.
2153: Logically Collective
2155: Input Parameters:
2156: + svd - singular value solver
2157: - scale - scale parameter
2159: Options Database Key:
2160: . -svd_trlanczos_scale <real> - scale factor/threshold
2162: Notes:
2163: This parameter is relevant for the GSVD case only. If the parameter is
2164: positive, it indicates the scale factor for B in matrix Z=[A;B]. If
2165: negative, its absolute value is the threshold for automatic scaling.
2166: In automatic scaling, whenever the largest approximate generalized singular
2167: value (or the inverse of the smallest value, if SVD_SMALLEST is used)
2168: exceeds the threshold, the computation is restarted with matrix B
2169: scaled by that value.
2171: Level: advanced
2173: .seealso: SVDTRLanczosGetScale()
2174: @*/
2175: PetscErrorCode SVDTRLanczosSetScale(SVD svd,PetscReal scale)
2176: {
2177: PetscFunctionBegin;
2180: PetscTryMethod(svd,"SVDTRLanczosSetScale_C",(SVD,PetscReal),(svd,scale));
2181: PetscFunctionReturn(PETSC_SUCCESS);
2182: }
2184: static PetscErrorCode SVDTRLanczosGetScale_TRLanczos(SVD svd,PetscReal *scale)
2185: {
2186: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
2188: PetscFunctionBegin;
2189: if (ctx->scaleth==0) *scale = ctx->scalef;
2190: else *scale = -ctx->scaleth;
2191: PetscFunctionReturn(PETSC_SUCCESS);
2192: }
2194: /*@
2195: SVDTRLanczosGetScale - Gets the scale parameter for the GSVD.
2197: Not Collective
2199: Input Parameter:
2200: . svd - the singular value solver
2202: Output Parameter:
2203: . scale - the scale parameter
2205: Notes:
2206: This parameter is relevant for the GSVD case only. If the parameter is
2207: positive, it indicates the scale factor for B in matrix Z=[A;B]. If
2208: negative, its absolute value is the threshold for automatic scaling.
2210: Level: advanced
2212: .seealso: SVDTRLanczosSetScale()
2213: @*/
2214: PetscErrorCode SVDTRLanczosGetScale(SVD svd,PetscReal *scale)
2215: {
2216: PetscFunctionBegin;
2218: PetscAssertPointer(scale,2);
2219: PetscUseMethod(svd,"SVDTRLanczosGetScale_C",(SVD,PetscReal*),(svd,scale));
2220: PetscFunctionReturn(PETSC_SUCCESS);
2221: }
2223: static PetscErrorCode SVDReset_TRLanczos(SVD svd)
2224: {
2225: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2227: PetscFunctionBegin;
2228: if (svd->isgeneralized || (!svd->problem_type && svd->OPb)) {
2229: PetscCall(KSPReset(lanczos->ksp));
2230: PetscCall(MatDestroy(&lanczos->Z));
2231: }
2232: PetscFunctionReturn(PETSC_SUCCESS);
2233: }
2235: static PetscErrorCode SVDDestroy_TRLanczos(SVD svd)
2236: {
2237: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2239: PetscFunctionBegin;
2240: if (svd->isgeneralized || (!svd->problem_type && svd->OPb)) PetscCall(KSPDestroy(&lanczos->ksp));
2241: PetscCall(PetscFree(svd->data));
2242: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",NULL));
2243: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",NULL));
2244: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetGBidiag_C",NULL));
2245: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetGBidiag_C",NULL));
2246: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetKSP_C",NULL));
2247: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetKSP_C",NULL));
2248: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetRestart_C",NULL));
2249: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetRestart_C",NULL));
2250: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetLocking_C",NULL));
2251: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetLocking_C",NULL));
2252: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetExplicitMatrix_C",NULL));
2253: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetExplicitMatrix_C",NULL));
2254: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetScale_C",NULL));
2255: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetScale_C",NULL));
2256: PetscFunctionReturn(PETSC_SUCCESS);
2257: }
2259: static PetscErrorCode SVDView_TRLanczos(SVD svd,PetscViewer viewer)
2260: {
2261: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2262: PetscBool isascii;
2264: PetscFunctionBegin;
2265: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
2266: if (isascii) {
2267: PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*lanczos->keep)));
2268: PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",lanczos->lock?"":"non-"));
2269: if (svd->isgeneralized) {
2270: const char *bidiag="";
2272: switch (lanczos->bidiag) {
2273: case SVD_TRLANCZOS_GBIDIAG_SINGLE: bidiag = "single"; break;
2274: case SVD_TRLANCZOS_GBIDIAG_UPPER: bidiag = "joint upper-upper"; break;
2275: case SVD_TRLANCZOS_GBIDIAG_LOWER: bidiag = "joint lower-upper"; break;
2276: }
2277: PetscCall(PetscViewerASCIIPrintf(viewer," bidiagonalization choice: %s\n",bidiag));
2278: PetscCall(PetscViewerASCIIPrintf(viewer," %s matrix\n",lanczos->explicitmatrix?"explicit":"implicit"));
2279: if (lanczos->scaleth==0) PetscCall(PetscViewerASCIIPrintf(viewer," scale factor for matrix B: %g\n",(double)lanczos->scalef));
2280: else PetscCall(PetscViewerASCIIPrintf(viewer," automatic scaling for matrix B with threshold: %g\n",(double)lanczos->scaleth));
2281: if (!lanczos->ksp) PetscCall(SVDTRLanczosGetKSP(svd,&lanczos->ksp));
2282: PetscCall(PetscViewerASCIIPushTab(viewer));
2283: PetscCall(KSPView(lanczos->ksp,viewer));
2284: PetscCall(PetscViewerASCIIPopTab(viewer));
2285: } else PetscCall(PetscViewerASCIIPrintf(viewer," %s-sided reorthogonalization\n",lanczos->oneside? "one": "two"));
2286: }
2287: PetscFunctionReturn(PETSC_SUCCESS);
2288: }
2290: static PetscErrorCode SVDSetDSType_TRLanczos(SVD svd)
2291: {
2292: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2293: DSType dstype;
2295: PetscFunctionBegin;
2296: dstype = svd->ishyperbolic? DSHSVD: DSSVD;
2297: if (svd->OPb && (lanczos->bidiag==SVD_TRLANCZOS_GBIDIAG_UPPER || lanczos->bidiag==SVD_TRLANCZOS_GBIDIAG_LOWER)) dstype = DSGSVD;
2298: PetscCall(DSSetType(svd->ds,dstype));
2299: PetscFunctionReturn(PETSC_SUCCESS);
2300: }
2302: SLEPC_EXTERN PetscErrorCode SVDCreate_TRLanczos(SVD svd)
2303: {
2304: SVD_TRLANCZOS *ctx;
2306: PetscFunctionBegin;
2307: PetscCall(PetscNew(&ctx));
2308: svd->data = (void*)ctx;
2310: ctx->lock = PETSC_TRUE;
2311: ctx->bidiag = SVD_TRLANCZOS_GBIDIAG_LOWER;
2312: ctx->scalef = 1.0;
2313: ctx->scaleth = 0.0;
2315: svd->ops->setup = SVDSetUp_TRLanczos;
2316: svd->ops->solve = SVDSolve_TRLanczos;
2317: svd->ops->solveg = SVDSolve_TRLanczos_GSVD;
2318: svd->ops->solveh = SVDSolve_TRLanczos_HSVD;
2319: svd->ops->destroy = SVDDestroy_TRLanczos;
2320: svd->ops->reset = SVDReset_TRLanczos;
2321: svd->ops->setfromoptions = SVDSetFromOptions_TRLanczos;
2322: svd->ops->view = SVDView_TRLanczos;
2323: svd->ops->setdstype = SVDSetDSType_TRLanczos;
2324: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",SVDTRLanczosSetOneSide_TRLanczos));
2325: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",SVDTRLanczosGetOneSide_TRLanczos));
2326: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetGBidiag_C",SVDTRLanczosSetGBidiag_TRLanczos));
2327: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetGBidiag_C",SVDTRLanczosGetGBidiag_TRLanczos));
2328: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetKSP_C",SVDTRLanczosSetKSP_TRLanczos));
2329: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetKSP_C",SVDTRLanczosGetKSP_TRLanczos));
2330: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetRestart_C",SVDTRLanczosSetRestart_TRLanczos));
2331: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetRestart_C",SVDTRLanczosGetRestart_TRLanczos));
2332: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetLocking_C",SVDTRLanczosSetLocking_TRLanczos));
2333: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetLocking_C",SVDTRLanczosGetLocking_TRLanczos));
2334: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetExplicitMatrix_C",SVDTRLanczosSetExplicitMatrix_TRLanczos));
2335: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetExplicitMatrix_C",SVDTRLanczosGetExplicitMatrix_TRLanczos));
2336: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetScale_C",SVDTRLanczosSetScale_TRLanczos));
2337: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetScale_C",SVDTRLanczosGetScale_TRLanczos));
2338: PetscFunctionReturn(PETSC_SUCCESS);
2339: }