Actual source code: cross.c

slepc-3.18.0 2022-10-01
Report Typos and Errors
  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>
 17: #include <slepc/private/epsimpl.h>

 19: typedef struct {
 20:   PetscBool explicitmatrix;
 21:   EPS       eps;
 22:   PetscBool usereps;
 23:   Mat       C,D;
 24: } SVD_CROSS;

 26: typedef struct {
 27:   Mat       A,AT;
 28:   Vec       w,diag,omega;
 29:   PetscBool swapped;
 30: } SVD_CROSS_SHELL;

 32: static PetscErrorCode MatMult_Cross(Mat B,Vec x,Vec y)
 33: {
 34:   SVD_CROSS_SHELL *ctx;

 36:   MatShellGetContext(B,&ctx);
 37:   MatMult(ctx->A,x,ctx->w);
 38:   if (ctx->omega && !ctx->swapped) VecPointwiseMult(ctx->w,ctx->w,ctx->omega);
 39:   MatMult(ctx->AT,ctx->w,y);
 40:   return 0;
 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:   MatShellGetContext(B,&ctx);
 53:   if (!ctx->diag) {
 54:     /* compute diagonal from rows and store in ctx->diag */
 55:     VecDuplicate(d,&ctx->diag);
 56:     MatGetSize(ctx->A,NULL,&N);
 57:     MatGetLocalSize(ctx->A,NULL,&n);
 58:     PetscCalloc2(N,&work1,N,&work2);
 59:     if (ctx->swapped) {
 60:       MatGetOwnershipRange(ctx->AT,&start,&end);
 61:       for (i=start;i<end;i++) {
 62:         MatGetRow(ctx->AT,i,&ncols,NULL,&vals);
 63:         for (j=0;j<ncols;j++) work1[i] += vals[j]*vals[j];
 64:         MatRestoreRow(ctx->AT,i,&ncols,NULL,&vals);
 65:       }
 66:     } else {
 67:       MatGetOwnershipRange(ctx->A,&start,&end);
 68:       for (i=start;i<end;i++) {
 69:         MatGetRow(ctx->A,i,&ncols,&cols,&vals);
 70:         for (j=0;j<ncols;j++) work1[cols[j]] += vals[j]*vals[j];
 71:         MatRestoreRow(ctx->A,i,&ncols,&cols,&vals);
 72:       }
 73:     }
 74:     PetscMPIIntCast(N,&len);
 75:     MPIU_Allreduce(work1,work2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)B));
 76:     VecGetOwnershipRange(ctx->diag,&start,&end);
 77:     VecGetArrayWrite(ctx->diag,&diag);
 78:     for (i=start;i<end;i++) diag[i-start] = work2[i];
 79:     VecRestoreArrayWrite(ctx->diag,&diag);
 80:     PetscFree2(work1,work2);
 81:   }
 82:   VecCopy(ctx->diag,d);
 83:   return 0;
 84: }

 86: static PetscErrorCode MatDestroy_Cross(Mat B)
 87: {
 88:   SVD_CROSS_SHELL *ctx;

 90:   MatShellGetContext(B,&ctx);
 91:   VecDestroy(&ctx->w);
 92:   VecDestroy(&ctx->diag);
 93:   PetscFree(ctx);
 94:   return 0;
 95: }

 97: static PetscErrorCode SVDCrossGetProductMat(SVD svd,Mat A,Mat AT,Mat *C)
 98: {
 99:   SVD_CROSS       *cross = (SVD_CROSS*)svd->data;
100:   SVD_CROSS_SHELL *ctx;
101:   PetscInt        n;
102:   VecType         vtype;
103:   Mat             B;

105:   if (cross->explicitmatrix) {
106:     if (!svd->ishyperbolic || svd->swapped) B = (!svd->expltrans && svd->swapped)? AT: A;
107:     else {  /* duplicate A and scale by signature */
108:       MatDuplicate(A,MAT_COPY_VALUES,&B);
109:       MatDiagonalScale(B,svd->omega,NULL);
110:     }
111:     if (svd->expltrans) {  /* explicit transpose */
112:       MatProductCreate(AT,B,NULL,C);
113:       MatProductSetType(*C,MATPRODUCT_AB);
114:     } else {  /* implicit transpose */
115: #if defined(PETSC_USE_COMPLEX)
116:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Must use explicit transpose with complex scalars");
117: #else
118:       if (!svd->swapped) {
119:         MatProductCreate(A,B,NULL,C);
120:         MatProductSetType(*C,MATPRODUCT_AtB);
121:       } else {
122:         MatProductCreate(B,AT,NULL,C);
123:         MatProductSetType(*C,MATPRODUCT_ABt);
124:       }
125: #endif
126:     }
127:     MatProductSetFromOptions(*C);
128:     MatProductSymbolic(*C);
129:     MatProductNumeric(*C);
130:     if (svd->ishyperbolic && !svd->swapped) MatDestroy(&B);
131:   } else {
132:     PetscNew(&ctx);
133:     ctx->A       = A;
134:     ctx->AT      = AT;
135:     ctx->omega   = svd->omega;
136:     ctx->swapped = svd->swapped;
137:     MatCreateVecs(A,NULL,&ctx->w);
138:     MatGetLocalSize(A,NULL,&n);
139:     MatCreateShell(PetscObjectComm((PetscObject)svd),n,n,PETSC_DETERMINE,PETSC_DETERMINE,(void*)ctx,C);
140:     MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cross);
141:     if (!svd->ishyperbolic || svd->swapped) MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Cross);
142:     MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_Cross);
143:     MatGetVecType(A,&vtype);
144:     MatSetVecType(*C,vtype);
145:   }
146:   return 0;
147: }

