Actual source code: cyclic.c
slepc-3.18.0 2022-10-01
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: "cyclic"
13: Method: Uses a Hermitian eigensolver for H(A) = [ 0 A ; A^T 0 ]
14: */
16: #include <slepc/private/svdimpl.h>
17: #include <slepc/private/epsimpl.h>
18: #include "cyclic.h"
20: static PetscErrorCode MatMult_Cyclic(Mat B,Vec x,Vec y)
21: {
22: SVD_CYCLIC_SHELL *ctx;
23: const PetscScalar *px;
24: PetscScalar *py;
25: PetscInt m;
27: MatShellGetContext(B,&ctx);
28: MatGetLocalSize(ctx->A,&m,NULL);
29: VecGetArrayRead(x,&px);
30: VecGetArrayWrite(y,&py);
31: VecPlaceArray(ctx->x1,px);
32: VecPlaceArray(ctx->x2,px+m);
33: VecPlaceArray(ctx->y1,py);
34: VecPlaceArray(ctx->y2,py+m);
35: MatMult(ctx->A,ctx->x2,ctx->y1);
36: MatMult(ctx->AT,ctx->x1,ctx->y2);
37: VecResetArray(ctx->x1);
38: VecResetArray(ctx->x2);
39: VecResetArray(ctx->y1);
40: VecResetArray(ctx->y2);
41: VecRestoreArrayRead(x,&px);
42: VecRestoreArrayWrite(y,&py);
43: return 0;
44: }
46: static PetscErrorCode MatGetDiagonal_Cyclic(Mat B,Vec diag)
47: {
48: VecSet(diag,0.0);
49: return 0;
50: }
52: static PetscErrorCode MatDestroy_Cyclic(Mat B)
53: {
54: SVD_CYCLIC_SHELL *ctx;
56: MatShellGetContext(B,&ctx);
57: VecDestroy(&ctx->x1);
58: VecDestroy(&ctx->x2);
59: VecDestroy(&ctx->y1);
60: VecDestroy(&ctx->y2);
61: PetscFree(ctx);
62: return 0;
63: }
65: /*
66: Builds cyclic matrix C = | 0 A |
67: | AT 0 |
68: */
69: static PetscErrorCode SVDCyclicGetCyclicMat(SVD svd,Mat A,Mat AT,Mat *C)
70: {
71: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
72: SVD_CYCLIC_SHELL *ctx;
73: PetscInt i,M,N,m,n,Istart,Iend;
74: VecType vtype;
75: Mat Zm,Zn;
76: #if defined(PETSC_HAVE_CUDA)
77: PetscBool cuda;
78: #endif
80: MatGetSize(A,&M,&N);
81: MatGetLocalSize(A,&m,&n);
83: if (cyclic->explicitmatrix) {
85: MatCreate(PetscObjectComm((PetscObject)svd),&Zm);
86: MatSetSizes(Zm,m,m,M,M);
87: MatSetFromOptions(Zm);
88: MatSetUp(Zm);
89: MatGetOwnershipRange(Zm,&Istart,&Iend);
90: for (i=Istart;i<Iend;i++) MatSetValue(Zm,i,i,0.0,INSERT_VALUES);
91: MatAssemblyBegin(Zm,MAT_FINAL_ASSEMBLY);
92: MatAssemblyEnd(Zm,MAT_FINAL_ASSEMBLY);
93: MatCreate(PetscObjectComm((PetscObject)svd),&Zn);
94: MatSetSizes(Zn,n,n,N,N);
95: MatSetFromOptions(Zn);
96: MatSetUp(Zn);
97: MatGetOwnershipRange(Zn,&Istart,&Iend);
98: for (i=Istart;i<Iend;i++) MatSetValue(Zn,i,i,0.0,INSERT_VALUES);
99: MatAssemblyBegin(Zn,MAT_FINAL_ASSEMBLY);
100: MatAssemblyEnd(Zn,MAT_FINAL_ASSEMBLY);
101: MatCreateTile(1.0,Zm,1.0,A,1.0,AT,1.0,Zn,C);
102: MatDestroy(&Zm);
103: MatDestroy(&Zn);
104: } else {
105: PetscNew(&ctx);
106: ctx->A = A;
107: ctx->AT = AT;
108: ctx->swapped = svd->swapped;
109: MatCreateVecsEmpty(A,&ctx->x2,&ctx->x1);
110: MatCreateVecsEmpty(A,&ctx->y2,&ctx->y1);
111: MatCreateShell(PetscObjectComm((PetscObject)svd),m+n,m+n,M+N,M+N,ctx,C);
112: MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Cyclic);
113: MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_Cyclic);
114: #if defined(PETSC_HAVE_CUDA)
115: PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&cuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,"");
116: if (cuda) MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic_CUDA);
117: else
118: #endif
119: MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic);
120: MatGetVecType(A,&vtype);
121: MatSetVecType(*C,vtype);
122: }
123: return 0;
124: }
126: static PetscErrorCode MatMult_ECross(Mat B,Vec x,Vec y)
127: {
128: SVD_CYCLIC_SHELL *ctx;
129: const PetscScalar *px;
130: PetscScalar *py;
131: PetscInt mn,m,n;
133: MatShellGetContext(B,&ctx);
134: MatGetLocalSize(ctx->A,NULL,&n);
135: VecGetLocalSize(y,&mn);
136: m = mn-n;
137: VecGetArrayRead(x,&px);
138: VecGetArrayWrite(y,&py);
139: VecPlaceArray(ctx->x1,px);
140: VecPlaceArray(ctx->x2,px+m);
141: VecPlaceArray(ctx->y1,py);
142: VecPlaceArray(ctx->y2,py+m);
143: VecCopy(ctx->x1,ctx->y1);
144: MatMult(ctx->A,ctx->x2,ctx->w);
145: MatMult(ctx->AT,ctx->w,ctx->y2);
146: VecResetArray(ctx->x1);
147: VecResetArray(ctx->x2);
148: VecResetArray(ctx->y1);
149: VecResetArray(ctx->y2);
150: VecRestoreArrayRead(x,&px);
151: VecRestoreArrayWrite(y,&py);
152: return 0;
153: }
155: static PetscErrorCode MatGetDiagonal_ECross(Mat B,Vec d)
156: {
157: SVD_CYCLIC_SHELL *ctx;
158: PetscScalar *pd;
159: PetscMPIInt len;
160: PetscInt mn,m,n,N,i,j,start,end,ncols;
161: PetscScalar *work1,*work2,*diag;
162: const PetscInt *cols;
163: const PetscScalar *vals;
165: MatShellGetContext(B,&ctx);
166: MatGetLocalSize(ctx->A,NULL,&n);
167: VecGetLocalSize(d,&mn);
168: m = mn-n;
169: VecGetArrayWrite(d,&pd);
170: VecPlaceArray(ctx->y1,pd);
171: VecSet(ctx->y1,1.0);
172: VecResetArray(ctx->y1);
173: VecPlaceArray(ctx->y2,pd+m);
174: if (!ctx->diag) {
175: /* compute diagonal from rows and store in ctx->diag */
176: VecDuplicate(ctx->y2,&ctx->diag);
177: MatGetSize(ctx->A,NULL,&N);
178: PetscCalloc2(N,&work1,N,&work2);
179: if (ctx->swapped) {
180: MatGetOwnershipRange(ctx->AT,&start,&end);
181: for (i=start;i<end;i++) {
182: MatGetRow(ctx->AT,i,&ncols,NULL,&vals);
183: for (j=0;j<ncols;j++) work1[i] += vals[j]*vals[j];
184: MatRestoreRow(ctx->AT,i,&ncols,NULL,&vals);
185: }
186: } else {
187: MatGetOwnershipRange(ctx->A,&start,&end);
188: for (i=start;i<end;i++) {
189: MatGetRow(ctx->A,i,&ncols,&cols,&vals);
190: for (j=0;j<ncols;j++) work1[cols[j]] += vals[j]*vals[j];
191: MatRestoreRow(ctx->A,i,&ncols,&cols,&vals);
192: }
193: }
194: PetscMPIIntCast(N,&len);
195: MPIU_Allreduce(work1,work2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)B));
196: VecGetOwnershipRange(ctx->diag,&start,&end);
197: VecGetArrayWrite(ctx->diag,&diag);
198: for (i=start;i<end;i++) diag[i-start] = work2[i];
199: VecRestoreArrayWrite(ctx->diag,&diag);
200: PetscFree2(work1,work2);
201: }
202: VecCopy(ctx->diag,ctx->y2);
203: VecResetArray(ctx->y2);
204: VecRestoreArrayWrite(d,&pd);
205: return 0;
206: }
208: static PetscErrorCode MatDestroy_ECross(Mat B)
209: {
210: SVD_CYCLIC_SHELL *ctx;
212: MatShellGetContext(B,&ctx);
213: VecDestroy(&ctx->x1);
214: VecDestroy(&ctx->x2);
215: VecDestroy(&ctx->y1);
216: VecDestroy(&ctx->y2);
217: VecDestroy(&ctx->diag);
218: VecDestroy(&ctx->w);
219: PetscFree(ctx);
220: return 0;
221: }
223: /*
224: Builds extended cross product matrix C = | I_m 0 |
225: | 0 AT*A |
226: t is an auxiliary Vec used to take the dimensions of the upper block
227: */
228: static PetscErrorCode SVDCyclicGetECrossMat(SVD svd,Mat A,Mat AT,Mat *C,Vec t)
229: {
230: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
231: SVD_CYCLIC_SHELL *ctx;
232: PetscInt i,M,N,m,n,Istart,Iend;
233: VecType vtype;
234: Mat Id,Zm,Zn,ATA;
235: #if defined(PETSC_HAVE_CUDA)
236: PetscBool cuda;
237: #endif
239: MatGetSize(A,NULL,&N);
240: MatGetLocalSize(A,NULL,&n);
241: VecGetSize(t,&M);
242: VecGetLocalSize(t,&m);
244: if (cyclic->explicitmatrix) {
246: MatCreateConstantDiagonal(PetscObjectComm((PetscObject)svd),m,m,M,M,1.0,&Id);
247: MatCreate(PetscObjectComm((PetscObject)svd),&Zm);
248: MatSetSizes(Zm,m,n,M,N);
249: MatSetFromOptions(Zm);
250: MatSetUp(Zm);
251: MatGetOwnershipRange(Zm,&Istart,&Iend);
252: for (i=Istart;i<Iend;i++) {
253: if (i<N) MatSetValue(Zm,i,i,0.0,INSERT_VALUES);
254: }
255: MatAssemblyBegin(Zm,MAT_FINAL_ASSEMBLY);
256: MatAssemblyEnd(Zm,MAT_FINAL_ASSEMBLY);
257: MatCreate(PetscObjectComm((PetscObject)svd),&Zn);
258: MatSetSizes(Zn,n,m,N,M);
259: MatSetFromOptions(Zn);
260: MatSetUp(Zn);
261: MatGetOwnershipRange(Zn,&Istart,&Iend);
262: for (i=Istart;i<Iend;i++) {
263: if (i<m) MatSetValue(Zn,i,i,0.0,INSERT_VALUES);
264: }
265: MatAssemblyBegin(Zn,MAT_FINAL_ASSEMBLY);
266: MatAssemblyEnd(Zn,MAT_FINAL_ASSEMBLY);
267: MatProductCreate(AT,A,NULL,&ATA);
268: MatProductSetType(ATA,MATPRODUCT_AB);
269: MatProductSetFromOptions(ATA);
270: MatProductSymbolic(ATA);
271: MatProductNumeric(ATA);
272: MatCreateTile(1.0,Id,1.0,Zm,1.0,Zn,1.0,ATA,C);
273: MatDestroy(&Id);
274: MatDestroy(&Zm);
275: MatDestroy(&Zn);
276: MatDestroy(&ATA);
277: } else {
278: PetscNew(&ctx);
279: ctx->A = A;
280: ctx->AT = AT;
281: ctx->swapped = svd->swapped;
282: VecDuplicateEmpty(t,&ctx->x1);
283: VecDuplicateEmpty(t,&ctx->y1);
284: MatCreateVecsEmpty(A,&ctx->x2,NULL);
285: MatCreateVecsEmpty(A,&ctx->y2,NULL);
286: MatCreateVecs(A,NULL,&ctx->w);
287: MatCreateShell(PetscObjectComm((PetscObject)svd),m+n,m+n,M+N,M+N,ctx,C);
288: MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ECross);
289: MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_ECross);
290: #if defined(PETSC_HAVE_CUDA)
291: PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&cuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,"");
292: if (cuda) MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross_CUDA);
293: else
294: #endif
295: MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross);
296: MatGetVecType(A,&vtype);
297: MatSetVecType(*C,vtype);
298: }
299: return 0;
300: }
302: /* Convergence test relative to the norm of R (used in GSVD only) */
303: static PetscErrorCode EPSConv_Cyclic(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
304: {
305: SVD svd = (SVD)ctx;
307: *errest = res/PetscMax(svd->nrma,svd->nrmb);
308: return 0;
309: }
311: PetscErrorCode SVDSetUp_Cyclic(SVD svd)
312: {
313: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
314: PetscInt M,N,m,n,p,k,i,isl,offset,nev,ncv,mpd,maxit;
315: PetscReal tol;
316: const PetscScalar *isa,*oa;
317: PetscScalar *va;
318: EPSProblemType ptype;
319: PetscBool trackall,issinv;
320: Vec v,t;
321: ST st;
322: Mat Omega;
323: MatType Atype;
325: MatGetSize(svd->A,&M,&N);
326: MatGetLocalSize(svd->A,&m,&n);
327: if (!cyclic->eps) SVDCyclicGetEPS(svd,&cyclic->eps);
328: MatDestroy(&cyclic->C);
329: MatDestroy(&cyclic->D);
330: if (svd->isgeneralized) {
331: if (svd->which==SVD_SMALLEST) { /* alternative pencil */
332: MatCreateVecs(svd->B,NULL,&t);
333: SVDCyclicGetCyclicMat(svd,svd->B,svd->BT,&cyclic->C);
334: SVDCyclicGetECrossMat(svd,svd->A,svd->AT,&cyclic->D,t);
335: } else {
336: MatCreateVecs(svd->A,NULL,&t);
337: SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C);
338: SVDCyclicGetECrossMat(svd,svd->B,svd->BT,&cyclic->D,t);
339: }
340: VecDestroy(&t);
341: EPSSetOperators(cyclic->eps,cyclic->C,cyclic->D);
342: EPSGetProblemType(cyclic->eps,&ptype);
343: if (!ptype) EPSSetProblemType(cyclic->eps,EPS_GHEP);
344: } else if (svd->ishyperbolic) {
345: SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C);
346: MatCreateVecs(cyclic->C,&v,NULL);
347: VecSet(v,1.0);
348: VecGetArrayRead(svd->omega,&oa);
349: VecGetArray(v,&va);
350: if (svd->swapped) PetscArraycpy(va+m,oa,n);
351: else PetscArraycpy(va,oa,m);
352: VecRestoreArrayRead(svd->omega,&oa);
353: VecRestoreArray(v,&va);
354: MatGetType(svd->OP,&Atype);
355: MatCreate(PetscObjectComm((PetscObject)svd),&Omega);
356: MatSetSizes(Omega,m+n,m+n,M+N,M+N);
357: MatSetType(Omega,Atype);
358: MatSetUp(Omega);
359: MatDiagonalSet(Omega,v,INSERT_VALUES);
360: EPSSetOperators(cyclic->eps,cyclic->C,Omega);
361: EPSSetProblemType(cyclic->eps,EPS_GHIEP);
362: MatDestroy(&Omega);
363: VecDestroy(&v);
364: } else {
365: SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C);
366: EPSSetOperators(cyclic->eps,cyclic->C,NULL);
367: EPSSetProblemType(cyclic->eps,EPS_HEP);
368: }
369: if (!cyclic->usereps) {
370: if (svd->which == SVD_LARGEST) {
371: EPSGetST(cyclic->eps,&st);
372: PetscObjectTypeCompare((PetscObject)st,STSINVERT,&issinv);
373: if (issinv) EPSSetWhichEigenpairs(cyclic->eps,EPS_TARGET_MAGNITUDE);
374: else if (svd->ishyperbolic) EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_MAGNITUDE);
375: else EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL);
376: } else {
377: if (svd->isgeneralized) { /* computes sigma^{-1} via alternative pencil */
378: EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL);
379: } else {
380: if (svd->ishyperbolic) EPSSetWhichEigenpairs(cyclic->eps,EPS_TARGET_MAGNITUDE);
381: else EPSSetEigenvalueComparison(cyclic->eps,SlepcCompareSmallestPosReal,NULL);
382: EPSSetTarget(cyclic->eps,0.0);
383: }
384: }
385: EPSGetDimensions(cyclic->eps,&nev,&ncv,&mpd);
387: nev = PetscMax(nev,2*svd->nsv);
388: if (ncv==PETSC_DEFAULT && svd->ncv!=PETSC_DEFAULT) ncv = PetscMax(3*svd->nsv,svd->ncv);
389: if (mpd==PETSC_DEFAULT && svd->mpd!=PETSC_DEFAULT) mpd = svd->mpd;
390: EPSSetDimensions(cyclic->eps,nev,ncv,mpd);
391: EPSGetTolerances(cyclic->eps,&tol,&maxit);
392: if (tol==PETSC_DEFAULT) tol = svd->tol==PETSC_DEFAULT? SLEPC_DEFAULT_TOL/10.0: svd->tol;
393: if (maxit==PETSC_DEFAULT && svd->max_it!=PETSC_DEFAULT) maxit = svd->max_it;
394: EPSSetTolerances(cyclic->eps,tol,maxit);
395: switch (svd->conv) {
396: case SVD_CONV_ABS:
397: EPSSetConvergenceTest(cyclic->eps,EPS_CONV_ABS);break;
398: case SVD_CONV_REL:
399: EPSSetConvergenceTest(cyclic->eps,EPS_CONV_REL);break;
400: case SVD_CONV_NORM:
401: if (svd->isgeneralized) {
402: if (!svd->nrma) MatNorm(svd->OP,NORM_INFINITY,&svd->nrma);
403: if (!svd->nrmb) MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb);
404: EPSSetConvergenceTestFunction(cyclic->eps,EPSConv_Cyclic,svd,NULL);
405: } else {
406: EPSSetConvergenceTest(cyclic->eps,EPS_CONV_NORM);break;
407: }
408: break;
409: case SVD_CONV_MAXIT:
410: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Maxit convergence test not supported in this solver");
411: case SVD_CONV_USER:
412: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
413: }
414: }
415: SVDCheckUnsupported(svd,SVD_FEATURE_STOPPING);
416: /* Transfer the trackall option from svd to eps */
417: SVDGetTrackAll(svd,&trackall);
418: EPSSetTrackAll(cyclic->eps,trackall);
419: /* Transfer the initial subspace from svd to eps */
420: if (svd->nini<0 || svd->ninil<0) {
421: for (i=0;i<-PetscMin(svd->nini,svd->ninil);i++) {
422: MatCreateVecs(cyclic->C,&v,NULL);
423: VecGetArrayWrite(v,&va);
424: if (svd->isgeneralized) MatGetLocalSize(svd->B,&p,NULL);
425: k = (svd->isgeneralized && svd->which==SVD_SMALLEST)? p: m; /* size of upper block row */
426: if (i<-svd->ninil) {
427: VecGetArrayRead(svd->ISL[i],&isa);
428: if (svd->isgeneralized) {
429: VecGetLocalSize(svd->ISL[i],&isl);
431: offset = (svd->which==SVD_SMALLEST)? m: 0;
432: PetscArraycpy(va,isa+offset,k);
433: } else {
434: VecGetLocalSize(svd->ISL[i],&isl);
436: PetscArraycpy(va,isa,k);
437: }
438: VecRestoreArrayRead(svd->IS[i],&isa);
439: } else PetscArrayzero(&va,k);
440: if (i<-svd->nini) {
441: VecGetLocalSize(svd->IS[i],&isl);
443: VecGetArrayRead(svd->IS[i],&isa);
444: PetscArraycpy(va+k,isa,n);
445: VecRestoreArrayRead(svd->IS[i],&isa);
446: } else PetscArrayzero(va+k,n);
447: VecRestoreArrayWrite(v,&va);
448: VecDestroy(&svd->IS[i]);
449: svd->IS[i] = v;
450: }
451: svd->nini = PetscMin(svd->nini,svd->ninil);
452: EPSSetInitialSpace(cyclic->eps,-svd->nini,svd->IS);
453: SlepcBasisDestroy_Private(&svd->nini,&svd->IS);
454: SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL);
455: }
456: EPSSetUp(cyclic->eps);
457: EPSGetDimensions(cyclic->eps,NULL,&svd->ncv,&svd->mpd);
458: svd->ncv = PetscMin(svd->ncv,PetscMin(M,N));
459: EPSGetTolerances(cyclic->eps,NULL,&svd->max_it);
460: if (svd->tol==PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;
462: svd->leftbasis = PETSC_TRUE;
463: SVDAllocateSolution(svd,0);
464: return 0;
465: }
467: PetscErrorCode SVDCyclicCheckEigenvalue(SVD svd,PetscScalar er,PetscScalar ei,PetscReal *sigma,PetscBool *isreal)
468: {
469: if (svd->ishyperbolic && PetscDefined(USE_COMPLEX) && PetscAbsReal(PetscImaginaryPart(er))>10*PetscAbsReal(PetscRealPart(er))) {
470: *sigma = PetscImaginaryPart(er);
471: if (isreal) *isreal = PETSC_FALSE;
472: } else if (svd->ishyperbolic && !PetscDefined(USE_COMPLEX) && PetscAbsScalar(ei)>10*PetscAbsScalar(er)) {
473: *sigma = PetscRealPart(ei);
474: if (isreal) *isreal = PETSC_FALSE;
475: } else {
476: *sigma = PetscRealPart(er);
477: if (isreal) *isreal = PETSC_TRUE;
478: }
479: return 0;
480: }
482: PetscErrorCode SVDSolve_Cyclic(SVD svd)
483: {
484: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
485: PetscInt i,j,nconv;
486: PetscScalar er,ei;
487: PetscReal sigma;
489: EPSSolve(cyclic->eps);
490: EPSGetConverged(cyclic->eps,&nconv);
491: EPSGetIterationNumber(cyclic->eps,&svd->its);
492: EPSGetConvergedReason(cyclic->eps,(EPSConvergedReason*)&svd->reason);
493: for (i=0,j=0;i<nconv;i++) {
494: EPSGetEigenvalue(cyclic->eps,i,&er,&ei);
495: SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL);
496: if (sigma>0.0) {
497: if (svd->isgeneralized && svd->which==SVD_SMALLEST) svd->sigma[j] = 1.0/sigma;
498: else svd->sigma[j] = sigma;
499: j++;
500: }
501: }
502: svd->nconv = j;
503: return 0;
504: }
506: static PetscErrorCode SVDComputeVectors_Cyclic_Standard(SVD svd)
507: {
508: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
509: PetscInt i,j,m,nconv;
510: PetscScalar er,ei;
511: PetscReal sigma;
512: const PetscScalar *px;
513: Vec x,x1,x2;
515: MatCreateVecs(cyclic->C,&x,NULL);
516: MatGetLocalSize(svd->A,&m,NULL);
517: MatCreateVecsEmpty(svd->A,&x2,&x1);
518: EPSGetConverged(cyclic->eps,&nconv);
519: for (i=0,j=0;i<nconv;i++) {
520: EPSGetEigenpair(cyclic->eps,i,&er,&ei,x,NULL);
521: SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL);
522: if (sigma<0.0) continue;
523: VecGetArrayRead(x,&px);
524: VecPlaceArray(x1,px);
525: VecPlaceArray(x2,px+m);
526: BVInsertVec(svd->U,j,x1);
527: BVScaleColumn(svd->U,j,PETSC_SQRT2);
528: BVInsertVec(svd->V,j,x2);
529: BVScaleColumn(svd->V,j,PETSC_SQRT2);
530: VecResetArray(x1);
531: VecResetArray(x2);
532: VecRestoreArrayRead(x,&px);
533: j++;
534: }
535: VecDestroy(&x);
536: VecDestroy(&x1);
537: VecDestroy(&x2);
538: return 0;
539: }
541: static PetscErrorCode SVDComputeVectors_Cyclic_Generalized(SVD svd)
542: {
543: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
544: PetscInt i,j,m,p,nconv;
545: PetscScalar *dst,er,ei;
546: PetscReal sigma;
547: const PetscScalar *src,*px;
548: Vec u,v,x,x1,x2,uv;
550: MatGetLocalSize(svd->A,&m,NULL);
551: MatGetLocalSize(svd->B,&p,NULL);
552: MatCreateVecs(cyclic->C,&x,NULL);
553: if (svd->which==SVD_SMALLEST) MatCreateVecsEmpty(svd->B,&x1,&x2);
554: else MatCreateVecsEmpty(svd->A,&x2,&x1);
555: MatCreateVecs(svd->A,NULL,&u);
556: MatCreateVecs(svd->B,NULL,&v);
557: EPSGetConverged(cyclic->eps,&nconv);
558: for (i=0,j=0;i<nconv;i++) {
559: EPSGetEigenpair(cyclic->eps,i,&er,&ei,x,NULL);
560: SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL);
561: if (sigma<0.0) continue;
562: if (svd->which==SVD_SMALLEST) {
563: /* evec_i = 1/sqrt(2)*[ v_i; w_i ], w_i = x_i/c_i */
564: VecGetArrayRead(x,&px);
565: VecPlaceArray(x2,px);
566: VecPlaceArray(x1,px+p);
567: VecCopy(x2,v);
568: VecScale(v,PETSC_SQRT2); /* v_i = sqrt(2)*evec_i_1 */
569: VecScale(x1,PETSC_SQRT2); /* w_i = sqrt(2)*evec_i_2 */
570: MatMult(svd->A,x1,u); /* A*w_i = u_i */
571: VecScale(x1,1.0/PetscSqrtScalar(1.0+sigma*sigma)); /* x_i = w_i*c_i */
572: BVInsertVec(svd->V,j,x1);
573: VecResetArray(x2);
574: VecResetArray(x1);
575: VecRestoreArrayRead(x,&px);
576: } else {
577: /* evec_i = 1/sqrt(2)*[ u_i; w_i ], w_i = x_i/s_i */
578: VecGetArrayRead(x,&px);
579: VecPlaceArray(x1,px);
580: VecPlaceArray(x2,px+m);
581: VecCopy(x1,u);
582: VecScale(u,PETSC_SQRT2); /* u_i = sqrt(2)*evec_i_1 */
583: VecScale(x2,PETSC_SQRT2); /* w_i = sqrt(2)*evec_i_2 */
584: MatMult(svd->B,x2,v); /* B*w_i = v_i */
585: VecScale(x2,1.0/PetscSqrtScalar(1.0+sigma*sigma)); /* x_i = w_i*s_i */
586: BVInsertVec(svd->V,j,x2);
587: VecResetArray(x1);
588: VecResetArray(x2);
589: VecRestoreArrayRead(x,&px);
590: }
591: /* copy [u;v] to U[j] */
592: BVGetColumn(svd->U,j,&uv);
593: VecGetArrayWrite(uv,&dst);
594: VecGetArrayRead(u,&src);
595: PetscArraycpy(dst,src,m);
596: VecRestoreArrayRead(u,&src);
597: VecGetArrayRead(v,&src);
598: PetscArraycpy(dst+m,src,p);
599: VecRestoreArrayRead(v,&src);
600: VecRestoreArrayWrite(uv,&dst);
601: BVRestoreColumn(svd->U,j,&uv);
602: j++;
603: }
604: VecDestroy(&x);
605: VecDestroy(&x1);
606: VecDestroy(&x2);
607: VecDestroy(&u);
608: VecDestroy(&v);
609: return 0;
610: }
612: #if defined(PETSC_USE_COMPLEX)
613: /* VecMaxAbs: returns the entry of x that has max(abs(x(i))), using w as a workspace vector */
614: static PetscErrorCode VecMaxAbs(Vec x,Vec w,PetscScalar *v)
615: {
616: PetscMPIInt size,rank,root;
617: const PetscScalar *xx;
618: const PetscInt *ranges;
619: PetscReal val;
620: PetscInt p;
622: VecCopy(x,w);
623: VecAbs(w);
624: VecMax(w,&p,&val);
625: MPI_Comm_size(PetscObjectComm((PetscObject)x),&size);
626: MPI_Comm_rank(PetscObjectComm((PetscObject)x),&rank);
627: VecGetOwnershipRanges(x,&ranges);
628: for (root=0;root<size;root++) if (p>=ranges[root] && p<ranges[root+1]) break;
629: if (rank==root) {
630: VecGetArrayRead(x,&xx);
631: *v = xx[p-ranges[root]];
632: VecRestoreArrayRead(x,&xx);
633: }
634: MPI_Bcast(v,1,MPIU_SCALAR,root,PetscObjectComm((PetscObject)x));
635: return 0;
636: }
637: #endif
639: static PetscErrorCode SVDComputeVectors_Cyclic_Hyperbolic(SVD svd)
640: {
641: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
642: PetscInt i,j,m,n,N,nconv;
643: PetscScalar er,ei;
644: PetscReal sigma,nrm;
645: PetscBool isreal;
646: const PetscScalar *px;
647: Vec u,x,xi=NULL,x1,x2,x1i=NULL,x2i;
648: Mat Omega;
649: MatType Atype;
650: BV U=NULL,V=NULL;
651: #if !defined(PETSC_USE_COMPLEX)
652: const PetscScalar *pxi;
653: PetscReal nrmr,nrmi;
654: #else
655: PetscScalar alpha;
656: #endif
658: MatCreateVecs(cyclic->C,&x,svd->ishyperbolic?&xi:NULL);
659: MatGetLocalSize(svd->A,&m,NULL);
660: MatCreateVecsEmpty(svd->OP,&x2,&x1);
661: #if defined(PETSC_USE_COMPLEX)
662: MatCreateVecs(svd->OP,&x2i,&x1i);
663: #else
664: MatCreateVecsEmpty(svd->OP,&x2i,&x1i);
665: #endif
666: /* set-up Omega-normalization of U */
667: U = svd->swapped? svd->V: svd->U;
668: V = svd->swapped? svd->U: svd->V;
669: MatGetType(svd->A,&Atype);
670: BVGetSizes(U,&n,&N,NULL);
671: MatCreate(PetscObjectComm((PetscObject)svd),&Omega);
672: MatSetSizes(Omega,n,n,N,N);
673: MatSetType(Omega,Atype);
674: MatSetUp(Omega);
675: MatDiagonalSet(Omega,svd->omega,INSERT_VALUES);
676: BVSetMatrix(U,Omega,PETSC_TRUE);
677: MatDestroy(&Omega);
679: EPSGetConverged(cyclic->eps,&nconv);
680: for (i=0,j=0;i<nconv;i++) {
681: EPSGetEigenpair(cyclic->eps,i,&er,&ei,x,xi);
682: SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,&isreal);
683: if (sigma<0.0) continue;
684: VecGetArrayRead(x,&px);
685: if (svd->swapped) {
686: VecPlaceArray(x2,px);
687: VecPlaceArray(x1,px+m);
688: } else {
689: VecPlaceArray(x1,px);
690: VecPlaceArray(x2,px+n);
691: }
692: #if defined(PETSC_USE_COMPLEX)
693: BVInsertVec(U,j,x1);
694: BVInsertVec(V,j,x2);
695: if (!isreal) {
696: VecMaxAbs(x1,x1i,&alpha);
697: BVScaleColumn(U,j,PetscAbsScalar(alpha)/alpha);
698: BVScaleColumn(V,j,PetscAbsScalar(alpha)/(alpha*PETSC_i));
699: }
700: #else
701: VecGetArrayRead(xi,&pxi);
702: if (svd->swapped) {
703: VecPlaceArray(x2i,pxi);
704: VecPlaceArray(x1i,pxi+m);
705: } else {
706: VecPlaceArray(x1i,pxi);
707: VecPlaceArray(x2i,pxi+n);
708: }
709: VecNorm(x2,NORM_2,&nrmr);
710: VecNorm(x2i,NORM_2,&nrmi);
711: if (nrmi>nrmr) {
712: if (isreal) {
713: BVInsertVec(U,j,x1i);
714: BVInsertVec(V,j,x2i);
715: } else {
716: BVInsertVec(U,j,x1);
717: BVInsertVec(V,j,x2i);
718: }
719: } else {
720: if (isreal) {
721: BVInsertVec(U,j,x1);
722: BVInsertVec(V,j,x2);
723: } else {
724: BVInsertVec(U,j,x1i);
725: BVScaleColumn(U,j,-1.0);
726: BVInsertVec(V,j,x2);
727: }
728: }
729: VecResetArray(x1i);
730: VecResetArray(x2i);
731: VecRestoreArrayRead(xi,&pxi);
732: #endif
733: VecResetArray(x1);
734: VecResetArray(x2);
735: VecRestoreArrayRead(x,&px);
736: BVGetColumn(U,j,&u);
737: VecPointwiseMult(u,u,svd->omega);
738: BVRestoreColumn(U,j,&u);
739: BVNormColumn(U,j,NORM_2,&nrm);
740: BVScaleColumn(U,j,1.0/PetscAbs(nrm));
741: svd->sign[j] = PetscSign(nrm);
742: BVNormColumn(V,j,NORM_2,&nrm);
743: BVScaleColumn(V,j,1.0/nrm);
744: j++;
745: }
746: VecDestroy(&x);
747: VecDestroy(&x1);
748: VecDestroy(&x2);
749: VecDestroy(&xi);
750: VecDestroy(&x1i);
751: VecDestroy(&x2i);
752: return 0;
753: }
755: PetscErrorCode SVDComputeVectors_Cyclic(SVD svd)
756: {
757: switch (svd->problem_type) {
758: case SVD_STANDARD:
759: SVDComputeVectors_Cyclic_Standard(svd);
760: break;
761: case SVD_GENERALIZED:
762: SVDComputeVectors_Cyclic_Generalized(svd);
763: break;
764: case SVD_HYPERBOLIC:
765: SVDComputeVectors_Cyclic_Hyperbolic(svd);
766: break;
767: default:
768: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONG,"Unknown singular value problem type");
769: }
770: return 0;
771: }
773: static PetscErrorCode EPSMonitor_Cyclic(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
774: {
775: PetscInt i,j;
776: SVD svd = (SVD)ctx;
777: PetscScalar er,ei;
778: PetscReal sigma;
780: nconv = 0;
781: for (i=0,j=0;i<PetscMin(nest,svd->ncv);i++) {
782: er = eigr[i]; ei = eigi[i];
783: STBackTransform(eps->st,1,&er,&ei);
784: SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL);
785: if (sigma>0.0) {
786: svd->sigma[j] = sigma;
787: svd->errest[j] = errest[i];
788: if (errest[i] && errest[i] < svd->tol) nconv++;
789: j++;
790: }
791: }
792: nest = j;
793: SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest);
794: return 0;
795: }
797: PetscErrorCode SVDSetFromOptions_Cyclic(SVD svd,PetscOptionItems *PetscOptionsObject)
798: {
799: PetscBool set,val;
800: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
801: ST st;
803: PetscOptionsHeadBegin(PetscOptionsObject,"SVD Cyclic Options");
805: PetscOptionsBool("-svd_cyclic_explicitmatrix","Use cyclic explicit matrix","SVDCyclicSetExplicitMatrix",cyclic->explicitmatrix,&val,&set);
806: if (set) SVDCyclicSetExplicitMatrix(svd,val);
808: PetscOptionsHeadEnd();
810: if (!cyclic->eps) SVDCyclicGetEPS(svd,&cyclic->eps);
811: if (!cyclic->explicitmatrix && !cyclic->usereps) {
812: /* use as default an ST with shell matrix and Jacobi */
813: EPSGetST(cyclic->eps,&st);
814: STSetMatMode(st,ST_MATMODE_SHELL);
815: }
816: EPSSetFromOptions(cyclic->eps);
817: return 0;
818: }
820: static PetscErrorCode SVDCyclicSetExplicitMatrix_Cyclic(SVD svd,PetscBool explicitmat)
821: {
822: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
824: if (cyclic->explicitmatrix != explicitmat) {
825: cyclic->explicitmatrix = explicitmat;
826: svd->state = SVD_STATE_INITIAL;
827: }
828: return 0;
829: }
831: /*@
832: SVDCyclicSetExplicitMatrix - Indicate if the eigensolver operator
833: H(A) = [ 0 A ; A^T 0 ] must be computed explicitly.
835: Logically Collective on svd
837: Input Parameters:
838: + svd - singular value solver
839: - explicitmat - boolean flag indicating if H(A) is built explicitly
841: Options Database Key:
842: . -svd_cyclic_explicitmatrix <boolean> - Indicates the boolean flag
844: Level: advanced
846: .seealso: SVDCyclicGetExplicitMatrix()
847: @*/
848: PetscErrorCode SVDCyclicSetExplicitMatrix(SVD svd,PetscBool explicitmat)
849: {
852: PetscTryMethod(svd,"SVDCyclicSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
853: return 0;
854: }
856: static PetscErrorCode SVDCyclicGetExplicitMatrix_Cyclic(SVD svd,PetscBool *explicitmat)
857: {
858: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
860: *explicitmat = cyclic->explicitmatrix;
861: return 0;
862: }
864: /*@
865: SVDCyclicGetExplicitMatrix - Returns the flag indicating if H(A) is built explicitly.
867: Not Collective
869: Input Parameter:
870: . svd - singular value solver
872: Output Parameter:
873: . explicitmat - the mode flag
875: Level: advanced
877: .seealso: SVDCyclicSetExplicitMatrix()
878: @*/
879: PetscErrorCode SVDCyclicGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
880: {
883: PetscUseMethod(svd,"SVDCyclicGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
884: return 0;
885: }
887: static PetscErrorCode SVDCyclicSetEPS_Cyclic(SVD svd,EPS eps)
888: {
889: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
891: PetscObjectReference((PetscObject)eps);
892: EPSDestroy(&cyclic->eps);
893: cyclic->eps = eps;
894: cyclic->usereps = PETSC_TRUE;
895: svd->state = SVD_STATE_INITIAL;
896: return 0;
897: }
899: /*@
900: SVDCyclicSetEPS - Associate an eigensolver object (EPS) to the
901: singular value solver.
903: Collective on svd
905: Input Parameters:
906: + svd - singular value solver
907: - eps - the eigensolver object
909: Level: advanced
911: .seealso: SVDCyclicGetEPS()
912: @*/
913: PetscErrorCode SVDCyclicSetEPS(SVD svd,EPS eps)
914: {
918: PetscTryMethod(svd,"SVDCyclicSetEPS_C",(SVD,EPS),(svd,eps));
919: return 0;
920: }
922: static PetscErrorCode SVDCyclicGetEPS_Cyclic(SVD svd,EPS *eps)
923: {
924: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
926: if (!cyclic->eps) {
927: EPSCreate(PetscObjectComm((PetscObject)svd),&cyclic->eps);
928: PetscObjectIncrementTabLevel((PetscObject)cyclic->eps,(PetscObject)svd,1);
929: EPSSetOptionsPrefix(cyclic->eps,((PetscObject)svd)->prefix);
930: EPSAppendOptionsPrefix(cyclic->eps,"svd_cyclic_");
931: PetscObjectSetOptions((PetscObject)cyclic->eps,((PetscObject)svd)->options);
932: EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL);
933: EPSMonitorSet(cyclic->eps,EPSMonitor_Cyclic,svd,NULL);
934: }
935: *eps = cyclic->eps;
936: return 0;
937: }
939: /*@
940: SVDCyclicGetEPS - Retrieve the eigensolver object (EPS) associated
941: to the singular value solver.
943: Not Collective
945: Input Parameter:
946: . svd - singular value solver
948: Output Parameter:
949: . eps - the eigensolver object
951: Level: advanced
953: .seealso: SVDCyclicSetEPS()
954: @*/
955: PetscErrorCode SVDCyclicGetEPS(SVD svd,EPS *eps)
956: {
959: PetscUseMethod(svd,"SVDCyclicGetEPS_C",(SVD,EPS*),(svd,eps));
960: return 0;
961: }
963: PetscErrorCode SVDView_Cyclic(SVD svd,PetscViewer viewer)
964: {
965: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
966: PetscBool isascii;
968: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
969: if (isascii) {
970: if (!cyclic->eps) SVDCyclicGetEPS(svd,&cyclic->eps);
971: PetscViewerASCIIPrintf(viewer," %s matrix\n",cyclic->explicitmatrix?"explicit":"implicit");
972: PetscViewerASCIIPushTab(viewer);
973: EPSView(cyclic->eps,viewer);
974: PetscViewerASCIIPopTab(viewer);
975: }
976: return 0;
977: }
979: PetscErrorCode SVDReset_Cyclic(SVD svd)
980: {
981: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
983: EPSReset(cyclic->eps);
984: MatDestroy(&cyclic->C);
985: MatDestroy(&cyclic->D);
986: return 0;
987: }
989: PetscErrorCode SVDDestroy_Cyclic(SVD svd)
990: {
991: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
993: EPSDestroy(&cyclic->eps);
994: PetscFree(svd->data);
995: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",NULL);
996: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",NULL);
997: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",NULL);
998: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",NULL);
999: return 0;
1000: }
1002: SLEPC_EXTERN PetscErrorCode SVDCreate_Cyclic(SVD svd)
1003: {
1004: SVD_CYCLIC *cyclic;
1006: PetscNew(&cyclic);
1007: svd->data = (void*)cyclic;
1008: svd->ops->solve = SVDSolve_Cyclic;
1009: svd->ops->solveg = SVDSolve_Cyclic;
1010: svd->ops->solveh = SVDSolve_Cyclic;
1011: svd->ops->setup = SVDSetUp_Cyclic;
1012: svd->ops->setfromoptions = SVDSetFromOptions_Cyclic;
1013: svd->ops->destroy = SVDDestroy_Cyclic;
1014: svd->ops->reset = SVDReset_Cyclic;
1015: svd->ops->view = SVDView_Cyclic;
1016: svd->ops->computevectors = SVDComputeVectors_Cyclic;
1017: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",SVDCyclicSetEPS_Cyclic);
1018: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",SVDCyclicGetEPS_Cyclic);
1019: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",SVDCyclicSetExplicitMatrix_Cyclic);
1020: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",SVDCyclicGetExplicitMatrix_Cyclic);
1021: return 0;
1022: }