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