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