Actual source code: trlanczos.c
slepc-main 2025-01-19
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)*((svd->stop==SVD_STOP_THRESHOLD)?10:1);
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,*aux;
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: SVDSetCtxThreshold(svd,svd->sigma,k);
504: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
506: /* update l */
507: if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
508: else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
509: if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
510: if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
512: if (svd->reason == SVD_CONVERGED_ITERATING) {
513: if (PetscUnlikely(breakdown || k==nv)) {
514: /* Start a new bidiagonalization */
515: PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
516: if (k<svd->nsv) {
517: PetscCall(BVSetRandomColumn(svd->V,k));
518: PetscCall(BVOrthonormalizeColumn(svd->V,k,PETSC_FALSE,NULL,&breakdown));
519: if (breakdown) {
520: svd->reason = SVD_DIVERGED_BREAKDOWN;
521: PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
522: }
523: }
524: } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
525: }
527: /* compute converged singular vectors and restart vectors */
528: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&V));
529: PetscCall(BVMultInPlace(svd->V,V,svd->nconv,k+l));
530: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&V));
531: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
532: PetscCall(BVMultInPlace(svd->U,U,svd->nconv,k+l));
533: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
535: if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) {
536: PetscCall(BVCopyColumn(svd->V,nv,k+l)); /* copy the last vector to be the next initial vector */
537: if (svd->stop==SVD_STOP_THRESHOLD && nv-k<5) { /* reallocate */
538: svd->ncv = svd->mpd+k;
539: PetscCall(SVDReallocateSolution(svd,svd->ncv+1));
540: for (i=nv;i<svd->ncv;i++) svd->perm[i] = i;
541: PetscCall(DSReallocate(svd->ds,svd->ncv+1));
542: aux = w;
543: PetscCall(PetscMalloc1(svd->ncv+1,&w));
544: PetscCall(PetscArraycpy(w,aux,ld));
545: PetscCall(PetscFree(aux));
546: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
547: if (lanczos->oneside) {
548: PetscCall(PetscFree(swork));
549: PetscCall(PetscMalloc1(svd->ncv+1,&swork));
550: }
551: }
552: }
554: svd->nconv = k;
555: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv));
556: }
558: /* orthonormalize U columns in one side method */
559: if (lanczos->oneside) {
560: for (i=0;i<svd->nconv;i++) PetscCall(BVOrthonormalizeColumn(svd->U,i,PETSC_FALSE,NULL,NULL));
561: }
563: /* free working space */
564: PetscCall(PetscFree(w));
565: if (swork) PetscCall(PetscFree(swork));
566: PetscCall(DSTruncate(svd->ds,svd->nconv,PETSC_TRUE));
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: 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)
571: {
572: PetscInt i;
573: Vec u,v,ou=svd->workl[0];
574: PetscBool lindep=PETSC_FALSE;
575: PetscReal norm;
577: PetscFunctionBegin;
578: for (i=k;i<*n;i++) {
579: PetscCall(BVGetColumn(V,i,&v));
580: PetscCall(BVGetColumn(U,i,&u));
581: PetscCall(MatMult(A,v,u));
582: PetscCall(BVRestoreColumn(V,i,&v));
583: PetscCall(BVRestoreColumn(U,i,&u));
584: PetscCall(BVOrthonormalizeColumn(U,i,PETSC_FALSE,alpha+i,&lindep));
585: omega[i] = PetscSign(alpha[i]);
586: if (PetscUnlikely(lindep)) {
587: *n = i;
588: break;
589: }
591: PetscCall(BVGetColumn(V,i+1,&v));
592: PetscCall(BVGetColumn(U,i,&u));
593: PetscCall(VecPointwiseMult(ou,svd->omega,u));
594: PetscCall(MatMult(AT,ou,v));
595: PetscCall(BVRestoreColumn(V,i+1,&v));
596: PetscCall(BVRestoreColumn(U,i,&u));
597: PetscCall(BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,&norm,&lindep));
598: beta[i] = omega[i]*norm;
599: if (PetscUnlikely(lindep)) {
600: *n = i+1;
601: break;
602: }
603: }
605: if (breakdown) *breakdown = lindep;
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }
609: static PetscErrorCode SVDSolve_TRLanczos_HSVD(SVD svd)
610: {
611: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
612: PetscReal *alpha,*beta,*omega;
613: PetscScalar *w,*aux;
614: PetscInt i,k,l,nv,ld,nini;
615: Mat UU,VV,D,A,AT;
616: BV U,V;
617: PetscBool breakdown=PETSC_FALSE;
618: BVOrthogType orthog;
619: Vec vomega;
621: PetscFunctionBegin;
622: /* undo the effect of swapping in this function */
623: if (svd->swapped) {
624: A = svd->AT;
625: AT = svd->A;
626: U = svd->V;
627: V = svd->U;
628: nini = svd->ninil;
629: } else {
630: A = svd->A;
631: AT = svd->AT;
632: U = svd->U;
633: V = svd->V;
634: nini = svd->nini;
635: }
636: /* allocate working space */
637: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
638: PetscCall(BVGetOrthogonalization(V,&orthog,NULL,NULL,NULL));
639: PetscCall(PetscMalloc1(ld,&w));
640: PetscCheck(!lanczos->oneside,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Oneside orthogonalization not supported for HSVD");
642: /* normalize start vector */
643: if (!nini) {
644: PetscCall(BVSetRandomColumn(V,0));
645: PetscCall(BVOrthonormalizeColumn(V,0,PETSC_TRUE,NULL,NULL));
646: }
648: l = 0;
649: while (svd->reason == SVD_CONVERGED_ITERATING) {
650: svd->its++;
652: /* inner loop */
653: nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
654: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
655: beta = alpha + ld;
656: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&omega));
657: PetscCall(SVDLanczosHSVD(svd,alpha,beta,omega,A,AT,V,U,svd->nconv+l,&nv,&breakdown));
658: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
659: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&omega));
660: PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
661: PetscCall(BVSetActiveColumns(U,svd->nconv,nv));
663: /* solve projected problem */
664: PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l));
665: PetscCall(DSHSVDSetDimensions(svd->ds,nv));
666: PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
667: PetscCall(DSSolve(svd->ds,w,NULL));
668: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
669: PetscCall(DSUpdateExtraRow(svd->ds));
670: PetscCall(DSSynchronize(svd->ds,w,NULL));
671: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&omega));
672: for (i=svd->nconv;i<nv;i++) {
673: svd->sigma[i] = PetscRealPart(w[i]);
674: svd->sign[i] = omega[i];
675: }
676: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&omega));
678: /* check convergence */
679: PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,1.0,&k));
680: SVDSetCtxThreshold(svd,svd->sigma,k);
681: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
683: /* update l */
684: if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
685: else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
686: if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
687: if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
689: if (svd->reason == SVD_CONVERGED_ITERATING) {
690: if (PetscUnlikely(breakdown || k==nv)) {
691: /* Start a new bidiagonalization */
692: PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
693: if (k<svd->nsv) {
694: PetscCall(BVSetRandomColumn(V,k));
695: PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,NULL,&breakdown));
696: if (breakdown) {
697: svd->reason = SVD_DIVERGED_BREAKDOWN;
698: PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
699: }
700: }
701: } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
702: }
704: /* compute converged singular vectors and restart vectors */
705: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&VV));
706: PetscCall(BVMultInPlace(V,VV,svd->nconv,k+l));
707: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&VV));
708: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&UU));
709: PetscCall(BVMultInPlace(U,UU,svd->nconv,k+l));
710: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&UU));
712: if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) {
713: /* copy the last vector of V to be the next initial vector
714: and change signature matrix of U */
715: PetscCall(BVCopyColumn(V,nv,k+l));
716: PetscCall(BVSetActiveColumns(U,0,k+l));
717: PetscCall(DSGetMatAndColumn(svd->ds,DS_MAT_D,0,&D,&vomega));
718: PetscCall(BVSetSignature(U,vomega));
719: PetscCall(DSRestoreMatAndColumn(svd->ds,DS_MAT_D,0,&D,&vomega));
721: if (svd->stop==SVD_STOP_THRESHOLD && nv-k<5) { /* reallocate */
722: svd->ncv = svd->mpd+k;
723: PetscCall(SVDReallocateSolution(svd,svd->ncv+1));
724: for (i=nv;i<svd->ncv;i++) svd->perm[i] = i;
725: PetscCall(DSReallocate(svd->ds,svd->ncv+1));
726: aux = w;
727: PetscCall(PetscMalloc1(svd->ncv+1,&w));
728: PetscCall(PetscArraycpy(w,aux,ld));
729: PetscCall(PetscFree(aux));
730: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
731: }
732: }
734: svd->nconv = k;
735: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv));
736: }
738: /* free working space */
739: PetscCall(PetscFree(w));
740: PetscCall(DSTruncate(svd->ds,svd->nconv,PETSC_TRUE));
741: PetscFunctionReturn(PETSC_SUCCESS);
742: }
744: /* Given n computed generalized singular values in sigmain, backtransform them
745: in sigmaout by undoing scaling and reciprocating if swapped=true. Also updates vectors V
746: if given. If sigmaout=NULL then the result overwrites sigmain. */
747: static PetscErrorCode SVDLanczosBackTransform(SVD svd,PetscInt n,PetscReal *sigmain,PetscReal *sigmaout,BV V)
748: {
749: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
750: PetscInt i;
751: PetscReal c,s,r,f,scalef;
753: PetscFunctionBegin;
754: scalef = svd->swapped? 1.0/lanczos->scalef: lanczos->scalef;
755: for (i=0;i<n;i++) {
756: if (V && scalef != 1.0) {
757: s = 1.0/PetscSqrtReal(1.0+sigmain[i]*sigmain[i]);
758: c = sigmain[i]*s;
759: r = s/scalef;
760: f = 1.0/PetscSqrtReal(c*c+r*r);
761: PetscCall(BVScaleColumn(V,i,f));
762: }
763: if (sigmaout) {
764: if (svd->swapped) sigmaout[i] = 1.0/(sigmain[i]*scalef);
765: else sigmaout[i] = sigmain[i]*scalef;
766: } else {
767: sigmain[i] *= scalef;
768: if (svd->swapped) sigmain[i] = 1.0/sigmain[i];
769: }
770: }
771: PetscFunctionReturn(PETSC_SUCCESS);
772: }
774: static PetscErrorCode SVDLanczosGSingle(SVD svd,PetscReal *alpha,PetscReal *beta,Mat Z,BV V,BV U,KSP ksp,PetscInt k,PetscInt *n,PetscBool *breakdown)
775: {
776: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
777: PetscInt i,j,m;
778: const PetscScalar *carr;
779: PetscScalar *arr;
780: Vec u,v,ut=svd->workl[0],x=svd->workr[0],v1,u1,u2;
781: PetscBool lindep=PETSC_FALSE;
783: PetscFunctionBegin;
784: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&v1));
785: PetscCall(BVGetColumn(V,k,&v));
786: PetscCall(BVGetColumn(U,k,&u));
788: /* Form ut=[u;0] */
789: PetscCall(VecZeroEntries(ut));
790: PetscCall(VecGetLocalSize(u,&m));
791: PetscCall(VecGetArrayRead(u,&carr));
792: PetscCall(VecGetArray(ut,&arr));
793: for (j=0; j<m; j++) arr[j] = carr[j];
794: PetscCall(VecRestoreArrayRead(u,&carr));
795: PetscCall(VecRestoreArray(ut,&arr));
797: /* Solve least squares problem */
798: PetscCall(KSPSolve(ksp,ut,x));
800: PetscCall(MatMult(Z,x,v));
802: PetscCall(BVRestoreColumn(U,k,&u));
803: PetscCall(BVRestoreColumn(V,k,&v));
804: PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,alpha+k,&lindep));
805: if (PetscUnlikely(lindep)) {
806: *n = k;
807: if (breakdown) *breakdown = lindep;
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: for (i=k+1; i<*n; i++) {
813: /* Compute vector i of BV U */
814: PetscCall(BVGetColumn(V,i-1,&v));
815: PetscCall(VecGetArray(v,&arr));
816: PetscCall(VecPlaceArray(v1,arr));
817: PetscCall(VecRestoreArray(v,&arr));
818: PetscCall(BVRestoreColumn(V,i-1,&v));
819: PetscCall(BVInsertVec(U,i,v1));
820: PetscCall(VecResetArray(v1));
821: PetscCall(BVOrthonormalizeColumn(U,i,PETSC_FALSE,beta+i-1,&lindep));
822: if (PetscUnlikely(lindep)) {
823: *n = i;
824: break;
825: }
827: /* Compute vector i of BV V */
829: PetscCall(BVGetColumn(V,i,&v));
830: PetscCall(BVGetColumn(U,i,&u));
832: /* Form ut=[u;0] */
833: PetscCall(VecGetArrayRead(u,&carr));
834: PetscCall(VecGetArray(ut,&arr));
835: for (j=0; j<m; j++) arr[j] = carr[j];
836: PetscCall(VecRestoreArrayRead(u,&carr));
837: PetscCall(VecRestoreArray(ut,&arr));
839: /* Solve least squares problem */
840: PetscCall(KSPSolve(ksp,ut,x));
842: PetscCall(MatMult(Z,x,v));
844: PetscCall(BVRestoreColumn(U,i,&u));
845: PetscCall(BVRestoreColumn(V,i,&v));
846: if (!lanczos->oneside || i==k+1) PetscCall(BVOrthonormalizeColumn(V,i,PETSC_FALSE,alpha+i,&lindep));
847: else { /* cheap computation of V[i], if restart (i==k+1) do a full reorthogonalization */
848: PetscCall(BVGetColumn(V,i,&u2));
849: PetscCall(BVGetColumn(V,i-1,&u1));
850: PetscCall(VecAXPY(u2,-beta[i-1],u1));
851: PetscCall(BVRestoreColumn(V,i-1,&u1));
852: PetscCall(VecNorm(u2,NORM_2,&alpha[i]));
853: if (alpha[i]==0.0) lindep = PETSC_TRUE;
854: else PetscCall(VecScale(u2,1.0/alpha[i]));
855: PetscCall(BVRestoreColumn(V,i,&u2));
856: }
857: if (PetscUnlikely(lindep)) {
858: *n = i;
859: break;
860: }
861: }
863: /* Compute vector n of BV U */
864: if (!lindep) {
865: PetscCall(BVGetColumn(V,*n-1,&v));
866: PetscCall(VecGetArray(v,&arr));
867: PetscCall(VecPlaceArray(v1,arr));
868: PetscCall(VecRestoreArray(v,&arr));
869: PetscCall(BVRestoreColumn(V,*n-1,&v));
870: PetscCall(BVInsertVec(U,*n,v1));
871: PetscCall(VecResetArray(v1));
872: PetscCall(BVOrthonormalizeColumn(U,*n,PETSC_FALSE,beta+*n-1,&lindep));
873: }
874: if (breakdown) *breakdown = lindep;
875: PetscCall(VecDestroy(&v1));
876: PetscFunctionReturn(PETSC_SUCCESS);
877: }
879: /* solve generalized problem with single bidiagonalization of Q_A */
880: static PetscErrorCode SVDSolve_TRLanczosGSingle(SVD svd,BV U1,BV V)
881: {
882: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
883: PetscReal *alpha,*beta,normr,scaleth,sigma0,*sigma,*aux2;
884: PetscScalar *w,*aux1;
885: PetscInt i,k,l,nv,ld;
886: Mat U,VV;
887: PetscBool breakdown=PETSC_FALSE;
889: PetscFunctionBegin;
890: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
891: PetscCall(PetscMalloc2(ld,&w,ld,&sigma));
892: normr = (svd->conv==SVD_CONV_ABS)? PetscMax(svd->nrma,svd->nrmb*lanczos->scalef): 1.0;
893: /* Convert scale threshold th=c/s to the corresponding c */
894: scaleth = (lanczos->scaleth!=0)? lanczos->scaleth/PetscSqrtReal(lanczos->scaleth*lanczos->scaleth+1): 0.0;
896: /* normalize start vector */
897: if (!svd->ninil) {
898: PetscCall(BVSetRandomColumn(U1,0));
899: PetscCall(BVOrthonormalizeColumn(U1,0,PETSC_TRUE,NULL,NULL));
900: }
902: l = 0;
903: while (svd->reason == SVD_CONVERGED_ITERATING) {
904: svd->its++;
906: /* inner loop */
907: nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
908: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
909: beta = alpha + ld;
910: PetscCall(SVDLanczosGSingle(svd,alpha,beta,lanczos->Z,V,U1,lanczos->ksp,svd->nconv+l,&nv,&breakdown));
911: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
912: PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
913: PetscCall(BVSetActiveColumns(U1,svd->nconv,nv));
915: /* solve projected problem */
916: PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l));
917: PetscCall(DSSVDSetDimensions(svd->ds,nv));
918: PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
919: PetscCall(DSSolve(svd->ds,w,NULL));
920: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
921: PetscCall(DSUpdateExtraRow(svd->ds));
922: PetscCall(DSSynchronize(svd->ds,w,NULL));
923: for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
925: /* check convergence */
926: PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,normr,&k));
927: PetscCall(SVDLanczosBackTransform(svd,nv,svd->sigma,sigma,NULL));
928: SVDSetCtxThreshold(svd,sigma,k);
929: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
931: sigma0 = svd->which==SVD_LARGEST? svd->sigma[0] : 1.0/svd->sigma[0];
932: if (scaleth!=0 && k==0 && sigma0>scaleth) {
934: /* Scale and start from scratch */
935: lanczos->scalef *= svd->sigma[0]/PetscSqrtReal(1-svd->sigma[0]*svd->sigma[0]);
936: PetscCall(PetscInfo(svd,"Scaling by factor %g and starting from scratch\n",(double)lanczos->scalef));
937: PetscCall(MatZUpdateScale(svd));
938: if (svd->conv==SVD_CONV_ABS) normr = PetscMax(svd->nrma,svd->nrmb*lanczos->scalef);
939: l = 0;
941: } else {
943: /* update l */
944: if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
945: else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
946: if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
947: if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
949: if (svd->reason == SVD_CONVERGED_ITERATING) {
950: if (PetscUnlikely(breakdown || k==nv)) {
951: /* Start a new bidiagonalization */
952: PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
953: if (k<svd->nsv) {
954: PetscCall(BVSetRandomColumn(U1,k));
955: PetscCall(BVOrthonormalizeColumn(U1,k,PETSC_FALSE,NULL,&breakdown));
956: if (breakdown) {
957: svd->reason = SVD_DIVERGED_BREAKDOWN;
958: PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
959: }
960: }
961: } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
962: }
964: /* compute converged singular vectors and restart vectors */
965: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
966: PetscCall(BVMultInPlace(V,U,svd->nconv,k+l));
967: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
968: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&VV));
969: PetscCall(BVMultInPlace(U1,VV,svd->nconv,k+l));
970: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&VV));
972: if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) {
973: PetscCall(BVCopyColumn(U1,nv,k+l)); /* copy the last vector to be the next initial vector */
974: if (svd->stop==SVD_STOP_THRESHOLD && nv-k<5) { /* reallocate */
975: svd->ncv = svd->mpd+k;
976: PetscCall(SVDReallocateSolution(svd,svd->ncv+1));
977: PetscCall(BVResize(V,svd->ncv+1,PETSC_TRUE));
978: PetscCall(BVResize(U1,svd->ncv+1,PETSC_TRUE));
979: for (i=nv;i<svd->ncv;i++) svd->perm[i] = i;
980: PetscCall(DSReallocate(svd->ds,svd->ncv+1));
981: aux1 = w;
982: aux2 = sigma;
983: PetscCall(PetscMalloc2(svd->ncv+1,&w,svd->ncv+1,&sigma));
984: PetscCall(PetscArraycpy(w,aux1,ld));
985: PetscCall(PetscArraycpy(sigma,aux2,ld));
986: PetscCall(PetscFree2(aux1,aux2));
987: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
988: }
989: }
990: }
992: svd->nconv = k;
993: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,sigma,svd->errest,nv));
994: }
996: PetscCall(PetscFree2(w,sigma));
997: PetscFunctionReturn(PETSC_SUCCESS);
998: }
1000: /* Move generalized left singular vectors (0..nconv) from U1 and U2 to its final destination svd->U (single variant) */
1001: static inline PetscErrorCode SVDLeftSingularVectors_Single(SVD svd,BV U1,BV U2)
1002: {
1003: PetscInt i,k,m,p;
1004: Vec u,u1,u2;
1005: PetscScalar *ua,*u2a;
1006: const PetscScalar *u1a;
1007: PetscReal s;
1009: PetscFunctionBegin;
1010: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1011: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1012: for (i=0;i<svd->nconv;i++) {
1013: PetscCall(BVGetColumn(U1,i,&u1));
1014: PetscCall(BVGetColumn(U2,i,&u2));
1015: PetscCall(BVGetColumn(svd->U,i,&u));
1016: PetscCall(VecGetArrayRead(u1,&u1a));
1017: PetscCall(VecGetArray(u,&ua));
1018: PetscCall(VecGetArray(u2,&u2a));
1019: /* Copy column from U1 to upper part of u */
1020: for (k=0;k<m;k++) ua[k] = u1a[k];
1021: /* Copy column from lower part of U to U2. Orthogonalize column in U2 and copy back to U */
1022: for (k=0;k<p;k++) u2a[k] = ua[m+k];
1023: PetscCall(VecRestoreArray(u2,&u2a));
1024: PetscCall(BVRestoreColumn(U2,i,&u2));
1025: PetscCall(BVOrthonormalizeColumn(U2,i,PETSC_FALSE,&s,NULL));
1026: PetscCall(BVGetColumn(U2,i,&u2));
1027: PetscCall(VecGetArray(u2,&u2a));
1028: for (k=0;k<p;k++) ua[m+k] = u2a[k];
1029: /* Update singular value */
1030: svd->sigma[i] /= s;
1031: PetscCall(VecRestoreArrayRead(u1,&u1a));
1032: PetscCall(VecRestoreArray(u,&ua));
1033: PetscCall(VecRestoreArray(u2,&u2a));
1034: PetscCall(BVRestoreColumn(U1,i,&u1));
1035: PetscCall(BVRestoreColumn(U2,i,&u2));
1036: PetscCall(BVRestoreColumn(svd->U,i,&u));
1037: }
1038: PetscFunctionReturn(PETSC_SUCCESS);
1039: }
1041: 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)
1042: {
1043: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1044: PetscInt i,j,m,p;
1045: const PetscScalar *carr;
1046: PetscScalar *arr,*u2arr;
1047: Vec u,v,ut=svd->workl[0],x=svd->workr[0],v1,u1,u2;
1048: PetscBool lindep=PETSC_FALSE,lindep1=PETSC_FALSE,lindep2=PETSC_FALSE;
1050: PetscFunctionBegin;
1051: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&v1));
1052: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1053: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1055: for (i=k; i<*n; i++) {
1056: /* Compute vector i of BV U1 */
1057: PetscCall(BVGetColumn(V,i,&v));
1058: PetscCall(VecGetArrayRead(v,&carr));
1059: PetscCall(VecPlaceArray(v1,carr));
1060: PetscCall(BVInsertVec(U1,i,v1));
1061: PetscCall(VecResetArray(v1));
1062: if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(U1,i,PETSC_FALSE,alpha+i,&lindep1));
1063: else { /* cheap computation of U1[i], if restart (i==k) do a full reorthogonalization */
1064: PetscCall(BVGetColumn(U1,i,&u2));
1065: if (i>0) {
1066: PetscCall(BVGetColumn(U1,i-1,&u1));
1067: PetscCall(VecAXPY(u2,-beta[i-1],u1));
1068: PetscCall(BVRestoreColumn(U1,i-1,&u1));
1069: }
1070: PetscCall(VecNorm(u2,NORM_2,&alpha[i]));
1071: if (alpha[i]==0.0) lindep = PETSC_TRUE;
1072: else PetscCall(VecScale(u2,1.0/alpha[i]));
1073: PetscCall(BVRestoreColumn(U1,i,&u2));
1074: }
1076: /* Compute vector i of BV U2 */
1077: PetscCall(BVGetColumn(U2,i,&u2));
1078: PetscCall(VecGetArray(u2,&u2arr));
1079: if (i%2) {
1080: for (j=0; j<p; j++) u2arr[j] = -carr[m+j];
1081: } else {
1082: for (j=0; j<p; j++) u2arr[j] = carr[m+j];
1083: }
1084: PetscCall(VecRestoreArray(u2,&u2arr));
1085: PetscCall(VecRestoreArrayRead(v,&carr));
1086: PetscCall(BVRestoreColumn(V,i,&v));
1087: if (lanczos->oneside && i>k) { /* cheap computation of U2[i], if restart (i==k) do a full reorthogonalization */
1088: if (i>0) {
1089: PetscCall(BVGetColumn(U2,i-1,&u1));
1090: PetscCall(VecAXPY(u2,(i%2)?betah[i-1]:-betah[i-1],u1));
1091: PetscCall(BVRestoreColumn(U2,i-1,&u1));
1092: }
1093: PetscCall(VecNorm(u2,NORM_2,&alphah[i]));
1094: if (alphah[i]==0.0) lindep = PETSC_TRUE;
1095: else PetscCall(VecScale(u2,1.0/alphah[i]));
1096: }
1097: PetscCall(BVRestoreColumn(U2,i,&u2));
1098: if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(U2,i,PETSC_FALSE,alphah+i,&lindep2));
1099: if (i%2) alphah[i] = -alphah[i];
1100: if (PetscUnlikely(lindep1 || lindep2)) {
1101: lindep = PETSC_TRUE;
1102: *n = i;
1103: break;
1104: }
1106: /* Compute vector i+1 of BV V */
1107: PetscCall(BVGetColumn(V,i+1,&v));
1108: /* Form ut=[u;0] */
1109: PetscCall(BVGetColumn(U1,i,&u));
1110: PetscCall(VecZeroEntries(ut));
1111: PetscCall(VecGetArrayRead(u,&carr));
1112: PetscCall(VecGetArray(ut,&arr));
1113: for (j=0; j<m; j++) arr[j] = carr[j];
1114: PetscCall(VecRestoreArrayRead(u,&carr));
1115: PetscCall(VecRestoreArray(ut,&arr));
1116: /* Solve least squares problem */
1117: PetscCall(KSPSolve(ksp,ut,x));
1118: PetscCall(MatMult(Z,x,v));
1119: PetscCall(BVRestoreColumn(U1,i,&u));
1120: PetscCall(BVRestoreColumn(V,i+1,&v));
1121: PetscCall(BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,beta+i,&lindep));
1122: betah[i] = -alpha[i]*beta[i]/alphah[i];
1123: if (PetscUnlikely(lindep)) {
1124: *n = i;
1125: break;
1126: }
1127: }
1128: if (breakdown) *breakdown = lindep;
1129: PetscCall(VecDestroy(&v1));
1130: PetscFunctionReturn(PETSC_SUCCESS);
1131: }
1133: /* generate random initial vector in column k for joint upper-upper bidiagonalization */
1134: static inline PetscErrorCode SVDInitialVectorGUpper(SVD svd,BV V,BV U1,PetscInt k,PetscBool *breakdown)
1135: {
1136: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1137: Vec u,v,ut=svd->workl[0],x=svd->workr[0];
1138: PetscInt m,j;
1139: PetscScalar *arr;
1140: const PetscScalar *carr;
1142: PetscFunctionBegin;
1143: /* Form ut=[u;0] where u is the k-th column of U1 */
1144: PetscCall(VecZeroEntries(ut));
1145: PetscCall(BVGetColumn(U1,k,&u));
1146: PetscCall(VecGetLocalSize(u,&m));
1147: PetscCall(VecGetArrayRead(u,&carr));
1148: PetscCall(VecGetArray(ut,&arr));
1149: for (j=0; j<m; j++) arr[j] = carr[j];
1150: PetscCall(VecRestoreArrayRead(u,&carr));
1151: PetscCall(VecRestoreArray(ut,&arr));
1152: PetscCall(BVRestoreColumn(U1,k,&u));
1153: /* Solve least squares problem Z*x=ut for x. Then set v=Z*x */
1154: PetscCall(KSPSolve(lanczos->ksp,ut,x));
1155: PetscCall(BVGetColumn(V,k,&v));
1156: PetscCall(MatMult(lanczos->Z,x,v));
1157: PetscCall(BVRestoreColumn(V,k,&v));
1158: if (breakdown) PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,NULL,breakdown));
1159: else PetscCall(BVOrthonormalizeColumn(V,k,PETSC_TRUE,NULL,NULL));
1160: PetscFunctionReturn(PETSC_SUCCESS);
1161: }
1163: /* solve generalized problem with joint upper-upper bidiagonalization */
1164: static PetscErrorCode SVDSolve_TRLanczosGUpper(SVD svd,BV U1,BV U2,BV V)
1165: {
1166: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1167: PetscReal *alpha,*beta,*alphah,*betah,normr,sigma0,*sigma,*aux2;
1168: PetscScalar *w,*aux1;
1169: PetscInt i,k,l,nv,ld;
1170: Mat U,Vmat,X;
1171: PetscBool breakdown=PETSC_FALSE;
1173: PetscFunctionBegin;
1174: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
1175: PetscCall(PetscMalloc2(ld,&w,ld,&sigma));
1176: normr = (svd->conv==SVD_CONV_ABS)? PetscMax(svd->nrma,svd->nrmb*lanczos->scalef): 1.0;
1178: /* normalize start vector */
1179: if (!svd->ninil) PetscCall(BVSetRandomColumn(U1,0));
1180: PetscCall(SVDInitialVectorGUpper(svd,V,U1,0,NULL));
1182: l = 0;
1183: while (svd->reason == SVD_CONVERGED_ITERATING) {
1184: svd->its++;
1186: /* inner loop */
1187: nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
1188: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
1189: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&alphah));
1190: beta = alpha + ld;
1191: betah = alpha + 2*ld;
1192: PetscCall(SVDLanczosGUpper(svd,alpha,beta,alphah,betah,lanczos->Z,U1,U2,V,lanczos->ksp,svd->nconv+l,&nv,&breakdown));
1193: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
1194: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&alphah));
1195: PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
1196: PetscCall(BVSetActiveColumns(U1,svd->nconv,nv));
1197: PetscCall(BVSetActiveColumns(U2,svd->nconv,nv));
1199: /* solve projected problem */
1200: PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l));
1201: PetscCall(DSGSVDSetDimensions(svd->ds,nv,nv));
1202: PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
1203: PetscCall(DSSolve(svd->ds,w,NULL));
1204: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
1205: PetscCall(DSUpdateExtraRow(svd->ds));
1206: PetscCall(DSSynchronize(svd->ds,w,NULL));
1207: for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
1209: /* check convergence */
1210: PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,normr,&k));
1211: PetscCall(SVDLanczosBackTransform(svd,nv,svd->sigma,sigma,NULL));
1212: SVDSetCtxThreshold(svd,sigma,k);
1213: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
1215: sigma0 = svd->which==SVD_LARGEST? svd->sigma[0] : 1.0/svd->sigma[0];
1216: if (lanczos->scaleth!=0 && k==0 && sigma0>lanczos->scaleth) {
1218: /* Scale and start from scratch */
1219: lanczos->scalef *= svd->sigma[0];
1220: PetscCall(PetscInfo(svd,"Scaling by factor %g and starting from scratch\n",(double)lanczos->scalef));
1221: PetscCall(MatZUpdateScale(svd));
1222: if (svd->conv==SVD_CONV_ABS) normr = PetscMax(svd->nrma,svd->nrmb*lanczos->scalef);
1223: l = 0;
1224: if (!svd->ninil) PetscCall(BVSetRandomColumn(U1,0));
1225: PetscCall(SVDInitialVectorGUpper(svd,V,U1,0,NULL));
1227: } else {
1229: /* update l */
1230: if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
1231: else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
1232: if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
1233: if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
1235: if (svd->reason == SVD_CONVERGED_ITERATING) {
1236: if (PetscUnlikely(breakdown || k==nv)) {
1237: /* Start a new bidiagonalization */
1238: PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
1239: if (k<svd->nsv) {
1240: PetscCall(BVSetRandomColumn(U1,k));
1241: PetscCall(SVDInitialVectorGUpper(svd,V,U1,k,&breakdown));
1242: if (breakdown) {
1243: svd->reason = SVD_DIVERGED_BREAKDOWN;
1244: PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
1245: }
1246: }
1247: } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
1248: }
1249: /* compute converged singular vectors and restart vectors */
1250: PetscCall(DSGetMat(svd->ds,DS_MAT_X,&X));
1251: PetscCall(BVMultInPlace(V,X,svd->nconv,k+l));
1252: PetscCall(DSRestoreMat(svd->ds,DS_MAT_X,&X));
1253: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
1254: PetscCall(BVMultInPlace(U1,U,svd->nconv,k+l));
1255: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
1256: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&Vmat));
1257: PetscCall(BVMultInPlace(U2,Vmat,svd->nconv,k+l));
1258: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&Vmat));
1260: if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) {
1261: PetscCall(BVCopyColumn(V,nv,k+l)); /* copy the last vector to be the next initial vector */
1262: if (svd->stop==SVD_STOP_THRESHOLD && nv-k<5) { /* reallocate */
1263: svd->ncv = svd->mpd+k;
1264: PetscCall(SVDReallocateSolution(svd,svd->ncv+1));
1265: PetscCall(BVResize(U1,svd->ncv+1,PETSC_TRUE));
1266: PetscCall(BVResize(U2,svd->ncv+1,PETSC_TRUE));
1267: for (i=nv;i<svd->ncv;i++) svd->perm[i] = i;
1268: PetscCall(DSReallocate(svd->ds,svd->ncv+1));
1269: aux1 = w;
1270: aux2 = sigma;
1271: PetscCall(PetscMalloc2(svd->ncv+1,&w,svd->ncv+1,&sigma));
1272: PetscCall(PetscArraycpy(w,aux1,ld));
1273: PetscCall(PetscArraycpy(sigma,aux2,ld));
1274: PetscCall(PetscFree2(aux1,aux2));
1275: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
1276: }
1277: }
1278: }
1280: svd->nconv = k;
1281: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,sigma,svd->errest,nv));
1282: }
1284: PetscCall(PetscFree2(w,sigma));
1285: PetscFunctionReturn(PETSC_SUCCESS);
1286: }
1288: /* Move generalized left singular vectors (0..nconv) from U1 and U2 to its final destination svd->U (upper and lower variants) */
1289: static inline PetscErrorCode SVDLeftSingularVectors(SVD svd,BV U1,BV U2)
1290: {
1291: PetscInt i,k,m,p;
1292: Vec u,u1,u2;
1293: PetscScalar *ua;
1294: const PetscScalar *u1a,*u2a;
1296: PetscFunctionBegin;
1297: PetscCall(BVGetSizes(U1,&m,NULL,NULL));
1298: PetscCall(BVGetSizes(U2,&p,NULL,NULL));
1299: for (i=0;i<svd->nconv;i++) {
1300: PetscCall(BVGetColumn(U1,i,&u1));
1301: PetscCall(BVGetColumn(U2,i,&u2));
1302: PetscCall(BVGetColumn(svd->U,i,&u));
1303: PetscCall(VecGetArrayRead(u1,&u1a));
1304: PetscCall(VecGetArrayRead(u2,&u2a));
1305: PetscCall(VecGetArray(u,&ua));
1306: /* Copy column from u1 to upper part of u */
1307: for (k=0;k<m;k++) ua[k] = u1a[k];
1308: /* Copy column from u2 to lower part of u */
1309: for (k=0;k<p;k++) ua[m+k] = u2a[k];
1310: PetscCall(VecRestoreArrayRead(u1,&u1a));
1311: PetscCall(VecRestoreArrayRead(u2,&u2a));
1312: PetscCall(VecRestoreArray(u,&ua));
1313: PetscCall(BVRestoreColumn(U1,i,&u1));
1314: PetscCall(BVRestoreColumn(U2,i,&u2));
1315: PetscCall(BVRestoreColumn(svd->U,i,&u));
1316: }
1317: PetscFunctionReturn(PETSC_SUCCESS);
1318: }
1320: 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)
1321: {
1322: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1323: PetscInt i,j,m,p;
1324: const PetscScalar *carr;
1325: PetscScalar *arr,*u2arr;
1326: Vec u,v,ut=svd->workl[0],x=svd->workr[0],v1,u1,u2;
1327: PetscBool lindep=PETSC_FALSE;
1329: PetscFunctionBegin;
1330: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&v1));
1331: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1332: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1334: for (i=k; i<*n; i++) {
1335: /* Compute vector i of BV U2 */
1336: PetscCall(BVGetColumn(V,i,&v));
1337: PetscCall(VecGetArrayRead(v,&carr));
1338: PetscCall(BVGetColumn(U2,i,&u2));
1339: PetscCall(VecGetArray(u2,&u2arr));
1340: if (i%2) {
1341: for (j=0; j<p; j++) u2arr[j] = -carr[m+j];
1342: } else {
1343: for (j=0; j<p; j++) u2arr[j] = carr[m+j];
1344: }
1345: PetscCall(VecRestoreArray(u2,&u2arr));
1346: if (lanczos->oneside && i>k) { /* cheap computation of U2[i], if restart (i==k) do a full reorthogonalization */
1347: if (i>0) {
1348: PetscCall(BVGetColumn(U2,i-1,&u1));
1349: PetscCall(VecAXPY(u2,(i%2)?betah[i-1]:-betah[i-1],u1));
1350: PetscCall(BVRestoreColumn(U2,i-1,&u1));
1351: }
1352: PetscCall(VecNorm(u2,NORM_2,&alphah[i]));
1353: if (alphah[i]==0.0) lindep = PETSC_TRUE;
1354: else PetscCall(VecScale(u2,1.0/alphah[i]));
1355: }
1356: PetscCall(BVRestoreColumn(U2,i,&u2));
1357: if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(U2,i,PETSC_FALSE,alphah+i,&lindep));
1358: if (i%2) alphah[i] = -alphah[i];
1359: if (PetscUnlikely(lindep)) {
1360: PetscCall(BVRestoreColumn(V,i,&v));
1361: *n = i;
1362: break;
1363: }
1365: /* Compute vector i+1 of BV U1 */
1366: PetscCall(VecPlaceArray(v1,carr));
1367: PetscCall(BVInsertVec(U1,i+1,v1));
1368: PetscCall(VecResetArray(v1));
1369: PetscCall(BVOrthonormalizeColumn(U1,i+1,PETSC_FALSE,beta+i,&lindep));
1370: PetscCall(VecRestoreArrayRead(v,&carr));
1371: PetscCall(BVRestoreColumn(V,i,&v));
1372: if (PetscUnlikely(lindep)) {
1373: *n = i+1;
1374: break;
1375: }
1377: /* Compute vector i+1 of BV V */
1378: PetscCall(BVGetColumn(V,i+1,&v));
1379: /* Form ut=[u;0] where u is column i+1 of BV U1 */
1380: PetscCall(BVGetColumn(U1,i+1,&u));
1381: PetscCall(VecZeroEntries(ut));
1382: PetscCall(VecGetArrayRead(u,&carr));
1383: PetscCall(VecGetArray(ut,&arr));
1384: for (j=0; j<m; j++) arr[j] = carr[j];
1385: PetscCall(VecRestoreArrayRead(u,&carr));
1386: PetscCall(VecRestoreArray(ut,&arr));
1387: /* Solve least squares problem */
1388: PetscCall(KSPSolve(ksp,ut,x));
1389: PetscCall(MatMult(Z,x,v));
1390: PetscCall(BVRestoreColumn(U1,i+1,&u));
1391: PetscCall(BVRestoreColumn(V,i+1,&v));
1392: if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,alpha+i+1,&lindep));
1393: else { /* cheap computation of V[i+1], if restart (i==k) do a full reorthogonalization */
1394: PetscCall(BVGetColumn(V,i+1,&u2));
1395: PetscCall(BVGetColumn(V,i,&u1));
1396: PetscCall(VecAXPY(u2,-beta[i],u1));
1397: PetscCall(BVRestoreColumn(V,i,&u1));
1398: PetscCall(VecNorm(u2,NORM_2,&alpha[i+1]));
1399: if (alpha[i+1]==0.0) lindep = PETSC_TRUE;
1400: else PetscCall(VecScale(u2,1.0/alpha[i+1]));
1401: PetscCall(BVRestoreColumn(V,i+1,&u2));
1402: }
1403: betah[i] = -alpha[i+1]*beta[i]/alphah[i];
1404: if (PetscUnlikely(lindep)) {
1405: *n = i+1;
1406: break;
1407: }
1408: }
1409: if (breakdown) *breakdown = lindep;
1410: PetscCall(VecDestroy(&v1));
1411: PetscFunctionReturn(PETSC_SUCCESS);
1412: }
1414: /* generate random initial vector in column k for joint lower-upper bidiagonalization */
1415: static inline PetscErrorCode SVDInitialVectorGLower(SVD svd,BV V,BV U1,BV U2,PetscInt k,PetscBool *breakdown)
1416: {
1417: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1418: const PetscScalar *carr;
1419: PetscScalar *arr;
1420: PetscReal *alpha;
1421: PetscInt j,m,p;
1422: Vec u,uh,v,ut=svd->workl[0],x=svd->workr[0];
1424: PetscFunctionBegin;
1425: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1426: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1427: /* Form ut=[0;uh], where uh is the k-th column of U2 */
1428: PetscCall(BVGetColumn(U2,k,&uh));
1429: PetscCall(VecZeroEntries(ut));
1430: PetscCall(VecGetArrayRead(uh,&carr));
1431: PetscCall(VecGetArray(ut,&arr));
1432: for (j=0; j<p; j++) arr[m+j] = carr[j];
1433: PetscCall(VecRestoreArrayRead(uh,&carr));
1434: PetscCall(VecRestoreArray(ut,&arr));
1435: PetscCall(BVRestoreColumn(U2,k,&uh));
1436: /* Solve least squares problem Z*x=ut for x. Then set ut=Z*x */
1437: PetscCall(KSPSolve(lanczos->ksp,ut,x));
1438: PetscCall(MatMult(lanczos->Z,x,ut));
1439: /* Form u, column k of BV U1, as the upper part of ut and orthonormalize */
1440: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&u));
1441: PetscCall(VecGetArrayRead(ut,&carr));
1442: PetscCall(VecPlaceArray(u,carr));
1443: PetscCall(BVInsertVec(U1,k,u));
1444: PetscCall(VecResetArray(u));
1445: PetscCall(VecRestoreArrayRead(ut,&carr));
1446: PetscCall(VecDestroy(&u));
1447: if (breakdown) PetscCall(BVOrthonormalizeColumn(U1,k,PETSC_FALSE,NULL,breakdown));
1448: else PetscCall(BVOrthonormalizeColumn(U1,k,PETSC_TRUE,NULL,NULL));
1450: if (!breakdown || !*breakdown) {
1451: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1452: /* Compute k-th vector of BV V */
1453: PetscCall(BVGetColumn(V,k,&v));
1454: /* Form ut=[u;0] where u is the 1st column of U1 */
1455: PetscCall(BVGetColumn(U1,k,&u));
1456: PetscCall(VecZeroEntries(ut));
1457: PetscCall(VecGetArrayRead(u,&carr));
1458: PetscCall(VecGetArray(ut,&arr));
1459: for (j=0; j<m; j++) arr[j] = carr[j];
1460: PetscCall(VecRestoreArrayRead(u,&carr));
1461: PetscCall(VecRestoreArray(ut,&arr));
1462: /* Solve least squares problem */
1463: PetscCall(KSPSolve(lanczos->ksp,ut,x));
1464: PetscCall(MatMult(lanczos->Z,x,v));
1465: PetscCall(BVRestoreColumn(U1,k,&u));
1466: PetscCall(BVRestoreColumn(V,k,&v));
1467: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
1468: if (breakdown) PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,alpha+k,breakdown));
1469: else PetscCall(BVOrthonormalizeColumn(V,k,PETSC_TRUE,alpha+k,NULL));
1470: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
1471: }
1472: PetscFunctionReturn(PETSC_SUCCESS);
1473: }
1475: /* solve generalized problem with joint lower-upper bidiagonalization */
1476: static PetscErrorCode SVDSolve_TRLanczosGLower(SVD svd,BV U1,BV U2,BV V)
1477: {
1478: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1479: PetscReal *alpha,*beta,*alphah,*betah,normr,scalef,*sigma,sigma0,*aux2;
1480: PetscScalar *w,*aux1;
1481: PetscInt i,k,l,nv,ld;
1482: Mat U,Vmat,X;
1483: PetscBool breakdown=PETSC_FALSE,inverted;
1485: PetscFunctionBegin;
1486: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
1487: PetscCall(PetscMalloc2(ld,&w,ld,&sigma));
1488: inverted = ((svd->which==SVD_LARGEST && svd->swapped) || (svd->which==SVD_SMALLEST && !svd->swapped))? PETSC_TRUE: PETSC_FALSE;
1489: scalef = svd->swapped? 1.0/lanczos->scalef : lanczos->scalef;
1490: normr = (svd->conv==SVD_CONV_ABS)? PetscMax(svd->nrma,svd->nrmb*scalef): 1.0;
1492: /* normalize start vector */
1493: if (!svd->ninil) PetscCall(BVSetRandomColumn(U2,0));
1494: PetscCall(SVDInitialVectorGLower(svd,V,U1,U2,0,NULL));
1496: l = 0;
1497: while (svd->reason == SVD_CONVERGED_ITERATING) {
1498: svd->its++;
1500: /* inner loop */
1501: nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
1502: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
1503: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&alphah));
1504: beta = alpha + ld;
1505: betah = alpha + 2*ld;
1506: PetscCall(SVDLanczosGLower(svd,alpha,beta,alphah,betah,lanczos->Z,U1,U2,V,lanczos->ksp,svd->nconv+l,&nv,&breakdown));
1507: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
1508: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&alphah));
1509: PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
1510: PetscCall(BVSetActiveColumns(U1,svd->nconv,nv+1));
1511: PetscCall(BVSetActiveColumns(U2,svd->nconv,nv));
1513: /* solve projected problem */
1514: PetscCall(DSSetDimensions(svd->ds,nv+1,svd->nconv,svd->nconv+l));
1515: PetscCall(DSGSVDSetDimensions(svd->ds,nv,nv));
1516: PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
1517: PetscCall(DSSolve(svd->ds,w,NULL));
1518: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
1519: PetscCall(DSUpdateExtraRow(svd->ds));
1520: PetscCall(DSSynchronize(svd->ds,w,NULL));
1521: for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
1523: /* check convergence */
1524: PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,normr,&k));
1525: PetscCall(SVDLanczosBackTransform(svd,nv,svd->sigma,sigma,NULL));
1526: SVDSetCtxThreshold(svd,sigma,k);
1527: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
1529: sigma0 = inverted? 1.0/svd->sigma[0] : svd->sigma[0];
1530: if (lanczos->scaleth!=0 && k==0 && sigma0>lanczos->scaleth) {
1532: /* Scale and start from scratch */
1533: lanczos->scalef *= svd->swapped? 1.0/svd->sigma[0] : svd->sigma[0];
1534: PetscCall(PetscInfo(svd,"Scaling by factor %g and starting from scratch\n",(double)lanczos->scalef));
1535: PetscCall(MatZUpdateScale(svd));
1536: scalef = svd->swapped? 1.0/lanczos->scalef : lanczos->scalef;
1537: if (svd->conv==SVD_CONV_ABS) normr = PetscMax(svd->nrma,svd->nrmb*scalef);
1538: l = 0;
1539: if (!svd->ninil) PetscCall(BVSetRandomColumn(U2,0));
1540: PetscCall(SVDInitialVectorGLower(svd,V,U1,U2,0,NULL));
1542: } else {
1544: /* update l */
1545: if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
1546: else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
1547: if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
1548: if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
1550: if (svd->reason == SVD_CONVERGED_ITERATING) {
1551: if (PetscUnlikely(breakdown || k==nv)) {
1552: /* Start a new bidiagonalization */
1553: PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
1554: if (k<svd->nsv) {
1555: PetscCall(BVSetRandomColumn(U2,k));
1556: PetscCall(SVDInitialVectorGLower(svd,V,U1,U2,k,&breakdown));
1557: if (breakdown) {
1558: svd->reason = SVD_DIVERGED_BREAKDOWN;
1559: PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
1560: }
1561: }
1562: } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
1563: }
1565: /* compute converged singular vectors and restart vectors */
1566: PetscCall(DSGetMat(svd->ds,DS_MAT_X,&X));
1567: PetscCall(BVMultInPlace(V,X,svd->nconv,k+l));
1568: PetscCall(DSRestoreMat(svd->ds,DS_MAT_X,&X));
1569: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
1570: PetscCall(BVMultInPlace(U1,U,svd->nconv,k+l+1));
1571: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
1572: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&Vmat));
1573: PetscCall(BVMultInPlace(U2,Vmat,svd->nconv,k+l));
1574: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&Vmat));
1576: if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) {
1577: PetscCall(BVCopyColumn(V,nv,k+l)); /* copy the last vector to be the next initial vector */
1578: if (svd->stop==SVD_STOP_THRESHOLD && nv-k<5) { /* reallocate */
1579: svd->ncv = svd->mpd+k;
1580: PetscCall(SVDReallocateSolution(svd,svd->ncv+1));
1581: PetscCall(BVResize(U1,svd->ncv+1,PETSC_TRUE));
1582: PetscCall(BVResize(U2,svd->ncv+1,PETSC_TRUE));
1583: for (i=nv;i<svd->ncv;i++) svd->perm[i] = i;
1584: PetscCall(DSReallocate(svd->ds,svd->ncv+1));
1585: aux1 = w;
1586: aux2 = sigma;
1587: PetscCall(PetscMalloc2(svd->ncv+1,&w,svd->ncv+1,&sigma));
1588: PetscCall(PetscArraycpy(w,aux1,ld));
1589: PetscCall(PetscArraycpy(sigma,aux2,ld));
1590: PetscCall(PetscFree2(aux1,aux2));
1591: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
1592: }
1593: }
1594: }
1596: svd->nconv = k;
1597: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,sigma,svd->errest,nv));
1598: }
1600: PetscCall(PetscFree2(w,sigma));
1601: PetscFunctionReturn(PETSC_SUCCESS);
1602: }
1604: static PetscErrorCode SVDSolve_TRLanczos_GSVD(SVD svd)
1605: {
1606: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1607: PetscInt k,m,p;
1608: PetscBool convchg=PETSC_FALSE;
1609: BV U1,U2,UU;
1610: BVType type;
1611: VecType vtype;
1612: Mat U,V;
1613: SlepcSC sc;
1615: PetscFunctionBegin;
1616: PetscCall(PetscCitationsRegister(citationg,&citedg));
1618: if (svd->swapped) {
1619: PetscCall(DSGetSlepcSC(svd->ds,&sc));
1620: if (svd->which==SVD_LARGEST) sc->comparison = SlepcCompareSmallestReal;
1621: else sc->comparison = SlepcCompareLargestReal;
1622: }
1623: if (svd->converged==SVDConvergedNorm) { /* override temporarily since computed residual is already relative to the norms */
1624: svd->converged = SVDConvergedAbsolute;
1625: convchg = PETSC_TRUE;
1626: }
1627: PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1628: PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1630: /* Create BV for U1 */
1631: PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&U1));
1632: PetscCall(BVGetType(svd->U,&type));
1633: PetscCall(BVSetType(U1,type));
1634: PetscCall(BVGetSizes(svd->U,NULL,NULL,&k));
1635: PetscCall(BVSetSizes(U1,m,PETSC_DECIDE,k));
1636: PetscCall(BVGetVecType(svd->U,&vtype));
1637: PetscCall(BVSetVecType(U1,vtype));
1639: /* Create BV for U2 */
1640: PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&U2));
1641: PetscCall(BVSetType(U2,type));
1642: PetscCall(BVSetSizes(U2,p,PETSC_DECIDE,k));
1643: PetscCall(BVSetVecType(U2,vtype));
1645: /* Copy initial vectors from svd->U to U1 and U2 */
1646: if (svd->ninil) {
1647: Vec u, uh, nest, aux[2];
1648: PetscCall(BVGetColumn(U1,0,&u));
1649: PetscCall(BVGetColumn(U2,0,&uh));
1650: aux[0] = u;
1651: aux[1] = uh;
1652: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,aux,&nest));
1653: PetscCall(BVCopyVec(svd->U,0,nest));
1654: PetscCall(BVRestoreColumn(U1,0,&u));
1655: PetscCall(BVRestoreColumn(U2,0,&uh));
1656: PetscCall(VecDestroy(&nest));
1657: }
1659: switch (lanczos->bidiag) {
1660: case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1661: PetscCall(SVDSolve_TRLanczosGSingle(svd,U1,svd->U));
1662: if (svd->stop==SVD_STOP_THRESHOLD) PetscCall(BVResize(U2,svd->ncv+1,PETSC_FALSE));
1663: break;
1664: case SVD_TRLANCZOS_GBIDIAG_UPPER:
1665: PetscCall(SVDSolve_TRLanczosGUpper(svd,U1,U2,svd->U));
1666: break;
1667: case SVD_TRLANCZOS_GBIDIAG_LOWER:
1668: PetscCall(SVDSolve_TRLanczosGLower(svd,U1,U2,svd->U));
1669: break;
1670: }
1672: /* Compute converged right singular vectors */
1673: PetscCall(BVSetActiveColumns(svd->U,0,svd->nconv));
1674: PetscCall(BVSetActiveColumns(svd->V,0,svd->nconv));
1675: PetscCall(BVGetMat(svd->U,&U));
1676: PetscCall(BVGetMat(svd->V,&V));
1677: PetscCall(KSPMatSolve(lanczos->ksp,U,V));
1678: PetscCall(BVRestoreMat(svd->U,&U));
1679: PetscCall(BVRestoreMat(svd->V,&V));
1681: /* Finish computing left singular vectors and move them to its place */
1682: if (svd->swapped) SlepcSwap(U1,U2,UU);
1683: switch (lanczos->bidiag) {
1684: case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1685: PetscCall(SVDLeftSingularVectors_Single(svd,U1,U2));
1686: break;
1687: case SVD_TRLANCZOS_GBIDIAG_UPPER:
1688: case SVD_TRLANCZOS_GBIDIAG_LOWER:
1689: PetscCall(SVDLeftSingularVectors(svd,U1,U2));
1690: break;
1691: }
1693: /* undo scaling and compute the reciprocals of sigma if matrices were swapped */
1694: PetscCall(SVDLanczosBackTransform(svd,svd->nconv,svd->sigma,NULL,svd->V));
1696: PetscCall(BVDestroy(&U1));
1697: PetscCall(BVDestroy(&U2));
1698: PetscCall(DSTruncate(svd->ds,svd->nconv,PETSC_TRUE));
1699: if (convchg) svd->converged = SVDConvergedNorm;
1700: PetscFunctionReturn(PETSC_SUCCESS);
1701: }
1703: static PetscErrorCode SVDSetFromOptions_TRLanczos(SVD svd,PetscOptionItems *PetscOptionsObject)
1704: {
1705: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1706: PetscBool flg,val,lock;
1707: PetscReal keep,scale;
1708: SVDTRLanczosGBidiag bidiag;
1710: PetscFunctionBegin;
1711: PetscOptionsHeadBegin(PetscOptionsObject,"SVD TRLanczos Options");
1713: PetscCall(PetscOptionsBool("-svd_trlanczos_oneside","Use one-side reorthogonalization","SVDTRLanczosSetOneSide",lanczos->oneside,&val,&flg));
1714: if (flg) PetscCall(SVDTRLanczosSetOneSide(svd,val));
1716: PetscCall(PetscOptionsReal("-svd_trlanczos_restart","Proportion of vectors kept after restart","SVDTRLanczosSetRestart",0.5,&keep,&flg));
1717: if (flg) PetscCall(SVDTRLanczosSetRestart(svd,keep));
1719: PetscCall(PetscOptionsBool("-svd_trlanczos_locking","Choose between locking and non-locking variants","SVDTRLanczosSetLocking",PETSC_TRUE,&lock,&flg));
1720: if (flg) PetscCall(SVDTRLanczosSetLocking(svd,lock));
1722: PetscCall(PetscOptionsEnum("-svd_trlanczos_gbidiag","Bidiagonalization choice for Generalized Problem","SVDTRLanczosSetGBidiag",SVDTRLanczosGBidiags,(PetscEnum)lanczos->bidiag,(PetscEnum*)&bidiag,&flg));
1723: if (flg) PetscCall(SVDTRLanczosSetGBidiag(svd,bidiag));
1725: PetscCall(PetscOptionsBool("-svd_trlanczos_explicitmatrix","Build explicit matrix for KSP solver","SVDTRLanczosSetExplicitMatrix",lanczos->explicitmatrix,&val,&flg));
1726: if (flg) PetscCall(SVDTRLanczosSetExplicitMatrix(svd,val));
1728: PetscCall(SVDTRLanczosGetScale(svd,&scale));
1729: PetscCall(PetscOptionsReal("-svd_trlanczos_scale","Scale parameter for matrix B","SVDTRLanczosSetScale",scale,&scale,&flg));
1730: if (flg) PetscCall(SVDTRLanczosSetScale(svd,scale));
1732: PetscOptionsHeadEnd();
1734: if (svd->OPb) {
1735: if (!lanczos->ksp) PetscCall(SVDTRLanczosGetKSP(svd,&lanczos->ksp));
1736: PetscCall(KSPSetFromOptions(lanczos->ksp));
1737: }
1738: PetscFunctionReturn(PETSC_SUCCESS);
1739: }
1741: static PetscErrorCode SVDTRLanczosSetOneSide_TRLanczos(SVD svd,PetscBool oneside)
1742: {
1743: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1745: PetscFunctionBegin;
1746: if (lanczos->oneside != oneside) {
1747: lanczos->oneside = oneside;
1748: svd->state = SVD_STATE_INITIAL;
1749: }
1750: PetscFunctionReturn(PETSC_SUCCESS);
1751: }
1753: /*@
1754: SVDTRLanczosSetOneSide - Indicate if the variant of the Lanczos method
1755: to be used is one-sided or two-sided.
1757: Logically Collective
1759: Input Parameters:
1760: + svd - singular value solver
1761: - oneside - boolean flag indicating if the method is one-sided or not
1763: Options Database Key:
1764: . -svd_trlanczos_oneside <boolean> - Indicates the boolean flag
1766: Notes:
1767: By default, a two-sided variant is selected, which is sometimes slightly
1768: more robust. However, the one-sided variant is faster because it avoids
1769: the orthogonalization associated to left singular vectors.
1771: One-sided orthogonalization is also available for the GSVD, in which case
1772: two orthogonalizations out of three are avoided.
1774: Level: advanced
1776: .seealso: SVDLanczosSetOneSide()
1777: @*/
1778: PetscErrorCode SVDTRLanczosSetOneSide(SVD svd,PetscBool oneside)
1779: {
1780: PetscFunctionBegin;
1783: PetscTryMethod(svd,"SVDTRLanczosSetOneSide_C",(SVD,PetscBool),(svd,oneside));
1784: PetscFunctionReturn(PETSC_SUCCESS);
1785: }
1787: static PetscErrorCode SVDTRLanczosGetOneSide_TRLanczos(SVD svd,PetscBool *oneside)
1788: {
1789: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1791: PetscFunctionBegin;
1792: *oneside = lanczos->oneside;
1793: PetscFunctionReturn(PETSC_SUCCESS);
1794: }
1796: /*@
1797: SVDTRLanczosGetOneSide - Gets if the variant of the Lanczos method
1798: to be used is one-sided or two-sided.
1800: Not Collective
1802: Input Parameters:
1803: . svd - singular value solver
1805: Output Parameters:
1806: . oneside - boolean flag indicating if the method is one-sided or not
1808: Level: advanced
1810: .seealso: SVDTRLanczosSetOneSide()
1811: @*/
1812: PetscErrorCode SVDTRLanczosGetOneSide(SVD svd,PetscBool *oneside)
1813: {
1814: PetscFunctionBegin;
1816: PetscAssertPointer(oneside,2);
1817: PetscUseMethod(svd,"SVDTRLanczosGetOneSide_C",(SVD,PetscBool*),(svd,oneside));
1818: PetscFunctionReturn(PETSC_SUCCESS);
1819: }
1821: static PetscErrorCode SVDTRLanczosSetGBidiag_TRLanczos(SVD svd,SVDTRLanczosGBidiag bidiag)
1822: {
1823: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1825: PetscFunctionBegin;
1826: switch (bidiag) {
1827: case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1828: case SVD_TRLANCZOS_GBIDIAG_UPPER:
1829: case SVD_TRLANCZOS_GBIDIAG_LOWER:
1830: if (lanczos->bidiag != bidiag) {
1831: lanczos->bidiag = bidiag;
1832: svd->state = SVD_STATE_INITIAL;
1833: }
1834: break;
1835: default:
1836: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid bidiagonalization choice");
1837: }
1838: PetscFunctionReturn(PETSC_SUCCESS);
1839: }
1841: /*@
1842: SVDTRLanczosSetGBidiag - Sets the bidiagonalization choice to use in
1843: the GSVD TRLanczos solver.
1845: Logically Collective
1847: Input Parameters:
1848: + svd - the singular value solver
1849: - bidiag - the bidiagonalization choice
1851: Options Database Key:
1852: . -svd_trlanczos_gbidiag - Sets the bidiagonalization choice (either 's' or 'juu'
1853: or 'jlu')
1855: Level: advanced
1857: .seealso: SVDTRLanczosGetGBidiag(), SVDTRLanczosGBidiag
1858: @*/
1859: PetscErrorCode SVDTRLanczosSetGBidiag(SVD svd,SVDTRLanczosGBidiag bidiag)
1860: {
1861: PetscFunctionBegin;
1864: PetscTryMethod(svd,"SVDTRLanczosSetGBidiag_C",(SVD,SVDTRLanczosGBidiag),(svd,bidiag));
1865: PetscFunctionReturn(PETSC_SUCCESS);
1866: }
1868: static PetscErrorCode SVDTRLanczosGetGBidiag_TRLanczos(SVD svd,SVDTRLanczosGBidiag *bidiag)
1869: {
1870: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1872: PetscFunctionBegin;
1873: *bidiag = lanczos->bidiag;
1874: PetscFunctionReturn(PETSC_SUCCESS);
1875: }
1877: /*@
1878: SVDTRLanczosGetGBidiag - Gets the bidiagonalization choice used in the GSVD
1879: TRLanczos solver.
1881: Not Collective
1883: Input Parameter:
1884: . svd - the singular value solver
1886: Output Parameter:
1887: . bidiag - the bidiagonalization choice
1889: Level: advanced
1891: .seealso: SVDTRLanczosSetGBidiag(), SVDTRLanczosGBidiag
1892: @*/
1893: PetscErrorCode SVDTRLanczosGetGBidiag(SVD svd,SVDTRLanczosGBidiag *bidiag)
1894: {
1895: PetscFunctionBegin;
1897: PetscAssertPointer(bidiag,2);
1898: PetscUseMethod(svd,"SVDTRLanczosGetGBidiag_C",(SVD,SVDTRLanczosGBidiag*),(svd,bidiag));
1899: PetscFunctionReturn(PETSC_SUCCESS);
1900: }
1902: static PetscErrorCode SVDTRLanczosSetKSP_TRLanczos(SVD svd,KSP ksp)
1903: {
1904: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1906: PetscFunctionBegin;
1907: PetscCall(PetscObjectReference((PetscObject)ksp));
1908: PetscCall(KSPDestroy(&ctx->ksp));
1909: ctx->ksp = ksp;
1910: svd->state = SVD_STATE_INITIAL;
1911: PetscFunctionReturn(PETSC_SUCCESS);
1912: }
1914: /*@
1915: SVDTRLanczosSetKSP - Associate a linear solver object (KSP) to the SVD solver.
1917: Collective
1919: Input Parameters:
1920: + svd - SVD solver
1921: - ksp - the linear solver object
1923: Note:
1924: Only used for the GSVD problem.
1926: Level: advanced
1928: .seealso: SVDTRLanczosGetKSP()
1929: @*/
1930: PetscErrorCode SVDTRLanczosSetKSP(SVD svd,KSP ksp)
1931: {
1932: PetscFunctionBegin;
1935: PetscCheckSameComm(svd,1,ksp,2);
1936: PetscTryMethod(svd,"SVDTRLanczosSetKSP_C",(SVD,KSP),(svd,ksp));
1937: PetscFunctionReturn(PETSC_SUCCESS);
1938: }
1940: static PetscErrorCode SVDTRLanczosGetKSP_TRLanczos(SVD svd,KSP *ksp)
1941: {
1942: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1943: PC pc;
1945: PetscFunctionBegin;
1946: if (!ctx->ksp) {
1947: /* Create linear solver */
1948: PetscCall(KSPCreate(PetscObjectComm((PetscObject)svd),&ctx->ksp));
1949: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)svd,1));
1950: PetscCall(KSPSetOptionsPrefix(ctx->ksp,((PetscObject)svd)->prefix));
1951: PetscCall(KSPAppendOptionsPrefix(ctx->ksp,"svd_trlanczos_"));
1952: PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)svd)->options));
1953: PetscCall(KSPSetType(ctx->ksp,KSPLSQR));
1954: PetscCall(KSPGetPC(ctx->ksp,&pc));
1955: PetscCall(PCSetType(pc,PCNONE));
1956: PetscCall(KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE));
1957: PetscCall(KSPSetTolerances(ctx->ksp,SlepcDefaultTol(svd->tol)/10.0,PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
1958: }
1959: *ksp = ctx->ksp;
1960: PetscFunctionReturn(PETSC_SUCCESS);
1961: }
1963: /*@
1964: SVDTRLanczosGetKSP - Retrieve the linear solver object (KSP) associated with
1965: the SVD solver.
1967: Collective
1969: Input Parameter:
1970: . svd - SVD solver
1972: Output Parameter:
1973: . ksp - the linear solver object
1975: Level: advanced
1977: .seealso: SVDTRLanczosSetKSP()
1978: @*/
1979: PetscErrorCode SVDTRLanczosGetKSP(SVD svd,KSP *ksp)
1980: {
1981: PetscFunctionBegin;
1983: PetscAssertPointer(ksp,2);
1984: PetscUseMethod(svd,"SVDTRLanczosGetKSP_C",(SVD,KSP*),(svd,ksp));
1985: PetscFunctionReturn(PETSC_SUCCESS);
1986: }
1988: static PetscErrorCode SVDTRLanczosSetRestart_TRLanczos(SVD svd,PetscReal keep)
1989: {
1990: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1992: PetscFunctionBegin;
1993: if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) ctx->keep = 0.5;
1994: else {
1995: 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);
1996: ctx->keep = keep;
1997: }
1998: PetscFunctionReturn(PETSC_SUCCESS);
1999: }
2001: /*@
2002: SVDTRLanczosSetRestart - Sets the restart parameter for the thick-restart
2003: Lanczos method, in particular the proportion of basis vectors that must be
2004: kept after restart.
2006: Logically Collective
2008: Input Parameters:
2009: + svd - the singular value solver
2010: - keep - the number of vectors to be kept at restart
2012: Options Database Key:
2013: . -svd_trlanczos_restart - Sets the restart parameter
2015: Notes:
2016: Allowed values are in the range [0.1,0.9]. The default is 0.5.
2018: Level: advanced
2020: .seealso: SVDTRLanczosGetRestart()
2021: @*/
2022: PetscErrorCode SVDTRLanczosSetRestart(SVD svd,PetscReal keep)
2023: {
2024: PetscFunctionBegin;
2027: PetscTryMethod(svd,"SVDTRLanczosSetRestart_C",(SVD,PetscReal),(svd,keep));
2028: PetscFunctionReturn(PETSC_SUCCESS);
2029: }
2031: static PetscErrorCode SVDTRLanczosGetRestart_TRLanczos(SVD svd,PetscReal *keep)
2032: {
2033: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
2035: PetscFunctionBegin;
2036: *keep = ctx->keep;
2037: PetscFunctionReturn(PETSC_SUCCESS);
2038: }
2040: /*@
2041: SVDTRLanczosGetRestart - Gets the restart parameter used in the thick-restart
2042: Lanczos method.
2044: Not Collective
2046: Input Parameter:
2047: . svd - the singular value solver
2049: Output Parameter:
2050: . keep - the restart parameter
2052: Level: advanced
2054: .seealso: SVDTRLanczosSetRestart()
2055: @*/
2056: PetscErrorCode SVDTRLanczosGetRestart(SVD svd,PetscReal *keep)
2057: {
2058: PetscFunctionBegin;
2060: PetscAssertPointer(keep,2);
2061: PetscUseMethod(svd,"SVDTRLanczosGetRestart_C",(SVD,PetscReal*),(svd,keep));
2062: PetscFunctionReturn(PETSC_SUCCESS);
2063: }
2065: static PetscErrorCode SVDTRLanczosSetLocking_TRLanczos(SVD svd,PetscBool lock)
2066: {
2067: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
2069: PetscFunctionBegin;
2070: ctx->lock = lock;
2071: PetscFunctionReturn(PETSC_SUCCESS);
2072: }
2074: /*@
2075: SVDTRLanczosSetLocking - Choose between locking and non-locking variants of
2076: the thick-restart Lanczos method.
2078: Logically Collective
2080: Input Parameters:
2081: + svd - the singular value solver
2082: - lock - true if the locking variant must be selected
2084: Options Database Key:
2085: . -svd_trlanczos_locking - Sets the locking flag
2087: Notes:
2088: The default is to lock converged singular triplets when the method restarts.
2089: This behaviour can be changed so that all directions are kept in the
2090: working subspace even if already converged to working accuracy (the
2091: non-locking variant).
2093: Level: advanced
2095: .seealso: SVDTRLanczosGetLocking()
2096: @*/
2097: PetscErrorCode SVDTRLanczosSetLocking(SVD svd,PetscBool lock)
2098: {
2099: PetscFunctionBegin;
2102: PetscTryMethod(svd,"SVDTRLanczosSetLocking_C",(SVD,PetscBool),(svd,lock));
2103: PetscFunctionReturn(PETSC_SUCCESS);
2104: }
2106: static PetscErrorCode SVDTRLanczosGetLocking_TRLanczos(SVD svd,PetscBool *lock)
2107: {
2108: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
2110: PetscFunctionBegin;
2111: *lock = ctx->lock;
2112: PetscFunctionReturn(PETSC_SUCCESS);
2113: }
2115: /*@
2116: SVDTRLanczosGetLocking - Gets the locking flag used in the thick-restart
2117: Lanczos method.
2119: Not Collective
2121: Input Parameter:
2122: . svd - the singular value solver
2124: Output Parameter:
2125: . lock - the locking flag
2127: Level: advanced
2129: .seealso: SVDTRLanczosSetLocking()
2130: @*/
2131: PetscErrorCode SVDTRLanczosGetLocking(SVD svd,PetscBool *lock)
2132: {
2133: PetscFunctionBegin;
2135: PetscAssertPointer(lock,2);
2136: PetscUseMethod(svd,"SVDTRLanczosGetLocking_C",(SVD,PetscBool*),(svd,lock));
2137: PetscFunctionReturn(PETSC_SUCCESS);
2138: }
2140: static PetscErrorCode SVDTRLanczosSetExplicitMatrix_TRLanczos(SVD svd,PetscBool explicitmat)
2141: {
2142: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2144: PetscFunctionBegin;
2145: if (lanczos->explicitmatrix != explicitmat) {
2146: lanczos->explicitmatrix = explicitmat;
2147: svd->state = SVD_STATE_INITIAL;
2148: }
2149: PetscFunctionReturn(PETSC_SUCCESS);
2150: }
2152: /*@
2153: SVDTRLanczosSetExplicitMatrix - Indicate if the matrix Z=[A;B] must
2154: be built explicitly.
2156: Logically Collective
2158: Input Parameters:
2159: + svd - singular value solver
2160: - explicitmat - Boolean flag indicating if Z=[A;B] is built explicitly
2162: Options Database Key:
2163: . -svd_trlanczos_explicitmatrix <boolean> - Indicates the boolean flag
2165: Notes:
2166: This option is relevant for the GSVD case only.
2167: Z is the coefficient matrix of the KSP solver used internally.
2169: Level: advanced
2171: .seealso: SVDTRLanczosGetExplicitMatrix()
2172: @*/
2173: PetscErrorCode SVDTRLanczosSetExplicitMatrix(SVD svd,PetscBool explicitmat)
2174: {
2175: PetscFunctionBegin;
2178: PetscTryMethod(svd,"SVDTRLanczosSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
2179: PetscFunctionReturn(PETSC_SUCCESS);
2180: }
2182: static PetscErrorCode SVDTRLanczosGetExplicitMatrix_TRLanczos(SVD svd,PetscBool *explicitmat)
2183: {
2184: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2186: PetscFunctionBegin;
2187: *explicitmat = lanczos->explicitmatrix;
2188: PetscFunctionReturn(PETSC_SUCCESS);
2189: }
2191: /*@
2192: SVDTRLanczosGetExplicitMatrix - Returns the flag indicating if Z=[A;B] is built explicitly.
2194: Not Collective
2196: Input Parameter:
2197: . svd - singular value solver
2199: Output Parameter:
2200: . explicitmat - the mode flag
2202: Level: advanced
2204: .seealso: SVDTRLanczosSetExplicitMatrix()
2205: @*/
2206: PetscErrorCode SVDTRLanczosGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
2207: {
2208: PetscFunctionBegin;
2210: PetscAssertPointer(explicitmat,2);
2211: PetscUseMethod(svd,"SVDTRLanczosGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
2212: PetscFunctionReturn(PETSC_SUCCESS);
2213: }
2215: static PetscErrorCode SVDTRLanczosSetScale_TRLanczos(SVD svd,PetscReal scale)
2216: {
2217: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
2219: PetscFunctionBegin;
2220: if (scale<0) {
2221: ctx->scalef = 1.0;
2222: ctx->scaleth = -scale;
2223: } else {
2224: ctx->scalef = scale;
2225: ctx->scaleth = 0.0;
2226: }
2227: PetscFunctionReturn(PETSC_SUCCESS);
2228: }
2230: /*@
2231: SVDTRLanczosSetScale - Sets the scale parameter for the GSVD.
2233: Logically Collective
2235: Input Parameters:
2236: + svd - singular value solver
2237: - scale - scale parameter
2239: Options Database Key:
2240: . -svd_trlanczos_scale <real> - scale factor/threshold
2242: Notes:
2243: This parameter is relevant for the GSVD case only. If the parameter is
2244: positive, it indicates the scale factor for B in matrix Z=[A;B]. If
2245: negative, its absolute value is the threshold for automatic scaling.
2246: In automatic scaling, whenever the largest approximate generalized singular
2247: value (or the inverse of the smallest value, if SVD_SMALLEST is used)
2248: exceeds the threshold, the computation is restarted with matrix B
2249: scaled by that value.
2251: Level: advanced
2253: .seealso: SVDTRLanczosGetScale()
2254: @*/
2255: PetscErrorCode SVDTRLanczosSetScale(SVD svd,PetscReal scale)
2256: {
2257: PetscFunctionBegin;
2260: PetscTryMethod(svd,"SVDTRLanczosSetScale_C",(SVD,PetscReal),(svd,scale));
2261: PetscFunctionReturn(PETSC_SUCCESS);
2262: }
2264: static PetscErrorCode SVDTRLanczosGetScale_TRLanczos(SVD svd,PetscReal *scale)
2265: {
2266: SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
2268: PetscFunctionBegin;
2269: if (ctx->scaleth==0) *scale = ctx->scalef;
2270: else *scale = -ctx->scaleth;
2271: PetscFunctionReturn(PETSC_SUCCESS);
2272: }
2274: /*@
2275: SVDTRLanczosGetScale - Gets the scale parameter for the GSVD.
2277: Not Collective
2279: Input Parameter:
2280: . svd - the singular value solver
2282: Output Parameter:
2283: . scale - the scale parameter
2285: Notes:
2286: This parameter is relevant for the GSVD case only. If the parameter is
2287: positive, it indicates the scale factor for B in matrix Z=[A;B]. If
2288: negative, its absolute value is the threshold for automatic scaling.
2290: Level: advanced
2292: .seealso: SVDTRLanczosSetScale()
2293: @*/
2294: PetscErrorCode SVDTRLanczosGetScale(SVD svd,PetscReal *scale)
2295: {
2296: PetscFunctionBegin;
2298: PetscAssertPointer(scale,2);
2299: PetscUseMethod(svd,"SVDTRLanczosGetScale_C",(SVD,PetscReal*),(svd,scale));
2300: PetscFunctionReturn(PETSC_SUCCESS);
2301: }
2303: static PetscErrorCode SVDReset_TRLanczos(SVD svd)
2304: {
2305: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2307: PetscFunctionBegin;
2308: if (svd->isgeneralized || (!svd->problem_type && svd->OPb)) {
2309: PetscCall(KSPReset(lanczos->ksp));
2310: PetscCall(MatDestroy(&lanczos->Z));
2311: }
2312: PetscFunctionReturn(PETSC_SUCCESS);
2313: }
2315: static PetscErrorCode SVDDestroy_TRLanczos(SVD svd)
2316: {
2317: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2319: PetscFunctionBegin;
2320: if (svd->isgeneralized || (!svd->problem_type && svd->OPb)) PetscCall(KSPDestroy(&lanczos->ksp));
2321: PetscCall(PetscFree(svd->data));
2322: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",NULL));
2323: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",NULL));
2324: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetGBidiag_C",NULL));
2325: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetGBidiag_C",NULL));
2326: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetKSP_C",NULL));
2327: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetKSP_C",NULL));
2328: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetRestart_C",NULL));
2329: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetRestart_C",NULL));
2330: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetLocking_C",NULL));
2331: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetLocking_C",NULL));
2332: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetExplicitMatrix_C",NULL));
2333: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetExplicitMatrix_C",NULL));
2334: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetScale_C",NULL));
2335: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetScale_C",NULL));
2336: PetscFunctionReturn(PETSC_SUCCESS);
2337: }
2339: static PetscErrorCode SVDView_TRLanczos(SVD svd,PetscViewer viewer)
2340: {
2341: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2342: PetscBool isascii;
2344: PetscFunctionBegin;
2345: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
2346: if (isascii) {
2347: PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*lanczos->keep)));
2348: PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",lanczos->lock?"":"non-"));
2349: if (svd->isgeneralized) {
2350: const char *bidiag="";
2352: switch (lanczos->bidiag) {
2353: case SVD_TRLANCZOS_GBIDIAG_SINGLE: bidiag = "single"; break;
2354: case SVD_TRLANCZOS_GBIDIAG_UPPER: bidiag = "joint upper-upper"; break;
2355: case SVD_TRLANCZOS_GBIDIAG_LOWER: bidiag = "joint lower-upper"; break;
2356: }
2357: PetscCall(PetscViewerASCIIPrintf(viewer," bidiagonalization choice: %s\n",bidiag));
2358: PetscCall(PetscViewerASCIIPrintf(viewer," %s matrix\n",lanczos->explicitmatrix?"explicit":"implicit"));
2359: if (lanczos->scaleth==0) PetscCall(PetscViewerASCIIPrintf(viewer," scale factor for matrix B: %g\n",(double)lanczos->scalef));
2360: else PetscCall(PetscViewerASCIIPrintf(viewer," automatic scaling for matrix B with threshold: %g\n",(double)lanczos->scaleth));
2361: if (!lanczos->ksp) PetscCall(SVDTRLanczosGetKSP(svd,&lanczos->ksp));
2362: PetscCall(PetscViewerASCIIPushTab(viewer));
2363: PetscCall(KSPView(lanczos->ksp,viewer));
2364: PetscCall(PetscViewerASCIIPopTab(viewer));
2365: } else PetscCall(PetscViewerASCIIPrintf(viewer," %s-sided reorthogonalization\n",lanczos->oneside? "one": "two"));
2366: }
2367: PetscFunctionReturn(PETSC_SUCCESS);
2368: }
2370: static PetscErrorCode SVDSetDSType_TRLanczos(SVD svd)
2371: {
2372: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2373: DSType dstype;
2375: PetscFunctionBegin;
2376: dstype = svd->ishyperbolic? DSHSVD: DSSVD;
2377: if (svd->OPb && (lanczos->bidiag==SVD_TRLANCZOS_GBIDIAG_UPPER || lanczos->bidiag==SVD_TRLANCZOS_GBIDIAG_LOWER)) dstype = DSGSVD;
2378: PetscCall(DSSetType(svd->ds,dstype));
2379: PetscFunctionReturn(PETSC_SUCCESS);
2380: }
2382: SLEPC_EXTERN PetscErrorCode SVDCreate_TRLanczos(SVD svd)
2383: {
2384: SVD_TRLANCZOS *ctx;
2386: PetscFunctionBegin;
2387: PetscCall(PetscNew(&ctx));
2388: svd->data = (void*)ctx;
2390: ctx->lock = PETSC_TRUE;
2391: ctx->bidiag = SVD_TRLANCZOS_GBIDIAG_LOWER;
2392: ctx->scalef = 1.0;
2393: ctx->scaleth = 0.0;
2395: svd->ops->setup = SVDSetUp_TRLanczos;
2396: svd->ops->solve = SVDSolve_TRLanczos;
2397: svd->ops->solveg = SVDSolve_TRLanczos_GSVD;
2398: svd->ops->solveh = SVDSolve_TRLanczos_HSVD;
2399: svd->ops->destroy = SVDDestroy_TRLanczos;
2400: svd->ops->reset = SVDReset_TRLanczos;
2401: svd->ops->setfromoptions = SVDSetFromOptions_TRLanczos;
2402: svd->ops->view = SVDView_TRLanczos;
2403: svd->ops->setdstype = SVDSetDSType_TRLanczos;
2404: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",SVDTRLanczosSetOneSide_TRLanczos));
2405: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",SVDTRLanczosGetOneSide_TRLanczos));
2406: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetGBidiag_C",SVDTRLanczosSetGBidiag_TRLanczos));
2407: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetGBidiag_C",SVDTRLanczosGetGBidiag_TRLanczos));
2408: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetKSP_C",SVDTRLanczosSetKSP_TRLanczos));
2409: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetKSP_C",SVDTRLanczosGetKSP_TRLanczos));
2410: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetRestart_C",SVDTRLanczosSetRestart_TRLanczos));
2411: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetRestart_C",SVDTRLanczosGetRestart_TRLanczos));
2412: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetLocking_C",SVDTRLanczosSetLocking_TRLanczos));
2413: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetLocking_C",SVDTRLanczosGetLocking_TRLanczos));
2414: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetExplicitMatrix_C",SVDTRLanczosSetExplicitMatrix_TRLanczos));
2415: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetExplicitMatrix_C",SVDTRLanczosGetExplicitMatrix_TRLanczos));
2416: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetScale_C",SVDTRLanczosSetScale_TRLanczos));
2417: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetScale_C",SVDTRLanczosGetScale_TRLanczos));
2418: PetscFunctionReturn(PETSC_SUCCESS);
2419: }