Actual source code: cyclic.c
slepc-3.17.2 2022-08-09
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: PetscFunctionReturn(0);
44: }
46: static PetscErrorCode MatGetDiagonal_Cyclic(Mat B,Vec diag)
47: {
48: VecSet(diag,0.0);
49: PetscFunctionReturn(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: PetscFunctionReturn(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: PetscLogObjectParent((PetscObject)svd,(PetscObject)ctx->x1);
112: PetscLogObjectParent((PetscObject)svd,(PetscObject)ctx->x2);
113: PetscLogObjectParent((PetscObject)svd,(PetscObject)ctx->y1);
114: PetscLogObjectParent((PetscObject)svd,(PetscObject)ctx->y2);
115: MatCreateShell(PetscObjectComm((PetscObject)svd),m+n,m+n,M+N,M+N,ctx,C);
116: MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Cyclic);
117: MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_Cyclic);
118: #if defined(PETSC_HAVE_CUDA)
119: PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&cuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,"");
120: if (cuda) MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic_CUDA);
121: else
122: #endif
123: MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic);
124: MatGetVecType(A,&vtype);
125: MatSetVecType(*C,vtype);
126: }
127: PetscLogObjectParent((PetscObject)svd,(PetscObject)*C);
128: PetscFunctionReturn(0);
129: }
131: static PetscErrorCode MatMult_ECross(Mat B,Vec x,Vec y)
132: {
133: SVD_CYCLIC_SHELL *ctx;
134: const PetscScalar *px;
135: PetscScalar *py;
136: PetscInt mn,m,n;
138: MatShellGetContext(B,&ctx);
139: MatGetLocalSize(ctx->A,NULL,&n);
140: VecGetLocalSize(y,&mn);
141: m = mn-n;
142: VecGetArrayRead(x,&px);
143: VecGetArrayWrite(y,&py);
144: VecPlaceArray(ctx->x1,px);
145: VecPlaceArray(ctx->x2,px+m);
146: VecPlaceArray(ctx->y1,py);
147: VecPlaceArray(ctx->y2,py+m);
148: VecCopy(ctx->x1,ctx->y1);
149: MatMult(ctx->A,ctx->x2,ctx->w);
150: MatMult(ctx->AT,ctx->w,ctx->y2);
151: VecResetArray(ctx->x1);
152: VecResetArray(ctx->x2);
153: VecResetArray(ctx->y1);
154: VecResetArray(ctx->y2);
155: VecRestoreArrayRead(x,&px);
156: VecRestoreArrayWrite(y,&py);
157: PetscFunctionReturn(0);
158: }
160: static PetscErrorCode MatGetDiagonal_ECross(Mat B,Vec d)
161: {
162: SVD_CYCLIC_SHELL *ctx;
163: PetscScalar *pd;
164: PetscMPIInt len;
165: PetscInt mn,m,n,N,i,j,start,end,ncols;
166: PetscScalar *work1,*work2,*diag;
167: const PetscInt *cols;
168: const PetscScalar *vals;
170: MatShellGetContext(B,&ctx);
171: MatGetLocalSize(ctx->A,NULL,&n);
172: VecGetLocalSize(d,&mn);
173: m = mn-n;
174: VecGetArrayWrite(d,&pd);
175: VecPlaceArray(ctx->y1,pd);
176: VecSet(ctx->y1,1.0);
177: VecResetArray(ctx->y1);
178: VecPlaceArray(ctx->y2,pd+m);
179: if (!ctx->diag) {
180: /* compute diagonal from rows and store in ctx->diag */
181: VecDuplicate(ctx->y2,&ctx->diag);
182: MatGetSize(ctx->A,NULL,&N);
183: PetscCalloc2(N,&work1,N,&work2);
184: if (ctx->swapped) {
185: MatGetOwnershipRange(ctx->AT,&start,&end);
186: for (i=start;i<end;i++) {
187: MatGetRow(ctx->AT,i,&ncols,NULL,&vals);
188: for (j=0;j<ncols;j++) work1[i] += vals[j]*vals[j];
189: MatRestoreRow(ctx->AT,i,&ncols,NULL,&vals);
190: }
191: } else {
192: MatGetOwnershipRange(ctx->A,&start,&end);
193: for (i=start;i<end;i++) {
194: MatGetRow(ctx->A,i,&ncols,&cols,&vals);
195: for (j=0;j<ncols;j++) work1[cols[j]] += vals[j]*vals[j];
196: MatRestoreRow(ctx->A,i,&ncols,&cols,&vals);
197: }
198: }
199: PetscMPIIntCast(N,&len);
200: MPIU_Allreduce(work1,work2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)B));
201: VecGetOwnershipRange(ctx->diag,&start,&end);
202: VecGetArrayWrite(ctx->diag,&diag);
203: for (i=start;i<end;i++) diag[i-start] = work2[i];
204: VecRestoreArrayWrite(ctx->diag,&diag);
205: PetscFree2(work1,work2);
206: }
207: VecCopy(ctx->diag,ctx->y2);
208: VecResetArray(ctx->y2);
209: VecRestoreArrayWrite(d,&pd);
210: PetscFunctionReturn(0);
211: }
213: static PetscErrorCode MatDestroy_ECross(Mat B)
214: {
215: SVD_CYCLIC_SHELL *ctx;
217: MatShellGetContext(B,&ctx);
218: VecDestroy(&ctx->x1);
219: VecDestroy(&ctx->x2);
220: VecDestroy(&ctx->y1);
221: VecDestroy(&ctx->y2);
222: VecDestroy(&ctx->diag);
223: VecDestroy(&ctx->w);
224: PetscFree(ctx);
225: PetscFunctionReturn(0);
226: }
228: /*
229: Builds extended cross product matrix C = | I_m 0 |
230: | 0 AT*A |
231: t is an auxiliary Vec used to take the dimensions of the upper block
232: */
233: static PetscErrorCode SVDCyclicGetECrossMat(SVD svd,Mat A,Mat AT,Mat *C,Vec t)
234: {
235: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
236: SVD_CYCLIC_SHELL *ctx;
237: PetscInt i,M,N,m,n,Istart,Iend;
238: VecType vtype;
239: Mat Id,Zm,Zn,ATA;
240: #if defined(PETSC_HAVE_CUDA)
241: PetscBool cuda;
242: #endif
244: MatGetSize(A,NULL,&N);
245: MatGetLocalSize(A,NULL,&n);
246: VecGetSize(t,&M);
247: VecGetLocalSize(t,&m);
249: if (cyclic->explicitmatrix) {
251: MatCreateConstantDiagonal(PetscObjectComm((PetscObject)svd),m,m,M,M,1.0,&Id);
252: MatCreate(PetscObjectComm((PetscObject)svd),&Zm);
253: MatSetSizes(Zm,m,n,M,N);
254: MatSetFromOptions(Zm);
255: MatSetUp(Zm);
256: MatGetOwnershipRange(Zm,&Istart,&Iend);
257: for (i=Istart;i<Iend;i++) {
258: if (i<N) MatSetValue(Zm,i,i,0.0,INSERT_VALUES);
259: }
260: MatAssemblyBegin(Zm,MAT_FINAL_ASSEMBLY);
261: MatAssemblyEnd(Zm,MAT_FINAL_ASSEMBLY);
262: MatCreate(PetscObjectComm((PetscObject)svd),&Zn);
263: MatSetSizes(Zn,n,m,N,M);
264: MatSetFromOptions(Zn);
265: MatSetUp(Zn);
266: MatGetOwnershipRange(Zn,&Istart,&Iend);
267: for (i=Istart;i<Iend;i++) {
268: if (i<m) MatSetValue(Zn,i,i,0.0,INSERT_VALUES);
269: }
270: MatAssemblyBegin(Zn,MAT_FINAL_ASSEMBLY);
271: MatAssemblyEnd(Zn,MAT_FINAL_ASSEMBLY);
272: MatProductCreate(AT,A,NULL,&ATA);
273: MatProductSetType(ATA,MATPRODUCT_AB);
274: MatProductSetFromOptions(ATA);
275: MatProductSymbolic(ATA);
276: MatProductNumeric(ATA);
277: MatCreateTile(1.0,Id,1.0,Zm,1.0,Zn,1.0,ATA,C);
278: MatDestroy(&Id);
279: MatDestroy(&Zm);
280: MatDestroy(&Zn);
281: MatDestroy(&ATA);
282: } else {
283: PetscNew(&ctx);
284: ctx->A = A;
285: ctx->AT = AT;
286: ctx->swapped = svd->swapped;
287: VecDuplicateEmpty(t,&ctx->x1);
288: VecDuplicateEmpty(t,&ctx->y1);
289: MatCreateVecsEmpty(A,&ctx->x2,NULL);
290: MatCreateVecsEmpty(A,&ctx->y2,NULL);
291: MatCreateVecs(A,NULL,&ctx->w);
292: PetscLogObjectParent((PetscObject)svd,(PetscObject)ctx->x1);
293: PetscLogObjectParent((PetscObject)svd,(PetscObject)ctx->x2);
294: PetscLogObjectParent((PetscObject)svd,(PetscObject)ctx->y1);
295: PetscLogObjectParent((PetscObject)svd,(PetscObject)ctx->y2);
296: MatCreateShell(PetscObjectComm((PetscObject)svd),m+n,m+n,M+N,M+N,ctx,C);
297: MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ECross);
298: MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_ECross);
299: #if defined(PETSC_HAVE_CUDA)
300: PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&cuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,"");
301: if (cuda) MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross_CUDA);
302: else
303: #endif
304: MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross);
305: MatGetVecType(A,&vtype);
306: MatSetVecType(*C,vtype);
307: }
308: PetscLogObjectParent((PetscObject)svd,(PetscObject)*C);
309: PetscFunctionReturn(0);
310: }
312: /* Convergence test relative to the norm of R (used in GSVD only) */
313: static PetscErrorCode EPSConv_Cyclic(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
314: {
315: SVD svd = (SVD)ctx;
317: *errest = res/PetscMax(svd->nrma,svd->nrmb);
318: PetscFunctionReturn(0);
319: }
321: PetscErrorCode SVDSetUp_Cyclic(SVD svd)
322: {
323: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
324: PetscInt M,N,m,n,p,k,i,isl,offset,nev,ncv,mpd,maxit;
325: PetscReal tol;
326: const PetscScalar *isa;
327: PetscScalar *va;
328: EPSProblemType ptype;
329: PetscBool trackall,issinv;
330: Vec v,t;
331: ST st;
333: MatGetSize(svd->A,&M,&N);
334: MatGetLocalSize(svd->A,&m,&n);
335: if (!cyclic->eps) SVDCyclicGetEPS(svd,&cyclic->eps);
336: MatDestroy(&cyclic->C);
337: MatDestroy(&cyclic->D);
338: if (svd->isgeneralized) {
339: if (svd->which==SVD_SMALLEST) { /* alternative pencil */
340: MatCreateVecs(svd->B,NULL,&t);
341: SVDCyclicGetCyclicMat(svd,svd->B,svd->BT,&cyclic->C);
342: SVDCyclicGetECrossMat(svd,svd->A,svd->AT,&cyclic->D,t);
343: } else {
344: MatCreateVecs(svd->A,NULL,&t);
345: SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C);
346: SVDCyclicGetECrossMat(svd,svd->B,svd->BT,&cyclic->D,t);
347: }
348: VecDestroy(&t);
349: EPSSetOperators(cyclic->eps,cyclic->C,cyclic->D);
350: EPSGetProblemType(cyclic->eps,&ptype);
351: if (!ptype) EPSSetProblemType(cyclic->eps,EPS_GHEP);
352: } else {
353: SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C);
354: EPSSetOperators(cyclic->eps,cyclic->C,NULL);
355: EPSSetProblemType(cyclic->eps,EPS_HEP);
356: }
357: if (!cyclic->usereps) {
358: if (svd->which == SVD_LARGEST) {
359: EPSGetST(cyclic->eps,&st);
360: PetscObjectTypeCompare((PetscObject)st,STSINVERT,&issinv);
361: if (issinv) EPSSetWhichEigenpairs(cyclic->eps,EPS_TARGET_MAGNITUDE);
362: else EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL);
363: } else {
364: if (svd->isgeneralized) { /* computes sigma^{-1} via alternative pencil */
365: EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL);
366: } else {
367: EPSSetEigenvalueComparison(cyclic->eps,SlepcCompareSmallestPosReal,NULL);
368: EPSSetTarget(cyclic->eps,0.0);
369: }
370: }
371: EPSGetDimensions(cyclic->eps,&nev,&ncv,&mpd);
373: nev = PetscMax(nev,2*svd->nsv);
374: if (ncv==PETSC_DEFAULT && svd->ncv!=PETSC_DEFAULT) ncv = PetscMax(3*svd->nsv,svd->ncv);
375: if (mpd==PETSC_DEFAULT && svd->mpd!=PETSC_DEFAULT) mpd = svd->mpd;
376: EPSSetDimensions(cyclic->eps,nev,ncv,mpd);
377: EPSGetTolerances(cyclic->eps,&tol,&maxit);
378: if (tol==PETSC_DEFAULT) tol = svd->tol==PETSC_DEFAULT? SLEPC_DEFAULT_TOL/10.0: svd->tol;
379: if (maxit==PETSC_DEFAULT && svd->max_it!=PETSC_DEFAULT) maxit = svd->max_it;
380: EPSSetTolerances(cyclic->eps,tol,maxit);
381: switch (svd->conv) {
382: case SVD_CONV_ABS:
383: EPSSetConvergenceTest(cyclic->eps,EPS_CONV_ABS);break;
384: case SVD_CONV_REL:
385: EPSSetConvergenceTest(cyclic->eps,EPS_CONV_REL);break;
386: case SVD_CONV_NORM:
387: if (svd->isgeneralized) {
388: if (!svd->nrma) MatNorm(svd->OP,NORM_INFINITY,&svd->nrma);
389: if (!svd->nrmb) MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb);
390: EPSSetConvergenceTestFunction(cyclic->eps,EPSConv_Cyclic,svd,NULL);
391: } else {
392: EPSSetConvergenceTest(cyclic->eps,EPS_CONV_NORM);break;
393: }
394: break;
395: case SVD_CONV_MAXIT:
396: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Maxit convergence test not supported in this solver");
397: case SVD_CONV_USER:
398: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
399: }
400: }
401: SVDCheckUnsupported(svd,SVD_FEATURE_STOPPING);
402: /* Transfer the trackall option from svd to eps */
403: SVDGetTrackAll(svd,&trackall);
404: EPSSetTrackAll(cyclic->eps,trackall);
405: /* Transfer the initial subspace from svd to eps */
406: if (svd->nini<0 || svd->ninil<0) {
407: for (i=0;i<-PetscMin(svd->nini,svd->ninil);i++) {
408: MatCreateVecs(cyclic->C,&v,NULL);
409: VecGetArrayWrite(v,&va);
410: if (svd->isgeneralized) MatGetLocalSize(svd->B,&p,NULL);
411: k = (svd->isgeneralized && svd->which==SVD_SMALLEST)? p: m; /* size of upper block row */
412: if (i<-svd->ninil) {
413: VecGetArrayRead(svd->ISL[i],&isa);
414: if (svd->isgeneralized) {
415: VecGetLocalSize(svd->ISL[i],&isl);
417: offset = (svd->which==SVD_SMALLEST)? m: 0;
418: PetscArraycpy(va,isa+offset,k);
419: } else {
420: VecGetLocalSize(svd->ISL[i],&isl);
422: PetscArraycpy(va,isa,k);
423: }
424: VecRestoreArrayRead(svd->IS[i],&isa);
425: } else PetscArrayzero(&va,k);
426: if (i<-svd->nini) {
427: VecGetLocalSize(svd->IS[i],&isl);
429: VecGetArrayRead(svd->IS[i],&isa);
430: PetscArraycpy(va+k,isa,n);
431: VecRestoreArrayRead(svd->IS[i],&isa);
432: } else PetscArrayzero(va+k,n);
433: VecRestoreArrayWrite(v,&va);
434: VecDestroy(&svd->IS[i]);
435: svd->IS[i] = v;
436: }
437: svd->nini = PetscMin(svd->nini,svd->ninil);
438: EPSSetInitialSpace(cyclic->eps,-svd->nini,svd->IS);
439: SlepcBasisDestroy_Private(&svd->nini,&svd->IS);
440: SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL);
441: }
442: EPSSetUp(cyclic->eps);
443: EPSGetDimensions(cyclic->eps,NULL,&svd->ncv,&svd->mpd);
444: svd->ncv = PetscMin(svd->ncv,PetscMin(M,N));
445: EPSGetTolerances(cyclic->eps,NULL,&svd->max_it);
446: if (svd->tol==PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;
448: svd->leftbasis = PETSC_TRUE;
449: SVDAllocateSolution(svd,0);
450: PetscFunctionReturn(0);
451: }
453: PetscErrorCode SVDSolve_Cyclic(SVD svd)
454: {
455: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
456: PetscInt i,j,nconv;
457: PetscScalar sigma;
459: EPSSolve(cyclic->eps);
460: EPSGetConverged(cyclic->eps,&nconv);
461: EPSGetIterationNumber(cyclic->eps,&svd->its);
462: EPSGetConvergedReason(cyclic->eps,(EPSConvergedReason*)&svd->reason);
463: for (i=0,j=0;i<nconv;i++) {
464: EPSGetEigenvalue(cyclic->eps,i,&sigma,NULL);
465: if (PetscRealPart(sigma) > 0.0) {
466: if (svd->isgeneralized && svd->which==SVD_SMALLEST) svd->sigma[j] = 1.0/PetscRealPart(sigma);
467: else svd->sigma[j] = PetscRealPart(sigma);
468: j++;
469: }
470: }
471: svd->nconv = j;
472: PetscFunctionReturn(0);
473: }
475: PetscErrorCode SVDComputeVectors_Cyclic(SVD svd)
476: {
477: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
478: PetscInt i,j,m,p,nconv;
479: PetscScalar *dst,sigma;
480: const PetscScalar *src,*px;
481: Vec u,v,x,x1,x2,uv;
483: EPSGetConverged(cyclic->eps,&nconv);
484: MatCreateVecs(cyclic->C,&x,NULL);
485: MatGetLocalSize(svd->A,&m,NULL);
486: if (svd->isgeneralized && svd->which==SVD_SMALLEST) MatCreateVecsEmpty(svd->B,&x1,&x2);
487: else MatCreateVecsEmpty(svd->A,&x2,&x1);
488: if (svd->isgeneralized) {
489: MatCreateVecs(svd->A,NULL,&u);
490: MatCreateVecs(svd->B,NULL,&v);
491: MatGetLocalSize(svd->B,&p,NULL);
492: }
493: for (i=0,j=0;i<nconv;i++) {
494: EPSGetEigenpair(cyclic->eps,i,&sigma,NULL,x,NULL);
495: if (PetscRealPart(sigma) > 0.0) {
496: if (svd->isgeneralized) {
497: if (svd->which==SVD_SMALLEST) {
498: /* evec_i = 1/sqrt(2)*[ v_i; w_i ], w_i = x_i/c_i */
499: VecGetArrayRead(x,&px);
500: VecPlaceArray(x2,px);
501: VecPlaceArray(x1,px+p);
502: VecCopy(x2,v);
503: VecScale(v,PETSC_SQRT2); /* v_i = sqrt(2)*evec_i_1 */
504: VecScale(x1,PETSC_SQRT2); /* w_i = sqrt(2)*evec_i_2 */
505: MatMult(svd->A,x1,u); /* A*w_i = u_i */
506: VecScale(x1,1.0/PetscSqrtScalar(1.0+sigma*sigma)); /* x_i = w_i*c_i */
507: BVInsertVec(svd->V,j,x1);
508: VecResetArray(x2);
509: VecResetArray(x1);
510: VecRestoreArrayRead(x,&px);
511: } else {
512: /* evec_i = 1/sqrt(2)*[ u_i; w_i ], w_i = x_i/s_i */
513: VecGetArrayRead(x,&px);
514: VecPlaceArray(x1,px);
515: VecPlaceArray(x2,px+m);
516: VecCopy(x1,u);
517: VecScale(u,PETSC_SQRT2); /* u_i = sqrt(2)*evec_i_1 */
518: VecScale(x2,PETSC_SQRT2); /* w_i = sqrt(2)*evec_i_2 */
519: MatMult(svd->B,x2,v); /* B*w_i = v_i */
520: VecScale(x2,1.0/PetscSqrtScalar(1.0+sigma*sigma)); /* x_i = w_i*s_i */
521: BVInsertVec(svd->V,j,x2);
522: VecResetArray(x1);
523: VecResetArray(x2);
524: VecRestoreArrayRead(x,&px);
525: }
526: /* copy [u;v] to U[j] */
527: BVGetColumn(svd->U,j,&uv);
528: VecGetArrayWrite(uv,&dst);
529: VecGetArrayRead(u,&src);
530: PetscArraycpy(dst,src,m);
531: VecRestoreArrayRead(u,&src);
532: VecGetArrayRead(v,&src);
533: PetscArraycpy(dst+m,src,p);
534: VecRestoreArrayRead(v,&src);
535: VecRestoreArrayWrite(uv,&dst);
536: BVRestoreColumn(svd->U,j,&uv);
537: } else {
538: VecGetArrayRead(x,&px);
539: VecPlaceArray(x1,px);
540: VecPlaceArray(x2,px+m);
541: BVInsertVec(svd->U,j,x1);
542: BVScaleColumn(svd->U,j,PETSC_SQRT2);
543: BVInsertVec(svd->V,j,x2);
544: BVScaleColumn(svd->V,j,PETSC_SQRT2);
545: VecResetArray(x1);
546: VecResetArray(x2);
547: VecRestoreArrayRead(x,&px);
548: }
549: j++;
550: }
551: }
552: VecDestroy(&x);
553: VecDestroy(&x1);
554: VecDestroy(&x2);
555: if (svd->isgeneralized) {
556: VecDestroy(&u);
557: VecDestroy(&v);
558: }
559: PetscFunctionReturn(0);
560: }
562: static PetscErrorCode EPSMonitor_Cyclic(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
563: {
564: PetscInt i,j;
565: SVD svd = (SVD)ctx;
566: PetscScalar er,ei;
568: nconv = 0;
569: for (i=0,j=0;i<PetscMin(nest,svd->ncv);i++) {
570: er = eigr[i]; ei = eigi[i];
571: STBackTransform(eps->st,1,&er,&ei);
572: if (PetscRealPart(er) > 0.0) {
573: svd->sigma[j] = PetscRealPart(er);
574: svd->errest[j] = errest[i];
575: if (errest[i] && errest[i] < svd->tol) nconv++;
576: j++;
577: }
578: }
579: nest = j;
580: SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest);
581: PetscFunctionReturn(0);
582: }
584: PetscErrorCode SVDSetFromOptions_Cyclic(PetscOptionItems *PetscOptionsObject,SVD svd)
585: {
586: PetscBool set,val;
587: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
588: ST st;
590: PetscOptionsHead(PetscOptionsObject,"SVD Cyclic Options");
592: PetscOptionsBool("-svd_cyclic_explicitmatrix","Use cyclic explicit matrix","SVDCyclicSetExplicitMatrix",cyclic->explicitmatrix,&val,&set);
593: if (set) SVDCyclicSetExplicitMatrix(svd,val);
595: PetscOptionsTail();
597: if (!cyclic->eps) SVDCyclicGetEPS(svd,&cyclic->eps);
598: if (!cyclic->explicitmatrix && !cyclic->usereps) {
599: /* use as default an ST with shell matrix and Jacobi */
600: EPSGetST(cyclic->eps,&st);
601: STSetMatMode(st,ST_MATMODE_SHELL);
602: }
603: EPSSetFromOptions(cyclic->eps);
604: PetscFunctionReturn(0);
605: }
607: static PetscErrorCode SVDCyclicSetExplicitMatrix_Cyclic(SVD svd,PetscBool explicitmat)
608: {
609: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
611: if (cyclic->explicitmatrix != explicitmat) {
612: cyclic->explicitmatrix = explicitmat;
613: svd->state = SVD_STATE_INITIAL;
614: }
615: PetscFunctionReturn(0);
616: }
618: /*@
619: SVDCyclicSetExplicitMatrix - Indicate if the eigensolver operator
620: H(A) = [ 0 A ; A^T 0 ] must be computed explicitly.
622: Logically Collective on svd
624: Input Parameters:
625: + svd - singular value solver
626: - explicitmat - boolean flag indicating if H(A) is built explicitly
628: Options Database Key:
629: . -svd_cyclic_explicitmatrix <boolean> - Indicates the boolean flag
631: Level: advanced
633: .seealso: SVDCyclicGetExplicitMatrix()
634: @*/
635: PetscErrorCode SVDCyclicSetExplicitMatrix(SVD svd,PetscBool explicitmat)
636: {
639: PetscTryMethod(svd,"SVDCyclicSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
640: PetscFunctionReturn(0);
641: }
643: static PetscErrorCode SVDCyclicGetExplicitMatrix_Cyclic(SVD svd,PetscBool *explicitmat)
644: {
645: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
647: *explicitmat = cyclic->explicitmatrix;
648: PetscFunctionReturn(0);
649: }
651: /*@
652: SVDCyclicGetExplicitMatrix - Returns the flag indicating if H(A) is built explicitly.
654: Not Collective
656: Input Parameter:
657: . svd - singular value solver
659: Output Parameter:
660: . explicitmat - the mode flag
662: Level: advanced
664: .seealso: SVDCyclicSetExplicitMatrix()
665: @*/
666: PetscErrorCode SVDCyclicGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
667: {
670: PetscUseMethod(svd,"SVDCyclicGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
671: PetscFunctionReturn(0);
672: }
674: static PetscErrorCode SVDCyclicSetEPS_Cyclic(SVD svd,EPS eps)
675: {
676: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
678: PetscObjectReference((PetscObject)eps);
679: EPSDestroy(&cyclic->eps);
680: cyclic->eps = eps;
681: cyclic->usereps = PETSC_TRUE;
682: PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->eps);
683: svd->state = SVD_STATE_INITIAL;
684: PetscFunctionReturn(0);
685: }
687: /*@
688: SVDCyclicSetEPS - Associate an eigensolver object (EPS) to the
689: singular value solver.
691: Collective on svd
693: Input Parameters:
694: + svd - singular value solver
695: - eps - the eigensolver object
697: Level: advanced
699: .seealso: SVDCyclicGetEPS()
700: @*/
701: PetscErrorCode SVDCyclicSetEPS(SVD svd,EPS eps)
702: {
706: PetscTryMethod(svd,"SVDCyclicSetEPS_C",(SVD,EPS),(svd,eps));
707: PetscFunctionReturn(0);
708: }
710: static PetscErrorCode SVDCyclicGetEPS_Cyclic(SVD svd,EPS *eps)
711: {
712: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
714: if (!cyclic->eps) {
715: EPSCreate(PetscObjectComm((PetscObject)svd),&cyclic->eps);
716: PetscObjectIncrementTabLevel((PetscObject)cyclic->eps,(PetscObject)svd,1);
717: EPSSetOptionsPrefix(cyclic->eps,((PetscObject)svd)->prefix);
718: EPSAppendOptionsPrefix(cyclic->eps,"svd_cyclic_");
719: PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->eps);
720: PetscObjectSetOptions((PetscObject)cyclic->eps,((PetscObject)svd)->options);
721: EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL);
722: EPSMonitorSet(cyclic->eps,EPSMonitor_Cyclic,svd,NULL);
723: }
724: *eps = cyclic->eps;
725: PetscFunctionReturn(0);
726: }
728: /*@
729: SVDCyclicGetEPS - Retrieve the eigensolver object (EPS) associated
730: to the singular value solver.
732: Not Collective
734: Input Parameter:
735: . svd - singular value solver
737: Output Parameter:
738: . eps - the eigensolver object
740: Level: advanced
742: .seealso: SVDCyclicSetEPS()
743: @*/
744: PetscErrorCode SVDCyclicGetEPS(SVD svd,EPS *eps)
745: {
748: PetscUseMethod(svd,"SVDCyclicGetEPS_C",(SVD,EPS*),(svd,eps));
749: PetscFunctionReturn(0);
750: }
752: PetscErrorCode SVDView_Cyclic(SVD svd,PetscViewer viewer)
753: {
754: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
755: PetscBool isascii;
757: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
758: if (isascii) {
759: if (!cyclic->eps) SVDCyclicGetEPS(svd,&cyclic->eps);
760: PetscViewerASCIIPrintf(viewer," %s matrix\n",cyclic->explicitmatrix?"explicit":"implicit");
761: PetscViewerASCIIPushTab(viewer);
762: EPSView(cyclic->eps,viewer);
763: PetscViewerASCIIPopTab(viewer);
764: }
765: PetscFunctionReturn(0);
766: }
768: PetscErrorCode SVDReset_Cyclic(SVD svd)
769: {
770: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
772: EPSReset(cyclic->eps);
773: MatDestroy(&cyclic->C);
774: MatDestroy(&cyclic->D);
775: PetscFunctionReturn(0);
776: }
778: PetscErrorCode SVDDestroy_Cyclic(SVD svd)
779: {
780: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
782: EPSDestroy(&cyclic->eps);
783: PetscFree(svd->data);
784: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",NULL);
785: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",NULL);
786: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",NULL);
787: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",NULL);
788: PetscFunctionReturn(0);
789: }
791: SLEPC_EXTERN PetscErrorCode SVDCreate_Cyclic(SVD svd)
792: {
793: SVD_CYCLIC *cyclic;
795: PetscNewLog(svd,&cyclic);
796: svd->data = (void*)cyclic;
797: svd->ops->solve = SVDSolve_Cyclic;
798: svd->ops->solveg = SVDSolve_Cyclic;
799: svd->ops->setup = SVDSetUp_Cyclic;
800: svd->ops->setfromoptions = SVDSetFromOptions_Cyclic;
801: svd->ops->destroy = SVDDestroy_Cyclic;
802: svd->ops->reset = SVDReset_Cyclic;
803: svd->ops->view = SVDView_Cyclic;
804: svd->ops->computevectors = SVDComputeVectors_Cyclic;
805: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",SVDCyclicSetEPS_Cyclic);
806: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",SVDCyclicGetEPS_Cyclic);
807: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",SVDCyclicSetExplicitMatrix_Cyclic);
808: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",SVDCyclicGetExplicitMatrix_Cyclic);
809: PetscFunctionReturn(0);
810: }