Actual source code: trlanczos.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepcsvd.h" I*/
33: static PetscBool cited = PETSC_FALSE;
34: static const char citation[] =
35: "@Article{slepc-svd,\n"
36: " author = \"V. Hern{\\'a}ndez and J. E. Rom{\\'a}n and A. Tom{\\'a}s\",\n"
37: " title = \"A robust and efficient parallel {SVD} solver based on restarted {Lanczos} bidiagonalization\",\n"
38: " journal = \"Electron. Trans. Numer. Anal.\",\n"
39: " volume = \"31\",\n"
40: " pages = \"68--85\",\n"
41: " year = \"2008\"\n"
42: "}\n";
44: typedef struct {
45: PetscBool oneside;
46: } SVD_TRLANCZOS;
48: PetscErrorCode SVDSetUp_TRLanczos(SVD svd)
49: {
51: PetscInt N;
54: SVDMatGetSize(svd,NULL,&N);
55: SVDSetDimensions_Default(svd);
56: if (svd->ncv>svd->nsv+svd->mpd) SETERRQ(PetscObjectComm((PetscObject)svd),1,"The value of ncv must not be larger than nev+mpd");
57: if (!svd->max_it) svd->max_it = PetscMax(N/svd->ncv,100);
58: svd->leftbasis = PETSC_TRUE;
59: SVDAllocateSolution(svd,1);
60: DSSetType(svd->ds,DSSVD);
61: DSSetCompact(svd->ds,PETSC_TRUE);
62: DSAllocate(svd->ds,svd->ncv);
63: return(0);
64: }
66: static PetscErrorCode SVDOneSideTRLanczosMGS(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,BV U,PetscInt nconv,PetscInt l,PetscInt n,PetscScalar* work)
67: {
69: PetscReal a,b;
70: PetscInt i,k=nconv+l;
71: Vec ui,ui1,vi;
74: BVGetColumn(V,k,&vi);
75: BVGetColumn(U,k,&ui);
76: SVDMatMult(svd,PETSC_FALSE,vi,ui);
77: BVRestoreColumn(V,k,&vi);
78: BVRestoreColumn(U,k,&ui);
79: if (l>0) {
80: for (i=0;i<l;i++) work[i]=beta[i+nconv];
81: BVMultColumn(U,-1.0,1.0,k,work);
82: }
83: BVNormColumn(U,k,NORM_2,&a);
84: BVScaleColumn(U,k,1.0/a);
85: alpha[k] = a;
87: for (i=k+1;i<n;i++) {
88: BVGetColumn(V,i,&vi);
89: BVGetColumn(U,i-1,&ui1);
90: SVDMatMult(svd,PETSC_TRUE,ui1,vi);
91: BVRestoreColumn(V,i,&vi);
92: BVRestoreColumn(U,i-1,&ui1);
93: BVOrthogonalizeColumn(V,i,NULL,&b,NULL);
94: BVScaleColumn(V,i,1.0/b);
95: beta[i-1] = b;
97: BVGetColumn(V,i,&vi);
98: BVGetColumn(U,i,&ui);
99: SVDMatMult(svd,PETSC_FALSE,vi,ui);
100: BVRestoreColumn(V,i,&vi);
101: BVGetColumn(U,i-1,&ui1);
102: VecAXPY(ui,-b,ui1);
103: BVRestoreColumn(U,i-1,&ui1);
104: BVRestoreColumn(U,i,&ui);
105: BVNormColumn(U,i,NORM_2,&a);
106: BVScaleColumn(U,i,1.0/a);
107: alpha[i] = a;
108: }
110: BVGetColumn(V,n,&vi);
111: BVGetColumn(U,n-1,&ui1);
112: SVDMatMult(svd,PETSC_TRUE,ui1,vi);
113: BVRestoreColumn(V,n,&vi);
114: BVRestoreColumn(U,n-1,&ui1);
115: BVOrthogonalizeColumn(V,n,NULL,&b,NULL);
116: beta[n-1] = b;
117: return(0);
118: }
120: /*
121: Custom CGS orthogonalization, preprocess after first orthogonalization
122: */
123: static PetscErrorCode SVDOrthogonalizeCGS(BV V,PetscInt i,PetscScalar* h,PetscReal a,BVOrthogRefineType refine,PetscReal eta,PetscReal *norm)
124: {
126: PetscReal sum,onorm;
127: PetscScalar dot;
128: PetscInt j;
131: switch (refine) {
132: case BV_ORTHOG_REFINE_NEVER:
133: BVNormColumn(V,i,NORM_2,norm);
134: break;
135: case BV_ORTHOG_REFINE_ALWAYS:
136: BVSetActiveColumns(V,0,i);
137: BVDotColumn(V,i,h);
138: BVMultColumn(V,-1.0,1.0,i,h);
139: BVNormColumn(V,i,NORM_2,norm);
140: break;
141: case BV_ORTHOG_REFINE_IFNEEDED:
142: dot = h[i];
143: onorm = PetscSqrtReal(PetscRealPart(dot)) / a;
144: sum = 0.0;
145: for (j=0;j<i;j++) {
146: sum += PetscRealPart(h[j] * PetscConj(h[j]));
147: }
148: *norm = PetscRealPart(dot)/(a*a) - sum;
149: if (*norm>0.0) *norm = PetscSqrtReal(*norm);
150: else {
151: BVNormColumn(V,i,NORM_2,norm);
152: }
153: if (*norm < eta*onorm) {
154: BVSetActiveColumns(V,0,i);
155: BVDotColumn(V,i,h);
156: BVMultColumn(V,-1.0,1.0,i,h);
157: BVNormColumn(V,i,NORM_2,norm);
158: }
159: break;
160: }
161: return(0);
162: }
164: static PetscErrorCode SVDOneSideTRLanczosCGS(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,BV U,PetscInt nconv,PetscInt l,PetscInt n,PetscScalar* work)
165: {
166: PetscErrorCode ierr;
167: PetscReal a,b,eta;
168: PetscInt i,j,k=nconv+l;
169: Vec ui,ui1,vi;
170: BVOrthogRefineType refine;
173: BVGetColumn(V,k,&vi);
174: BVGetColumn(U,k,&ui);
175: SVDMatMult(svd,PETSC_FALSE,vi,ui);
176: BVRestoreColumn(V,k,&vi);
177: BVRestoreColumn(U,k,&ui);
178: if (l>0) {
179: for (i=0;i<l;i++) work[i]=beta[i+nconv];
180: BVMultColumn(U,-1.0,1.0,k,work);
181: }
182: BVGetOrthogonalization(V,NULL,&refine,&eta,NULL);
184: for (i=k+1;i<n;i++) {
185: BVGetColumn(V,i,&vi);
186: BVGetColumn(U,i-1,&ui1);
187: SVDMatMult(svd,PETSC_TRUE,ui1,vi);
188: BVRestoreColumn(V,i,&vi);
189: BVRestoreColumn(U,i-1,&ui1);
190: BVNormColumnBegin(U,i-1,NORM_2,&a);
191: if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
192: BVSetActiveColumns(V,0,i+1);
193: BVGetColumn(V,i,&vi);
194: BVDotVecBegin(V,vi,work);
195: } else {
196: BVSetActiveColumns(V,0,i);
197: BVDotColumnBegin(V,i,work);
198: }
199: BVNormColumnEnd(U,i-1,NORM_2,&a);
200: if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
201: BVDotVecEnd(V,vi,work);
202: BVRestoreColumn(V,i,&vi);
203: BVSetActiveColumns(V,0,i);
204: } else {
205: BVDotColumnEnd(V,i,work);
206: }
208: BVScaleColumn(U,i-1,1.0/a);
209: for (j=0;j<i;j++) work[j] = work[j] / a;
210: BVMultColumn(V,-1.0,1.0/a,i,work);
211: SVDOrthogonalizeCGS(V,i,work,a,refine,eta,&b);
212: BVScaleColumn(V,i,1.0/b);
213: if (PetscAbsReal(b)<10*PETSC_MACHINE_EPSILON) SETERRQ(PetscObjectComm((PetscObject)svd),1,"Recurrence generated a zero vector; use a two-sided variant");
215: BVGetColumn(V,i,&vi);
216: BVGetColumn(U,i,&ui);
217: BVGetColumn(U,i-1,&ui1);
218: SVDMatMult(svd,PETSC_FALSE,vi,ui);
219: VecAXPY(ui,-b,ui1);
220: BVRestoreColumn(V,i,&vi);
221: BVRestoreColumn(U,i,&ui);
222: BVRestoreColumn(U,i-1,&ui1);
224: alpha[i-1] = a;
225: beta[i-1] = b;
226: }
228: BVGetColumn(V,n,&vi);
229: BVGetColumn(U,n-1,&ui1);
230: SVDMatMult(svd,PETSC_TRUE,ui1,vi);
231: BVRestoreColumn(V,n,&vi);
232: BVRestoreColumn(U,n-1,&ui1);
234: BVNormColumnBegin(svd->U,n-1,NORM_2,&a);
235: if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
236: BVSetActiveColumns(V,0,n+1);
237: BVGetColumn(V,n,&vi);
238: BVDotVecBegin(V,vi,work);
239: } else {
240: BVSetActiveColumns(V,0,n);
241: BVDotColumnBegin(V,n,work);
242: }
243: BVNormColumnEnd(svd->U,n-1,NORM_2,&a);
244: if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
245: BVDotVecEnd(V,vi,work);
246: BVRestoreColumn(V,n,&vi);
247: } else {
248: BVDotColumnEnd(V,n,work);
249: }
251: BVScaleColumn(U,n-1,1.0/a);
252: for (j=0;j<n;j++) work[j] = work[j] / a;
253: BVMultColumn(V,-1.0,1.0/a,n,work);
254: SVDOrthogonalizeCGS(V,n,work,a,refine,eta,&b);
255: BVSetActiveColumns(V,nconv,n);
256: alpha[n-1] = a;
257: beta[n-1] = b;
258: return(0);
259: }
261: PetscErrorCode SVDSolve_TRLanczos(SVD svd)
262: {
264: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
265: PetscReal *alpha,*beta,lastbeta,resnorm;
266: PetscScalar *Q,*swork=NULL,*w;
267: PetscInt i,k,l,nv,ld;
268: Mat U,VT;
269: PetscBool conv;
270: BVOrthogType orthog;
273: PetscCitationsRegister(citation,&cited);
274: /* allocate working space */
275: DSGetLeadingDimension(svd->ds,&ld);
276: BVGetOrthogonalization(svd->V,&orthog,NULL,NULL,NULL);
277: PetscMalloc1(ld,&w);
278: if (lanczos->oneside) {
279: PetscMalloc1(svd->ncv+1,&swork);
280: }
282: /* normalize start vector */
283: if (!svd->nini) {
284: BVSetRandomColumn(svd->V,0);
285: BVOrthonormalizeColumn(svd->V,0,PETSC_TRUE,NULL,NULL);
286: }
288: l = 0;
289: while (svd->reason == SVD_CONVERGED_ITERATING) {
290: svd->its++;
292: /* inner loop */
293: nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
294: BVSetActiveColumns(svd->V,svd->nconv,nv);
295: BVSetActiveColumns(svd->U,svd->nconv,nv);
296: DSGetArrayReal(svd->ds,DS_MAT_T,&alpha);
297: beta = alpha + ld;
298: if (lanczos->oneside) {
299: if (orthog == BV_ORTHOG_MGS) {
300: SVDOneSideTRLanczosMGS(svd,alpha,beta,svd->V,svd->U,svd->nconv,l,nv,swork);
301: } else {
302: SVDOneSideTRLanczosCGS(svd,alpha,beta,svd->V,svd->U,svd->nconv,l,nv,swork);
303: }
304: } else {
305: SVDTwoSideLanczos(svd,alpha,beta,svd->V,svd->U,svd->nconv+l,nv);
306: }
307: lastbeta = beta[nv-1];
308: DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha);
309: BVScaleColumn(svd->V,nv,1.0/lastbeta);
311: /* compute SVD of general matrix */
312: DSSetDimensions(svd->ds,nv,nv,svd->nconv,svd->nconv+l);
313: if (l==0) {
314: DSSetState(svd->ds,DS_STATE_INTERMEDIATE);
315: } else {
316: DSSetState(svd->ds,DS_STATE_RAW);
317: }
318: DSSolve(svd->ds,w,NULL);
319: DSSort(svd->ds,w,NULL,NULL,NULL,NULL);
320: DSSynchronize(svd->ds,w,NULL);
322: /* compute error estimates */
323: k = 0;
324: conv = PETSC_TRUE;
325: DSGetArray(svd->ds,DS_MAT_U,&Q);
326: DSGetArrayReal(svd->ds,DS_MAT_T,&alpha);
327: beta = alpha + ld;
328: for (i=svd->nconv;i<nv;i++) {
329: svd->sigma[i] = PetscRealPart(w[i]);
330: beta[i] = PetscRealPart(Q[nv-1+i*ld])*lastbeta;
331: resnorm = PetscAbsReal(beta[i]);
332: (*svd->converged)(svd,svd->sigma[i],resnorm,&svd->errest[i],svd->convergedctx);
333: if (conv) {
334: if (svd->errest[i] < svd->tol) k++;
335: else conv = PETSC_FALSE;
336: }
337: }
338: DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha);
339: DSRestoreArray(svd->ds,DS_MAT_U,&Q);
341: /* check convergence and update l */
342: (*svd->stopping)(svd,svd->its,svd->max_it,svd->nconv+k,svd->nsv,&svd->reason,svd->stoppingctx);
343: if (svd->reason != SVD_CONVERGED_ITERATING) l = 0;
344: else l = PetscMax((nv-svd->nconv-k)/2,0);
346: /* compute converged singular vectors and restart vectors */
347: DSGetMat(svd->ds,DS_MAT_VT,&VT);
348: BVMultInPlaceTranspose(svd->V,VT,svd->nconv,svd->nconv+k+l);
349: MatDestroy(&VT);
350: DSGetMat(svd->ds,DS_MAT_U,&U);
351: BVMultInPlace(svd->U,U,svd->nconv,svd->nconv+k+l);
352: MatDestroy(&U);
354: /* copy the last vector to be the next initial vector */
355: if (svd->reason == SVD_CONVERGED_ITERATING) {
356: BVCopyColumn(svd->V,nv,svd->nconv+k+l);
357: }
359: svd->nconv += k;
360: SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv);
361: }
363: /* orthonormalize U columns in one side method */
364: if (lanczos->oneside) {
365: for (i=0;i<svd->nconv;i++) {
366: BVOrthonormalizeColumn(svd->U,i,PETSC_FALSE,NULL,NULL);
367: }
368: }
370: /* free working space */
371: PetscFree(w);
372: if (swork) { PetscFree(swork); }
373: return(0);
374: }
376: PetscErrorCode SVDSetFromOptions_TRLanczos(PetscOptionItems *PetscOptionsObject,SVD svd)
377: {
379: PetscBool set,val;
380: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
383: PetscOptionsHead(PetscOptionsObject,"SVD TRLanczos Options");
385: PetscOptionsBool("-svd_trlanczos_oneside","Use one-side reorthogonalization","SVDTRLanczosSetOneSide",lanczos->oneside,&val,&set);
386: if (set) { SVDTRLanczosSetOneSide(svd,val); }
388: PetscOptionsTail();
389: return(0);
390: }
392: static PetscErrorCode SVDTRLanczosSetOneSide_TRLanczos(SVD svd,PetscBool oneside)
393: {
394: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
397: if (lanczos->oneside != oneside) {
398: lanczos->oneside = oneside;
399: svd->state = SVD_STATE_INITIAL;
400: }
401: return(0);
402: }
404: /*@
405: SVDTRLanczosSetOneSide - Indicate if the variant of the Lanczos method
406: to be used is one-sided or two-sided.
408: Logically Collective on SVD
410: Input Parameters:
411: + svd - singular value solver
412: - oneside - boolean flag indicating if the method is one-sided or not
414: Options Database Key:
415: . -svd_trlanczos_oneside <boolean> - Indicates the boolean flag
417: Note:
418: By default, a two-sided variant is selected, which is sometimes slightly
419: more robust. However, the one-sided variant is faster because it avoids
420: the orthogonalization associated to left singular vectors.
422: Level: advanced
424: .seealso: SVDLanczosSetOneSide()
425: @*/
426: PetscErrorCode SVDTRLanczosSetOneSide(SVD svd,PetscBool oneside)
427: {
433: PetscTryMethod(svd,"SVDTRLanczosSetOneSide_C",(SVD,PetscBool),(svd,oneside));
434: return(0);
435: }
437: static PetscErrorCode SVDTRLanczosGetOneSide_TRLanczos(SVD svd,PetscBool *oneside)
438: {
439: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
442: *oneside = lanczos->oneside;
443: return(0);
444: }
446: /*@
447: SVDTRLanczosGetOneSide - Gets if the variant of the Lanczos method
448: to be used is one-sided or two-sided.
450: Not Collective
452: Input Parameters:
453: . svd - singular value solver
455: Output Parameters:
456: . oneside - boolean flag indicating if the method is one-sided or not
458: Level: advanced
460: .seealso: SVDTRLanczosSetOneSide()
461: @*/
462: PetscErrorCode SVDTRLanczosGetOneSide(SVD svd,PetscBool *oneside)
463: {
469: PetscUseMethod(svd,"SVDTRLanczosGetOneSide_C",(SVD,PetscBool*),(svd,oneside));
470: return(0);
471: }
473: PetscErrorCode SVDDestroy_TRLanczos(SVD svd)
474: {
478: PetscFree(svd->data);
479: PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",NULL);
480: PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",NULL);
481: return(0);
482: }
484: PetscErrorCode SVDView_TRLanczos(SVD svd,PetscViewer viewer)
485: {
487: SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
488: PetscBool isascii;
491: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
492: if (isascii) {
493: PetscViewerASCIIPrintf(viewer," %s-sided reorthogonalization\n",lanczos->oneside? "one": "two");
494: }
495: return(0);
496: }
498: SLEPC_EXTERN PetscErrorCode SVDCreate_TRLanczos(SVD svd)
499: {
501: SVD_TRLANCZOS *ctx;
504: PetscNewLog(svd,&ctx);
505: svd->data = (void*)ctx;
507: svd->ops->setup = SVDSetUp_TRLanczos;
508: svd->ops->solve = SVDSolve_TRLanczos;
509: svd->ops->destroy = SVDDestroy_TRLanczos;
510: svd->ops->setfromoptions = SVDSetFromOptions_TRLanczos;
511: svd->ops->view = SVDView_TRLanczos;
512: PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",SVDTRLanczosSetOneSide_TRLanczos);
513: PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",SVDTRLanczosGetOneSide_TRLanczos);
514: return(0);
515: }