Actual source code: cross.c

slepc-3.18.1 2022-11-02
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>

 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:   MatShellGetContext(B,&ctx);
 36:   MatMult(ctx->A,x,ctx->w);
 37:   if (ctx->omega && !ctx->swapped) VecPointwiseMult(ctx->w,ctx->w,ctx->omega);
 38:   MatMult(ctx->AT,ctx->w,y);
 39:   return 0;
 40: }

 42: static PetscErrorCode MatGetDiagonal_Cross(Mat B,Vec d)
 43: {
 44:   SVD_CROSS_SHELL   *ctx;
 45:   PetscMPIInt       len;
 46:   PetscInt          N,n,i,j,start,end,ncols;
 47:   PetscScalar       *work1,*work2,*diag;
 48:   const PetscInt    *cols;
 49:   const PetscScalar *vals;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

378:   PetscOptionsHeadEnd();

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

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

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

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

405:    Logically Collective on svd

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

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

414:    Level: advanced

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

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

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

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

437:    Not Collective

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

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

445:    Level: advanced

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

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

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

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

473:    Collective on svd

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

479:    Level: advanced

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

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

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

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

513:    Not Collective

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

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

521:    Level: advanced

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

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

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

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

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

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

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

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

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

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