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