Actual source code: cross.c
slepc-main 2024-12-17
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: "cross"
13: Method: Uses a Hermitian eigensolver for A^T*A
14: */
16: #include <slepc/private/svdimpl.h>
18: typedef struct {
19: PetscBool explicitmatrix;
20: EPS eps;
21: PetscBool usereps;
22: Mat C,D;
23: } SVD_CROSS;
25: typedef struct {
26: Mat A,AT;
27: Vec w,diag,omega;
28: PetscBool swapped;
29: } SVD_CROSS_SHELL;
31: static PetscErrorCode MatMult_Cross(Mat B,Vec x,Vec y)
32: {
33: SVD_CROSS_SHELL *ctx;
35: PetscFunctionBegin;
36: PetscCall(MatShellGetContext(B,&ctx));
37: PetscCall(MatMult(ctx->A,x,ctx->w));
38: if (ctx->omega && !ctx->swapped) PetscCall(VecPointwiseMult(ctx->w,ctx->w,ctx->omega));
39: PetscCall(MatMult(ctx->AT,ctx->w,y));
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: static PetscErrorCode MatGetDiagonal_Cross(Mat B,Vec d)
44: {
45: SVD_CROSS_SHELL *ctx;
46: PetscMPIInt len;
47: PetscInt N,n,i,j,start,end,ncols;
48: PetscScalar *work1,*work2,*diag;
49: const PetscInt *cols;
50: const PetscScalar *vals;
52: PetscFunctionBegin;
53: PetscCall(MatShellGetContext(B,&ctx));
54: if (!ctx->diag) {
55: /* compute diagonal from rows and store in ctx->diag */
56: PetscCall(VecDuplicate(d,&ctx->diag));
57: PetscCall(MatGetSize(ctx->A,NULL,&N));
58: PetscCall(MatGetLocalSize(ctx->A,NULL,&n));
59: PetscCall(PetscCalloc2(N,&work1,N,&work2));
60: if (ctx->swapped) {
61: PetscCall(MatGetOwnershipRange(ctx->AT,&start,&end));
62: for (i=start;i<end;i++) {
63: PetscCall(MatGetRow(ctx->AT,i,&ncols,NULL,&vals));
64: for (j=0;j<ncols;j++) work1[i] += vals[j]*vals[j];
65: PetscCall(MatRestoreRow(ctx->AT,i,&ncols,NULL,&vals));
66: }
67: } else {
68: PetscCall(MatGetOwnershipRange(ctx->A,&start,&end));
69: for (i=start;i<end;i++) {
70: PetscCall(MatGetRow(ctx->A,i,&ncols,&cols,&vals));
71: for (j=0;j<ncols;j++) work1[cols[j]] += vals[j]*vals[j];
72: PetscCall(MatRestoreRow(ctx->A,i,&ncols,&cols,&vals));
73: }
74: }
75: PetscCall(PetscMPIIntCast(N,&len));
76: PetscCallMPI(MPIU_Allreduce(work1,work2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)B)));
77: PetscCall(VecGetOwnershipRange(ctx->diag,&start,&end));
78: PetscCall(VecGetArrayWrite(ctx->diag,&diag));
79: for (i=start;i<end;i++) diag[i-start] = work2[i];
80: PetscCall(VecRestoreArrayWrite(ctx->diag,&diag));
81: PetscCall(PetscFree2(work1,work2));
82: }
83: PetscCall(VecCopy(ctx->diag,d));
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: static PetscErrorCode MatDestroy_Cross(Mat B)
88: {
89: SVD_CROSS_SHELL *ctx;
91: PetscFunctionBegin;
92: PetscCall(MatShellGetContext(B,&ctx));
93: PetscCall(VecDestroy(&ctx->w));
94: PetscCall(VecDestroy(&ctx->diag));
95: PetscCall(PetscFree(ctx));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: static PetscErrorCode SVDCrossGetProductMat(SVD svd,Mat A,Mat AT,Mat *C)
100: {
101: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
102: SVD_CROSS_SHELL *ctx;
103: PetscInt n;
104: VecType vtype;
105: Mat B;
107: PetscFunctionBegin;
108: if (cross->explicitmatrix) {
109: if (!svd->ishyperbolic || svd->swapped) B = (!svd->expltrans && svd->swapped)? AT: A;
110: else { /* duplicate A and scale by signature */
111: PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&B));
112: PetscCall(MatDiagonalScale(B,svd->omega,NULL));
113: }
114: if (svd->expltrans) { /* explicit transpose */
115: PetscCall(MatProductCreate(AT,B,NULL,C));
116: PetscCall(MatProductSetType(*C,MATPRODUCT_AB));
117: } else { /* implicit transpose */
118: #if defined(PETSC_USE_COMPLEX)
119: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Must use explicit transpose with complex scalars");
120: #else
121: if (!svd->swapped) {
122: PetscCall(MatProductCreate(A,B,NULL,C));
123: PetscCall(MatProductSetType(*C,MATPRODUCT_AtB));
124: } else {
125: PetscCall(MatProductCreate(B,AT,NULL,C));
126: PetscCall(MatProductSetType(*C,MATPRODUCT_ABt));
127: }
128: #endif
129: }
130: PetscCall(MatProductSetFromOptions(*C));
131: PetscCall(MatProductSymbolic(*C));
132: PetscCall(MatProductNumeric(*C));
133: if (svd->ishyperbolic && !svd->swapped) PetscCall(MatDestroy(&B));
134: } else {
135: PetscCall(PetscNew(&ctx));
136: ctx->A = A;
137: ctx->AT = AT;
138: ctx->omega = svd->omega;
139: ctx->swapped = svd->swapped;
140: PetscCall(MatCreateVecs(A,NULL,&ctx->w));
141: PetscCall(MatGetLocalSize(A,NULL,&n));
142: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)svd),n,n,PETSC_DETERMINE,PETSC_DETERMINE,(void*)ctx,C));
143: PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cross));
144: if (!svd->ishyperbolic || svd->swapped) PetscCall(MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Cross));
145: PetscCall(MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_Cross));
146: PetscCall(MatGetVecType(A,&vtype));
147: PetscCall(MatSetVecType(*C,vtype));
148: }
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: /* Convergence test relative to the norm of R (used in GSVD only) */
153: static PetscErrorCode EPSConv_Cross(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
154: {
155: SVD svd = (SVD)ctx;
157: PetscFunctionBegin;
158: *errest = res/PetscMax(svd->nrma,svd->nrmb);
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: static PetscErrorCode SVDSetUp_Cross(SVD svd)
163: {
164: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
165: ST st;
166: PetscBool trackall,issinv,isks;
167: EPSProblemType ptype;
168: EPSWhich which;
169: Mat Omega;
170: MatType Atype;
171: PetscInt n,N;
173: PetscFunctionBegin;
174: if (svd->nsv==0 && svd->stop!=SVD_STOP_THRESHOLD) svd->nsv = 1;
175: if (!cross->eps) PetscCall(SVDCrossGetEPS(svd,&cross->eps));
176: PetscCall(MatDestroy(&cross->C));
177: PetscCall(MatDestroy(&cross->D));
178: PetscCall(SVDCrossGetProductMat(svd,svd->A,svd->AT,&cross->C));
179: if (svd->isgeneralized) {
180: PetscCall(SVDCrossGetProductMat(svd,svd->B,svd->BT,&cross->D));
181: PetscCall(EPSSetOperators(cross->eps,cross->C,cross->D));
182: PetscCall(EPSGetProblemType(cross->eps,&ptype));
183: if (!ptype) PetscCall(EPSSetProblemType(cross->eps,EPS_GHEP));
184: } else if (svd->ishyperbolic && svd->swapped) {
185: PetscCall(MatGetType(svd->OP,&Atype));
186: PetscCall(MatGetSize(svd->A,NULL,&N));
187: PetscCall(MatGetLocalSize(svd->A,NULL,&n));
188: PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Omega));
189: PetscCall(MatSetSizes(Omega,n,n,N,N));
190: PetscCall(MatSetType(Omega,Atype));
191: PetscCall(MatDiagonalSet(Omega,svd->omega,INSERT_VALUES));
192: PetscCall(EPSSetOperators(cross->eps,cross->C,Omega));
193: PetscCall(EPSSetProblemType(cross->eps,EPS_GHIEP));
194: PetscCall(MatDestroy(&Omega));
195: } else {
196: PetscCall(EPSSetOperators(cross->eps,cross->C,NULL));
197: PetscCall(EPSSetProblemType(cross->eps,EPS_HEP));
198: }
199: if (!cross->usereps) {
200: PetscCall(EPSGetST(cross->eps,&st));
201: PetscCall(PetscObjectTypeCompare((PetscObject)st,STSINVERT,&issinv));
202: PetscCall(PetscObjectTypeCompare((PetscObject)cross->eps,EPSKRYLOVSCHUR,&isks));
203: if (svd->isgeneralized && svd->which==SVD_SMALLEST) {
204: if (cross->explicitmatrix && isks && !issinv) { /* default to shift-and-invert */
205: PetscCall(STSetType(st,STSINVERT));
206: PetscCall(EPSSetTarget(cross->eps,0.0));
207: which = EPS_TARGET_REAL;
208: } else which = issinv?EPS_TARGET_REAL:EPS_SMALLEST_REAL;
209: } else {
210: if (issinv) which = EPS_TARGET_MAGNITUDE;
211: else if (svd->ishyperbolic) which = svd->which==SVD_LARGEST?EPS_LARGEST_MAGNITUDE:EPS_SMALLEST_MAGNITUDE;
212: else which = svd->which==SVD_LARGEST?EPS_LARGEST_MAGNITUDE:EPS_SMALLEST_MAGNITUDE;
213: }
214: PetscCall(EPSSetWhichEigenpairs(cross->eps,which));
215: PetscCall(EPSSetDimensions(cross->eps,svd->nsv?svd->nsv:PETSC_CURRENT,svd->ncv,svd->mpd));
216: if (svd->stop==SVD_STOP_THRESHOLD) PetscCall(EPSSetThreshold(cross->eps,svd->thres*svd->thres,svd->threlative));
217: PetscCall(EPSSetTolerances(cross->eps,svd->tol==(PetscReal)PETSC_DETERMINE?SLEPC_DEFAULT_TOL/10.0:svd->tol,svd->max_it));
218: switch (svd->conv) {
219: case SVD_CONV_ABS:
220: PetscCall(EPSSetConvergenceTest(cross->eps,EPS_CONV_ABS));break;
221: case SVD_CONV_REL:
222: PetscCall(EPSSetConvergenceTest(cross->eps,EPS_CONV_REL));break;
223: case SVD_CONV_NORM:
224: if (svd->isgeneralized) {
225: if (!svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
226: if (!svd->nrmb) PetscCall(MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb));
227: PetscCall(EPSSetConvergenceTestFunction(cross->eps,EPSConv_Cross,svd,NULL));
228: } else {
229: PetscCall(EPSSetConvergenceTest(cross->eps,EPS_CONV_NORM));break;
230: }
231: break;
232: case SVD_CONV_MAXIT:
233: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Maxit convergence test not supported in this solver");
234: case SVD_CONV_USER:
235: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
236: }
237: }
238: /* Transfer the trackall option from svd to eps */
239: PetscCall(SVDGetTrackAll(svd,&trackall));
240: PetscCall(EPSSetTrackAll(cross->eps,trackall));
241: /* Transfer the initial space from svd to eps */
242: if (svd->nini<0) {
243: PetscCall(EPSSetInitialSpace(cross->eps,-svd->nini,svd->IS));
244: PetscCall(SlepcBasisDestroy_Private(&svd->nini,&svd->IS));
245: }
246: PetscCall(EPSSetUp(cross->eps));
247: PetscCall(EPSGetDimensions(cross->eps,NULL,&svd->ncv,&svd->mpd));
248: PetscCall(EPSGetTolerances(cross->eps,NULL,&svd->max_it));
249: if (svd->tol==(PetscReal)PETSC_DETERMINE) svd->tol = SLEPC_DEFAULT_TOL;
251: svd->leftbasis = PETSC_FALSE;
252: PetscCall(SVDAllocateSolution(svd,0));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: static PetscErrorCode SVDSolve_Cross(SVD svd)
257: {
258: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
259: PetscInt i;
260: PetscScalar lambda;
261: PetscReal sigma;
263: PetscFunctionBegin;
264: PetscCall(EPSSolve(cross->eps));
265: PetscCall(EPSGetConverged(cross->eps,&svd->nconv));
266: PetscCall(EPSGetIterationNumber(cross->eps,&svd->its));
267: PetscCall(EPSGetConvergedReason(cross->eps,(EPSConvergedReason*)&svd->reason));
268: for (i=0;i<svd->nconv;i++) {
269: PetscCall(EPSGetEigenvalue(cross->eps,i,&lambda,NULL));
270: sigma = PetscRealPart(lambda);
271: if (svd->ishyperbolic) svd->sigma[i] = PetscSqrtReal(PetscAbsReal(sigma));
272: else {
273: PetscCheck(sigma>-10*PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)svd),PETSC_ERR_FP,"Negative eigenvalue computed by EPS: %g",(double)sigma);
274: if (sigma<0.0) {
275: PetscCall(PetscInfo(svd,"Negative eigenvalue computed by EPS: %g, resetting to 0\n",(double)sigma));
276: sigma = 0.0;
277: }
278: svd->sigma[i] = PetscSqrtReal(sigma);
279: }
280: }
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: static PetscErrorCode SVDComputeVectors_Cross(SVD svd)
285: {
286: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
287: PetscInt i,mloc,ploc;
288: Vec u,v,x,uv,w,omega2=NULL;
289: Mat Omega;
290: PetscScalar *dst,alpha,lambda,*varray;
291: const PetscScalar *src;
292: PetscReal nrm;
294: PetscFunctionBegin;
295: if (svd->isgeneralized) {
296: PetscCall(MatCreateVecs(svd->A,NULL,&u));
297: PetscCall(VecGetLocalSize(u,&mloc));
298: PetscCall(MatCreateVecs(svd->B,NULL,&v));
299: PetscCall(VecGetLocalSize(v,&ploc));
300: for (i=0;i<svd->nconv;i++) {
301: PetscCall(BVGetColumn(svd->V,i,&x));
302: PetscCall(EPSGetEigenpair(cross->eps,i,&lambda,NULL,x,NULL));
303: PetscCall(MatMult(svd->A,x,u)); /* u_i*c_i/alpha = A*x_i */
304: PetscCall(VecNormalize(u,NULL));
305: PetscCall(MatMult(svd->B,x,v)); /* v_i*s_i/alpha = B*x_i */
306: PetscCall(VecNormalize(v,&nrm)); /* ||v||_2 = s_i/alpha */
307: alpha = 1.0/(PetscSqrtReal(1.0+PetscRealPart(lambda))*nrm); /* alpha=s_i/||v||_2 */
308: PetscCall(VecScale(x,alpha));
309: PetscCall(BVRestoreColumn(svd->V,i,&x));
310: /* copy [u;v] to U[i] */
311: PetscCall(BVGetColumn(svd->U,i,&uv));
312: PetscCall(VecGetArrayWrite(uv,&dst));
313: PetscCall(VecGetArrayRead(u,&src));
314: PetscCall(PetscArraycpy(dst,src,mloc));
315: PetscCall(VecRestoreArrayRead(u,&src));
316: PetscCall(VecGetArrayRead(v,&src));
317: PetscCall(PetscArraycpy(dst+mloc,src,ploc));
318: PetscCall(VecRestoreArrayRead(v,&src));
319: PetscCall(VecRestoreArrayWrite(uv,&dst));
320: PetscCall(BVRestoreColumn(svd->U,i,&uv));
321: }
322: PetscCall(VecDestroy(&v));
323: PetscCall(VecDestroy(&u));
324: } else if (svd->ishyperbolic && svd->swapped) { /* was solved as GHIEP, set u=Omega*u and normalize */
325: PetscCall(EPSGetOperators(cross->eps,NULL,&Omega));
326: PetscCall(MatCreateVecs(Omega,&w,NULL));
327: PetscCall(VecCreateSeq(PETSC_COMM_SELF,svd->ncv,&omega2));
328: PetscCall(VecGetArrayWrite(omega2,&varray));
329: for (i=0;i<svd->nconv;i++) {
330: PetscCall(BVGetColumn(svd->V,i,&v));
331: PetscCall(EPSGetEigenvector(cross->eps,i,v,NULL));
332: PetscCall(MatMult(Omega,v,w));
333: PetscCall(VecDot(v,w,&alpha));
334: svd->sign[i] = PetscSign(PetscRealPart(alpha));
335: varray[i] = svd->sign[i];
336: alpha = 1.0/PetscSqrtScalar(PetscAbsScalar(alpha));
337: PetscCall(VecScale(w,alpha));
338: PetscCall(VecCopy(w,v));
339: PetscCall(BVRestoreColumn(svd->V,i,&v));
340: }
341: PetscCall(BVSetSignature(svd->V,omega2));
342: PetscCall(VecRestoreArrayWrite(omega2,&varray));
343: PetscCall(VecDestroy(&omega2));
344: PetscCall(VecDestroy(&w));
345: PetscCall(SVDComputeVectors_Left(svd));
346: } else {
347: for (i=0;i<svd->nconv;i++) {
348: PetscCall(BVGetColumn(svd->V,i,&v));
349: PetscCall(EPSGetEigenvector(cross->eps,i,v,NULL));
350: PetscCall(BVRestoreColumn(svd->V,i,&v));
351: }
352: PetscCall(SVDComputeVectors_Left(svd));
353: }
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: static PetscErrorCode EPSMonitor_Cross(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
358: {
359: PetscInt i,ncv;
360: SVD svd = (SVD)ctx;
361: SVD_CROSS *cross;
362: PetscScalar er,ei;
363: ST st;
365: PetscFunctionBegin;
366: if (svd->stop==SVD_STOP_THRESHOLD) {
367: cross = (SVD_CROSS*)svd->data;
368: PetscCall(EPSGetDimensions(cross->eps,NULL,&ncv,NULL));
369: if (ncv!=svd->ncv) { /* reallocate */
370: PetscCall(SVDReallocateSolution(svd,ncv));
371: for (i=svd->ncv;i<ncv;i++) svd->perm[i] = i;
372: svd->ncv = ncv;
373: }
374: }
375: PetscCall(EPSGetST(eps,&st));
376: for (i=0;i<PetscMin(nest,svd->ncv);i++) {
377: er = eigr[i]; ei = eigi[i];
378: PetscCall(STBackTransform(st,1,&er,&ei));
379: svd->sigma[i] = PetscSqrtReal(PetscAbsReal(PetscRealPart(er)));
380: svd->errest[i] = errest[i];
381: }
382: PetscCall(SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest));
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: static PetscErrorCode SVDSetFromOptions_Cross(SVD svd,PetscOptionItems *PetscOptionsObject)
387: {
388: PetscBool set,val;
389: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
390: ST st;
392: PetscFunctionBegin;
393: PetscOptionsHeadBegin(PetscOptionsObject,"SVD Cross Options");
395: PetscCall(PetscOptionsBool("-svd_cross_explicitmatrix","Use cross explicit matrix","SVDCrossSetExplicitMatrix",cross->explicitmatrix,&val,&set));
396: if (set) PetscCall(SVDCrossSetExplicitMatrix(svd,val));
398: PetscOptionsHeadEnd();
400: if (!cross->eps) PetscCall(SVDCrossGetEPS(svd,&cross->eps));
401: if (!cross->explicitmatrix && !cross->usereps) {
402: /* use as default an ST with shell matrix and Jacobi */
403: PetscCall(EPSGetST(cross->eps,&st));
404: PetscCall(STSetMatMode(st,ST_MATMODE_SHELL));
405: }
406: PetscCall(EPSSetFromOptions(cross->eps));
407: PetscFunctionReturn(PETSC_SUCCESS);
408: }
410: static PetscErrorCode SVDCrossSetExplicitMatrix_Cross(SVD svd,PetscBool explicitmatrix)
411: {
412: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
414: PetscFunctionBegin;
415: if (cross->explicitmatrix != explicitmatrix) {
416: cross->explicitmatrix = explicitmatrix;
417: svd->state = SVD_STATE_INITIAL;
418: }
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: /*@
423: SVDCrossSetExplicitMatrix - Indicate if the eigensolver operator A^T*A must
424: be computed explicitly.
426: Logically Collective
428: Input Parameters:
429: + svd - singular value solver
430: - explicitmat - boolean flag indicating if A^T*A is built explicitly
432: Options Database Key:
433: . -svd_cross_explicitmatrix <boolean> - Indicates the boolean flag
435: Level: advanced
437: .seealso: SVDCrossGetExplicitMatrix()
438: @*/
439: PetscErrorCode SVDCrossSetExplicitMatrix(SVD svd,PetscBool explicitmat)
440: {
441: PetscFunctionBegin;
444: PetscTryMethod(svd,"SVDCrossSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: static PetscErrorCode SVDCrossGetExplicitMatrix_Cross(SVD svd,PetscBool *explicitmat)
449: {
450: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
452: PetscFunctionBegin;
453: *explicitmat = cross->explicitmatrix;
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: /*@
458: SVDCrossGetExplicitMatrix - Returns the flag indicating if A^T*A is built explicitly.
460: Not Collective
462: Input Parameter:
463: . svd - singular value solver
465: Output Parameter:
466: . explicitmat - the mode flag
468: Level: advanced
470: .seealso: SVDCrossSetExplicitMatrix()
471: @*/
472: PetscErrorCode SVDCrossGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
473: {
474: PetscFunctionBegin;
476: PetscAssertPointer(explicitmat,2);
477: PetscUseMethod(svd,"SVDCrossGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: static PetscErrorCode SVDCrossSetEPS_Cross(SVD svd,EPS eps)
482: {
483: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
485: PetscFunctionBegin;
486: PetscCall(PetscObjectReference((PetscObject)eps));
487: PetscCall(EPSDestroy(&cross->eps));
488: cross->eps = eps;
489: cross->usereps = PETSC_TRUE;
490: svd->state = SVD_STATE_INITIAL;
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: /*@
495: SVDCrossSetEPS - Associate an eigensolver object (EPS) to the
496: singular value solver.
498: Collective
500: Input Parameters:
501: + svd - singular value solver
502: - eps - the eigensolver object
504: Level: advanced
506: .seealso: SVDCrossGetEPS()
507: @*/
508: PetscErrorCode SVDCrossSetEPS(SVD svd,EPS eps)
509: {
510: PetscFunctionBegin;
513: PetscCheckSameComm(svd,1,eps,2);
514: PetscTryMethod(svd,"SVDCrossSetEPS_C",(SVD,EPS),(svd,eps));
515: PetscFunctionReturn(PETSC_SUCCESS);
516: }
518: static PetscErrorCode SVDCrossGetEPS_Cross(SVD svd,EPS *eps)
519: {
520: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
522: PetscFunctionBegin;
523: if (!cross->eps) {
524: PetscCall(EPSCreate(PetscObjectComm((PetscObject)svd),&cross->eps));
525: PetscCall(PetscObjectIncrementTabLevel((PetscObject)cross->eps,(PetscObject)svd,1));
526: PetscCall(EPSSetOptionsPrefix(cross->eps,((PetscObject)svd)->prefix));
527: PetscCall(EPSAppendOptionsPrefix(cross->eps,"svd_cross_"));
528: PetscCall(PetscObjectSetOptions((PetscObject)cross->eps,((PetscObject)svd)->options));
529: PetscCall(EPSSetWhichEigenpairs(cross->eps,EPS_LARGEST_MAGNITUDE));
530: PetscCall(EPSMonitorSet(cross->eps,EPSMonitor_Cross,svd,NULL));
531: }
532: *eps = cross->eps;
533: PetscFunctionReturn(PETSC_SUCCESS);
534: }
536: /*@
537: SVDCrossGetEPS - Retrieve the eigensolver object (EPS) associated
538: to the singular value solver.
540: Collective
542: Input Parameter:
543: . svd - singular value solver
545: Output Parameter:
546: . eps - the eigensolver object
548: Level: advanced
550: .seealso: SVDCrossSetEPS()
551: @*/
552: PetscErrorCode SVDCrossGetEPS(SVD svd,EPS *eps)
553: {
554: PetscFunctionBegin;
556: PetscAssertPointer(eps,2);
557: PetscUseMethod(svd,"SVDCrossGetEPS_C",(SVD,EPS*),(svd,eps));
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: static PetscErrorCode SVDView_Cross(SVD svd,PetscViewer viewer)
562: {
563: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
564: PetscBool isascii;
566: PetscFunctionBegin;
567: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
568: if (isascii) {
569: if (!cross->eps) PetscCall(SVDCrossGetEPS(svd,&cross->eps));
570: PetscCall(PetscViewerASCIIPrintf(viewer," %s matrix\n",cross->explicitmatrix?"explicit":"implicit"));
571: PetscCall(PetscViewerASCIIPushTab(viewer));
572: PetscCall(EPSView(cross->eps,viewer));
573: PetscCall(PetscViewerASCIIPopTab(viewer));
574: }
575: PetscFunctionReturn(PETSC_SUCCESS);
576: }
578: static PetscErrorCode SVDReset_Cross(SVD svd)
579: {
580: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
582: PetscFunctionBegin;
583: PetscCall(EPSReset(cross->eps));
584: PetscCall(MatDestroy(&cross->C));
585: PetscCall(MatDestroy(&cross->D));
586: PetscFunctionReturn(PETSC_SUCCESS);
587: }
589: static PetscErrorCode SVDDestroy_Cross(SVD svd)
590: {
591: SVD_CROSS *cross = (SVD_CROSS*)svd->data;
593: PetscFunctionBegin;
594: PetscCall(EPSDestroy(&cross->eps));
595: PetscCall(PetscFree(svd->data));
596: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",NULL));
597: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",NULL));
598: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",NULL));
599: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",NULL));
600: PetscFunctionReturn(PETSC_SUCCESS);
601: }
603: SLEPC_EXTERN PetscErrorCode SVDCreate_Cross(SVD svd)
604: {
605: SVD_CROSS *cross;
607: PetscFunctionBegin;
608: PetscCall(PetscNew(&cross));
609: svd->data = (void*)cross;
611: svd->ops->solve = SVDSolve_Cross;
612: svd->ops->solveg = SVDSolve_Cross;
613: svd->ops->solveh = SVDSolve_Cross;
614: svd->ops->setup = SVDSetUp_Cross;
615: svd->ops->setfromoptions = SVDSetFromOptions_Cross;
616: svd->ops->destroy = SVDDestroy_Cross;
617: svd->ops->reset = SVDReset_Cross;
618: svd->ops->view = SVDView_Cross;
619: svd->ops->computevectors = SVDComputeVectors_Cross;
620: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",SVDCrossSetEPS_Cross));
621: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",SVDCrossGetEPS_Cross));
622: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",SVDCrossSetExplicitMatrix_Cross));
623: PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",SVDCrossGetExplicitMatrix_Cross));
624: PetscFunctionReturn(PETSC_SUCCESS);
625: }