149: /* Convergence test relative to the norm of R (used in GSVD only) */
150: static PetscErrorCode EPSConv_Cross(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
151: {
152:   SVD svd = (SVD)ctx;

154:   *errest = res/PetscMax(svd->nrma,svd->nrmb);
155:   return 0;
156: }

158: PetscErrorCode SVDSetUp_Cross(SVD svd)
159: {
160:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
161:   ST             st;
162:   PetscBool      trackall,issinv,isks;
163:   EPSProblemType ptype;
164:   EPSWhich       which;
165:   Mat            Omega;
166:   MatType        Atype;
167:   PetscInt       n,N;

169:   if (!cross->eps) SVDCrossGetEPS(svd,&cross->eps);
170:   MatDestroy(&cross->C);
171:   MatDestroy(&cross->D);
172:   SVDCrossGetProductMat(svd,svd->A,svd->AT,&cross->C);
173:   if (svd->isgeneralized) {
174:     SVDCrossGetProductMat(svd,svd->B,svd->BT,&cross->D);
175:     EPSSetOperators(cross->eps,cross->C,cross->D);
176:     EPSGetProblemType(cross->eps,&ptype);
177:     if (!ptype) EPSSetProblemType(cross->eps,EPS_GHEP);
178:   } else if (svd->ishyperbolic && svd->swapped) {
179:     MatGetType(svd->OP,&Atype);
180:     MatGetSize(svd->A,NULL,&N);
181:     MatGetLocalSize(svd->A,NULL,&n);
182:     MatCreate(PetscObjectComm((PetscObject)svd),&Omega);
183:     MatSetSizes(Omega,n,n,N,N);
184:     MatSetType(Omega,Atype);
185:     MatSetUp(Omega);
186:     MatDiagonalSet(Omega,svd->omega,INSERT_VALUES);
187:     EPSSetOperators(cross->eps,cross->C,Omega);
188:     EPSSetProblemType(cross->eps,EPS_GHIEP);
189:     MatDestroy(&Omega);
190:   } else {
191:     EPSSetOperators(cross->eps,cross->C,NULL);
192:     EPSSetProblemType(cross->eps,EPS_HEP);
193:   }
194:   if (!cross->usereps) {
195:     EPSGetST(cross->eps,&st);
196:     PetscObjectTypeCompare((PetscObject)st,STSINVERT,&issinv);
197:     PetscObjectTypeCompare((PetscObject)cross->eps,EPSKRYLOVSCHUR,&isks);
198:     if (svd->isgeneralized && svd->which==SVD_SMALLEST) {
199:       if (cross->explicitmatrix && isks && !issinv) {  /* default to shift-and-invert */
200:         STSetType(st,STSINVERT);
201:         EPSSetTarget(cross->eps,0.0);
202:         which = EPS_TARGET_REAL;
203:       } else which = issinv?EPS_TARGET_REAL:EPS_SMALLEST_REAL;
204:     } else {
205:       if (issinv) which = EPS_TARGET_MAGNITUDE;
206:       else if (svd->ishyperbolic) which = svd->which==SVD_LARGEST?EPS_LARGEST_MAGNITUDE:EPS_SMALLEST_MAGNITUDE;
207:       else which = svd->which==SVD_LARGEST?EPS_LARGEST_REAL:EPS_SMALLEST_REAL;
208:     }
209:     EPSSetWhichEigenpairs(cross->eps,which);
210:     EPSSetDimensions(cross->eps,svd->nsv,svd->ncv,svd->mpd);
211:     EPSSetTolerances(cross->eps,svd->tol==PETSC_DEFAULT?SLEPC_DEFAULT_TOL/10.0:svd->tol,svd->max_it);
212:     switch (svd->conv) {
213:     case SVD_CONV_ABS:
214:       EPSSetConvergenceTest(cross->eps,EPS_CONV_ABS);break;
215:     case SVD_CONV_REL:
216:       EPSSetConvergenceTest(cross->eps,EPS_CONV_REL);break;
217:     case SVD_CONV_NORM:
218:       if (svd->isgeneralized) {
219:         if (!svd->nrma) MatNorm(svd->OP,NORM_INFINITY,&svd->nrma);
220:         if (!svd->nrmb) MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb);
221:         EPSSetConvergenceTestFunction(cross->eps,EPSConv_Cross,svd,NULL);
222:       } else {
223:         EPSSetConvergenceTest(cross->eps,EPS_CONV_NORM);break;
224:       }
225:       break;
226:     case SVD_CONV_MAXIT:
227:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Maxit convergence test not supported in this solver");
228:     case SVD_CONV_USER:
229:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
230:     }
231:   }
232:   SVDCheckUnsupported(svd,SVD_FEATURE_STOPPING);
233:   /* Transfer the trackall option from svd to eps */
234:   SVDGetTrackAll(svd,&trackall);
235:   EPSSetTrackAll(cross->eps,trackall);
236:   /* Transfer the initial space from svd to eps */
237:   if (svd->nini<0) {
238:     EPSSetInitialSpace(cross->eps,-svd->nini,svd->IS);
239:     SlepcBasisDestroy_Private(&svd->nini,&svd->IS);
240:   }
241:   EPSSetUp(cross->eps);
242:   EPSGetDimensions(cross->eps,NULL,&svd->ncv,&svd->mpd);
243:   EPSGetTolerances(cross->eps,NULL,&svd->max_it);
244:   if (svd->tol==PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;

246:   svd->leftbasis = PETSC_FALSE;
247:   SVDAllocateSolution(svd,0);
248:   return 0;
249: }

251: PetscErrorCode SVDSolve_Cross(SVD svd)
252: {
253:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
254:   PetscInt       i;
255:   PetscScalar    lambda;
256:   PetscReal      sigma;

258:   EPSSolve(cross->eps);
259:   EPSGetConverged(cross->eps,&svd->nconv);
260:   EPSGetIterationNumber(cross->eps,&svd->its);
261:   EPSGetConvergedReason(cross->eps,(EPSConvergedReason*)&svd->reason);
262:   for (i=0;i<svd->nconv;i++) {
263:     EPSGetEigenvalue(cross->eps,i,&lambda,NULL);
264:     sigma = PetscRealPart(lambda);
265:     if (svd->ishyperbolic) svd->sigma[i] = PetscSqrtReal(PetscAbsReal(sigma));
266:     else {
268:       if (sigma<0.0) {
269:         PetscInfo(svd,"Negative eigenvalue computed by EPS: %g, resetting to 0\n",(double)sigma);
270:         sigma = 0.0;
271:       }
272:       svd->sigma[i] = PetscSqrtReal(sigma);
273:     }
274:   }
275:   return 0;
276: }

278: PetscErrorCode SVDComputeVectors_Cross(SVD svd)
279: {
280:   SVD_CROSS         *cross = (SVD_CROSS*)svd->data;
281:   PetscInt          i,mloc,ploc;
282:   Vec               u,v,x,uv,w,omega2=NULL;
283:   Mat               Omega;
284:   PetscScalar       *dst,alpha,lambda,*varray;
285:   const PetscScalar *src;
286:   PetscReal         nrm;

288:   if (svd->isgeneralized) {
289:     MatCreateVecs(svd->A,NULL,&u);
290:     VecGetLocalSize(u,&mloc);
291:     MatCreateVecs(svd->B,NULL,&v);
292:     VecGetLocalSize(v,&ploc);
293:     for (i=0;i<svd->nconv;i++) {
294:       BVGetColumn(svd->V,i,&x);
295:       EPSGetEigenpair(cross->eps,i,&lambda,NULL,x,NULL);
296:       MatMult(svd->A,x,u);     /* u_i*c_i/alpha = A*x_i */
297:       VecNormalize(u,NULL);
298:       MatMult(svd->B,x,v);     /* v_i*s_i/alpha = B*x_i */
299:       VecNormalize(v,&nrm);    /* ||v||_2 = s_i/alpha   */
300:       alpha = 1.0/(PetscSqrtReal(1.0+PetscRealPart(lambda))*nrm);    /* alpha=s_i/||v||_2 */
301:       VecScale(x,alpha);
302:       BVRestoreColumn(svd->V,i,&x);
303:       /* copy [u;v] to U[i] */
304:       BVGetColumn(svd->U,i,&uv);
305:       VecGetArrayWrite(uv,&dst);
306:       VecGetArrayRead(u,&src);
307:       PetscArraycpy(dst,src,mloc);
308:       VecRestoreArrayRead(u,&src);
309:       VecGetArrayRead(v,&src);
310:       PetscArraycpy(dst+mloc,src,ploc);
311:       VecRestoreArrayRead(v,&src);
312:       VecRestoreArrayWrite(uv,&dst);
313:       BVRestoreColumn(svd->U,i,&uv);
314:     }
315:     VecDestroy(&v);
316:     VecDestroy(&u);
317:   } else if (svd->ishyperbolic && svd->swapped) {  /* was solved as GHIEP, set u=Omega*u and normalize */
318:     EPSGetOperators(cross->eps,NULL,&Omega);
319:     MatCreateVecs(Omega,&w,NULL);
320:     VecCreateSeq(PETSC_COMM_SELF,svd->ncv,&omega2);
321:     VecGetArrayWrite(omega2,&varray);
322:     for (i=0;i<svd->nconv;i++) {
323:       BVGetColumn(svd->V,i,&v);
324:       EPSGetEigenvector(cross->eps,i,v,NULL);
325:       MatMult(Omega,v,w);
326:       VecDot(v,w,&alpha);
327:       svd->sign[i] = PetscSign(PetscRealPart(alpha));
328:       varray[i] = svd->sign[i];
329:       alpha = 1.0/PetscSqrtScalar(PetscAbsScalar(alpha));
330:       VecScale(w,alpha);
331:       VecCopy(w,v);
332:       BVRestoreColumn(svd->V,i,&v);
333:     }
334:     BVSetSignature(svd->V,omega2);
335:     VecRestoreArrayWrite(omega2,&varray);
336:     VecDestroy(&omega2);
337:     VecDestroy(&w);
338:     SVDComputeVectors_Left(svd);
339:   } else {
340:     for (i=0;i<svd->nconv;i++) {
341:       BVGetColumn(svd->V,i,&v);
342:       EPSGetEigenvector(cross->eps,i,v,NULL);
343:       BVRestoreColumn(svd->V,i,&v);
344:     }
345:     SVDComputeVectors_Left(svd);
346:   }
347:   return 0;
348: }

350: static PetscErrorCode EPSMonitor_Cross(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
351: {
352:   PetscInt       i;
353:   SVD            svd = (SVD)ctx;
354:   PetscScalar    er,ei;

356:   for (i=0;i<PetscMin(nest,svd->ncv);i++) {
357:     er = eigr[i]; ei = eigi[i];
358:     STBackTransform(eps->st,1,&er,&ei);
359:     svd->sigma[i] = PetscSqrtReal(PetscAbsReal(PetscRealPart(er)));
360:     svd->errest[i] = errest[i];
361:   }
362:   SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest);
363:   return 0;
364: }

366: PetscErrorCode SVDSetFromOptions_Cross(SVD svd,PetscOptionItems *PetscOptionsObject)
367: {
368:   PetscBool      set,val;
369:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
370:   ST             st;

372:   PetscOptionsHeadBegin(PetscOptionsObject,"SVD Cross Options");

374:     PetscOptionsBool("-svd_cross_explicitmatrix","Use cross explicit matrix","SVDCrossSetExplicitMatrix",cross->explicitmatrix,&val,&set);
375:     if (set) SVDCrossSetExplicitMatrix(svd,val);

377:   PetscOptionsHeadEnd();

379:   if (!cross->eps) SVDCrossGetEPS(svd,&cross->eps);
380:   if (!cross->explicitmatrix && !cross->usereps) {
381:     /* use as default an ST with shell matrix and Jacobi */
382:     EPSGetST(cross->eps,&st);
383:     STSetMatMode(st,ST_MATMODE_SHELL);
384:   }
385:   EPSSetFromOptions(cross->eps);
386:   return 0;
387: }

389: static PetscErrorCode SVDCrossSetExplicitMatrix_Cross(SVD svd,PetscBool explicitmatrix)
390: {
391:   SVD_CROSS *cross = (SVD_CROSS*)svd->data;

393:   if (cross->explicitmatrix != explicitmatrix) {
394:     cross->explicitmatrix = explicitmatrix;
395:     svd->state = SVD_STATE_INITIAL;
396:   }
397:   return 0;
398: }

400: /*@
401:    SVDCrossSetExplicitMatrix - Indicate if the eigensolver operator A^T*A must
402:    be computed explicitly.

404:    Logically Collective on svd

406:    Input Parameters:
407: +  svd         - singular value solver
408: -  explicitmat - boolean flag indicating if A^T*A is built explicitly

410:    Options Database Key:
411: .  -svd_cross_explicitmatrix <boolean> - Indicates the boolean flag

413:    Level: advanced

415: .seealso: SVDCrossGetExplicitMatrix()
416: @*/
417: PetscErrorCode SVDCrossSetExplicitMatrix(SVD svd,PetscBool explicitmat)
418: {
421:   PetscTryMethod(svd,"SVDCrossSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
422:   return 0;
423: }

425: static PetscErrorCode SVDCrossGetExplicitMatrix_Cross(SVD svd,PetscBool *explicitmat)
426: {
427:   SVD_CROSS *cross = (SVD_CROSS*)svd->data;

429:   *explicitmat = cross->explicitmatrix;
430:   return 0;
431: }

433: /*@
434:    SVDCrossGetExplicitMatrix - Returns the flag indicating if A^T*A is built explicitly.

436:    Not Collective

438:    Input Parameter:
439: .  svd  - singular value solver

441:    Output Parameter:
442: .  explicitmat - the mode flag

444:    Level: advanced

446: .seealso: SVDCrossSetExplicitMatrix()
447: @*/
448: PetscErrorCode SVDCrossGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
449: {
452:   PetscUseMethod(svd,"SVDCrossGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
453:   return 0;
454: }

456: static PetscErrorCode SVDCrossSetEPS_Cross(SVD svd,EPS eps)
457: {
458:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;

460:   PetscObjectReference((PetscObject)eps);
461:   EPSDestroy(&cross->eps);
462:   cross->eps     = eps;
463:   cross->usereps = PETSC_TRUE;
464:   svd->state     = SVD_STATE_INITIAL;
465:   return 0;
466: }

468: /*@
469:    SVDCrossSetEPS - Associate an eigensolver object (EPS) to the
470:    singular value solver.

472:    Collective on svd

474:    Input Parameters:
475: +  svd - singular value solver
476: -  eps - the eigensolver object

478:    Level: advanced

480: .seealso: SVDCrossGetEPS()
481: @*/
482: PetscErrorCode SVDCrossSetEPS(SVD svd,EPS eps)
483: {
487:   PetscTryMethod(svd,"SVDCrossSetEPS_C",(SVD,EPS),(svd,eps));
488:   return 0;
489: }

491: static PetscErrorCode SVDCrossGetEPS_Cross(SVD svd,EPS *eps)
492: {
493:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;

495:   if (!cross->eps) {
496:     EPSCreate(PetscObjectComm((PetscObject)svd),&cross->eps);
497:     PetscObjectIncrementTabLevel((PetscObject)cross->eps,(PetscObject)svd,1);
498:     EPSSetOptionsPrefix(cross->eps,((PetscObject)svd)->prefix);
499:     EPSAppendOptionsPrefix(cross->eps,"svd_cross_");
500:     PetscObjectSetOptions((PetscObject)cross->eps,((PetscObject)svd)->options);
501:     EPSSetWhichEigenpairs(cross->eps,EPS_LARGEST_REAL);
502:     EPSMonitorSet(cross->eps,EPSMonitor_Cross,svd,NULL);
503:   }
504:   *eps = cross->eps;
505:   return 0;
506: }

508: /*@
509:    SVDCrossGetEPS - Retrieve the eigensolver object (EPS) associated
510:    to the singular value solver.

512:    Not Collective

514:    Input Parameter:
515: .  svd - singular value solver

517:    Output Parameter:
518: .  eps - the eigensolver object

520:    Level: advanced

522: .seealso: SVDCrossSetEPS()
523: @*/
524: PetscErrorCode SVDCrossGetEPS(SVD svd,EPS *eps)
525: {
528:   PetscUseMethod(svd,"SVDCrossGetEPS_C",(SVD,EPS*),(svd,eps));
529:   return 0;
530: }

532: PetscErrorCode SVDView_Cross(SVD svd,PetscViewer viewer)
533: {
534:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;
535:   PetscBool      isascii;

537:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
538:   if (isascii) {
539:     if (!cross->eps) SVDCrossGetEPS(svd,&cross->eps);
540:     PetscViewerASCIIPrintf(viewer,"  %s matrix\n",cross->explicitmatrix?"explicit":"implicit");
541:     PetscViewerASCIIPushTab(viewer);
542:     EPSView(cross->eps,viewer);
543:     PetscViewerASCIIPopTab(viewer);
544:   }
545:   return 0;
546: }

548: PetscErrorCode SVDReset_Cross(SVD svd)
549: {
550:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;

552:   EPSReset(cross->eps);
553:   MatDestroy(&cross->C);
554:   MatDestroy(&cross->D);
555:   return 0;
556: }

558: PetscErrorCode SVDDestroy_Cross(SVD svd)
559: {
560:   SVD_CROSS      *cross = (SVD_CROSS*)svd->data;

562:   EPSDestroy(&cross->eps);
563:   PetscFree(svd->data);
564:   PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",NULL);
565:   PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",NULL);
566:   PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",NULL);
567:   PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",NULL);
568:   return 0;
569: }

571: SLEPC_EXTERN PetscErrorCode SVDCreate_Cross(SVD svd)
572: {
573:   SVD_CROSS      *cross;

575:   PetscNew(&cross);
576:   svd->data = (void*)cross;

578:   svd->ops->solve          = SVDSolve_Cross;
579:   svd->ops->solveg         = SVDSolve_Cross;
580:   svd->ops->solveh         = SVDSolve_Cross;
581:   svd->ops->setup          = SVDSetUp_Cross;
582:   svd->ops->setfromoptions = SVDSetFromOptions_Cross;
583:   svd->ops->destroy        = SVDDestroy_Cross;
584:   svd->ops->reset          = SVDReset_Cross;
585:   svd->ops->view           = SVDView_Cross;
586:   svd->ops->computevectors = SVDComputeVectors_Cross;
587:   PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",SVDCrossSetEPS_Cross);
588:   PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",SVDCrossGetEPS_Cross);
589:   PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",SVDCrossSetExplicitMatrix_Cross);
590:   PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",SVDCrossGetExplicitMatrix_Cross);
591:   return 0;
592: }