Actual source code: cyclic.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: "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/epsimpl.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:   MatShellGetContext(B,&ctx);
 28:   MatGetLocalSize(ctx->A,&m,NULL);
 29:   VecGetArrayRead(x,&px);
 30:   VecGetArrayWrite(y,&py);
 31:   VecPlaceArray(ctx->x1,px);
 32:   VecPlaceArray(ctx->x2,px+m);
 33:   VecPlaceArray(ctx->y1,py);
 34:   VecPlaceArray(ctx->y2,py+m);
 35:   MatMult(ctx->A,ctx->x2,ctx->y1);
 36:   MatMult(ctx->AT,ctx->x1,ctx->y2);
 37:   VecResetArray(ctx->x1);
 38:   VecResetArray(ctx->x2);
 39:   VecResetArray(ctx->y1);
 40:   VecResetArray(ctx->y2);
 41:   VecRestoreArrayRead(x,&px);
 42:   VecRestoreArrayWrite(y,&py);
 43:   return 0;
 44: }

 46: static PetscErrorCode MatGetDiagonal_Cyclic(Mat B,Vec diag)
 47: {
 48:   VecSet(diag,0.0);
 49:   return 0;
 50: }

 52: static PetscErrorCode MatDestroy_Cyclic(Mat B)
 53: {
 54:   SVD_CYCLIC_SHELL *ctx;

 56:   MatShellGetContext(B,&ctx);
 57:   VecDestroy(&ctx->x1);
 58:   VecDestroy(&ctx->x2);
 59:   VecDestroy(&ctx->y1);
 60:   VecDestroy(&ctx->y2);
 61:   PetscFree(ctx);
 62:   return 0;
 63: }

 65: /*
 66:    Builds cyclic matrix   C = | 0   A |
 67:                               | AT  0 |
 68: */
 69: static PetscErrorCode SVDCyclicGetCyclicMat(SVD svd,Mat A,Mat AT,Mat *C)
 70: {
 71:   SVD_CYCLIC       *cyclic = (SVD_CYCLIC*)svd->data;
 72:   SVD_CYCLIC_SHELL *ctx;
 73:   PetscInt         i,M,N,m,n,Istart,Iend;
 74:   VecType          vtype;
 75:   Mat              Zm,Zn;
 76: #if defined(PETSC_HAVE_CUDA)
 77:   PetscBool        cuda;
 78: #endif

 80:   MatGetSize(A,&M,&N);
 81:   MatGetLocalSize(A,&m,&n);

 83:   if (cyclic->explicitmatrix) {
 85:     MatCreate(PetscObjectComm((PetscObject)svd),&Zm);
 86:     MatSetSizes(Zm,m,m,M,M);
 87:     MatSetFromOptions(Zm);
 88:     MatSetUp(Zm);
 89:     MatGetOwnershipRange(Zm,&Istart,&Iend);
 90:     for (i=Istart;i<Iend;i++) MatSetValue(Zm,i,i,0.0,INSERT_VALUES);
 91:     MatAssemblyBegin(Zm,MAT_FINAL_ASSEMBLY);
 92:     MatAssemblyEnd(Zm,MAT_FINAL_ASSEMBLY);
 93:     MatCreate(PetscObjectComm((PetscObject)svd),&Zn);
 94:     MatSetSizes(Zn,n,n,N,N);
 95:     MatSetFromOptions(Zn);
 96:     MatSetUp(Zn);
 97:     MatGetOwnershipRange(Zn,&Istart,&Iend);
 98:     for (i=Istart;i<Iend;i++) MatSetValue(Zn,i,i,0.0,INSERT_VALUES);
 99:     MatAssemblyBegin(Zn,MAT_FINAL_ASSEMBLY);
100:     MatAssemblyEnd(Zn,MAT_FINAL_ASSEMBLY);
101:     MatCreateTile(1.0,Zm,1.0,A,1.0,AT,1.0,Zn,C);
102:     MatDestroy(&Zm);
103:     MatDestroy(&Zn);
104:   } else {
105:     PetscNew(&ctx);
106:     ctx->A       = A;
107:     ctx->AT      = AT;
108:     ctx->swapped = svd->swapped;
109:     MatCreateVecsEmpty(A,&ctx->x2,&ctx->x1);
110:     MatCreateVecsEmpty(A,&ctx->y2,&ctx->y1);
111:     MatCreateShell(PetscObjectComm((PetscObject)svd),m+n,m+n,M+N,M+N,ctx,C);
112:     MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Cyclic);
113:     MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_Cyclic);
114: #if defined(PETSC_HAVE_CUDA)
115:     PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&cuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,"");
116:     if (cuda) MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic_CUDA);
117:     else
118: #endif
119:       MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic);
120:     MatGetVecType(A,&vtype);
121:     MatSetVecType(*C,vtype);
122:   }
123:   return 0;
124: }

126: static PetscErrorCode MatMult_ECross(Mat B,Vec x,Vec y)
127: {
128:   SVD_CYCLIC_SHELL  *ctx;
129:   const PetscScalar *px;
130:   PetscScalar       *py;
131:   PetscInt          mn,m,n;

133:   MatShellGetContext(B,&ctx);
134:   MatGetLocalSize(ctx->A,NULL,&n);
135:   VecGetLocalSize(y,&mn);
136:   m = mn-n;
137:   VecGetArrayRead(x,&px);
138:   VecGetArrayWrite(y,&py);
139:   VecPlaceArray(ctx->x1,px);
140:   VecPlaceArray(ctx->x2,px+m);
141:   VecPlaceArray(ctx->y1,py);
142:   VecPlaceArray(ctx->y2,py+m);
143:   VecCopy(ctx->x1,ctx->y1);
144:   MatMult(ctx->A,ctx->x2,ctx->w);
145:   MatMult(ctx->AT,ctx->w,ctx->y2);
146:   VecResetArray(ctx->x1);
147:   VecResetArray(ctx->x2);
148:   VecResetArray(ctx->y1);
149:   VecResetArray(ctx->y2);
150:   VecRestoreArrayRead(x,&px);
151:   VecRestoreArrayWrite(y,&py);
152:   return 0;
153: }

155: static PetscErrorCode MatGetDiagonal_ECross(Mat B,Vec d)
156: {
157:   SVD_CYCLIC_SHELL  *ctx;
158:   PetscScalar       *pd;
159:   PetscMPIInt       len;
160:   PetscInt          mn,m,n,N,i,j,start,end,ncols;
161:   PetscScalar       *work1,*work2,*diag;
162:   const PetscInt    *cols;
163:   const PetscScalar *vals;

165:   MatShellGetContext(B,&ctx);
166:   MatGetLocalSize(ctx->A,NULL,&n);
167:   VecGetLocalSize(d,&mn);
168:   m = mn-n;
169:   VecGetArrayWrite(d,&pd);
170:   VecPlaceArray(ctx->y1,pd);
171:   VecSet(ctx->y1,1.0);
172:   VecResetArray(ctx->y1);
173:   VecPlaceArray(ctx->y2,pd+m);
174:   if (!ctx->diag) {
175:     /* compute diagonal from rows and store in ctx->diag */
176:     VecDuplicate(ctx->y2,&ctx->diag);
177:     MatGetSize(ctx->A,NULL,&N);
178:     PetscCalloc2(N,&work1,N,&work2);
179:     if (ctx->swapped) {
180:       MatGetOwnershipRange(ctx->AT,&start,&end);
181:       for (i=start;i<end;i++) {
182:         MatGetRow(ctx->AT,i,&ncols,NULL,&vals);
183:         for (j=0;j<ncols;j++) work1[i] += vals[j]*vals[j];
184:         MatRestoreRow(ctx->AT,i,&ncols,NULL,&vals);
185:       }
186:     } else {
187:       MatGetOwnershipRange(ctx->A,&start,&end);
188:       for (i=start;i<end;i++) {
189:         MatGetRow(ctx->A,i,&ncols,&cols,&vals);
190:         for (j=0;j<ncols;j++) work1[cols[j]] += vals[j]*vals[j];
191:         MatRestoreRow(ctx->A,i,&ncols,&cols,&vals);
192:       }
193:     }
194:     PetscMPIIntCast(N,&len);
195:     MPIU_Allreduce(work1,work2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)B));
196:     VecGetOwnershipRange(ctx->diag,&start,&end);
197:     VecGetArrayWrite(ctx->diag,&diag);
198:     for (i=start;i<end;i++) diag[i-start] = work2[i];
199:     VecRestoreArrayWrite(ctx->diag,&diag);
200:     PetscFree2(work1,work2);
201:   }
202:   VecCopy(ctx->diag,ctx->y2);
203:   VecResetArray(ctx->y2);
204:   VecRestoreArrayWrite(d,&pd);
205:   return 0;
206: }

208: static PetscErrorCode MatDestroy_ECross(Mat B)
209: {
210:   SVD_CYCLIC_SHELL *ctx;

212:   MatShellGetContext(B,&ctx);
213:   VecDestroy(&ctx->x1);
214:   VecDestroy(&ctx->x2);
215:   VecDestroy(&ctx->y1);
216:   VecDestroy(&ctx->y2);
217:   VecDestroy(&ctx->diag);
218:   VecDestroy(&ctx->w);
219:   PetscFree(ctx);
220:   return 0;
221: }

223: /*
224:    Builds extended cross product matrix   C = | I_m   0  |
225:                                               |  0  AT*A |
226:    t is an auxiliary Vec used to take the dimensions of the upper block
227: */
228: static PetscErrorCode SVDCyclicGetECrossMat(SVD svd,Mat A,Mat AT,Mat *C,Vec t)
229: {
230:   SVD_CYCLIC       *cyclic = (SVD_CYCLIC*)svd->data;
231:   SVD_CYCLIC_SHELL *ctx;
232:   PetscInt         i,M,N,m,n,Istart,Iend;
233:   VecType          vtype;
234:   Mat              Id,Zm,Zn,ATA;
235: #if defined(PETSC_HAVE_CUDA)
236:   PetscBool        cuda;
237: #endif

239:   MatGetSize(A,NULL,&N);
240:   MatGetLocalSize(A,NULL,&n);
241:   VecGetSize(t,&M);
242:   VecGetLocalSize(t,&m);

244:   if (cyclic->explicitmatrix) {
246:     MatCreateConstantDiagonal(PetscObjectComm((PetscObject)svd),m,m,M,M,1.0,&Id);
247:     MatCreate(PetscObjectComm((PetscObject)svd),&Zm);
248:     MatSetSizes(Zm,m,n,M,N);
249:     MatSetFromOptions(Zm);
250:     MatSetUp(Zm);
251:     MatGetOwnershipRange(Zm,&Istart,&Iend);
252:     for (i=Istart;i<Iend;i++) {
253:       if (i<N) MatSetValue(Zm,i,i,0.0,INSERT_VALUES);
254:     }
255:     MatAssemblyBegin(Zm,MAT_FINAL_ASSEMBLY);
256:     MatAssemblyEnd(Zm,MAT_FINAL_ASSEMBLY);
257:     MatCreate(PetscObjectComm((PetscObject)svd),&Zn);
258:     MatSetSizes(Zn,n,m,N,M);
259:     MatSetFromOptions(Zn);
260:     MatSetUp(Zn);
261:     MatGetOwnershipRange(Zn,&Istart,&Iend);
262:     for (i=Istart;i<Iend;i++) {
263:       if (i<m) MatSetValue(Zn,i,i,0.0,INSERT_VALUES);
264:     }
265:     MatAssemblyBegin(Zn,MAT_FINAL_ASSEMBLY);
266:     MatAssemblyEnd(Zn,MAT_FINAL_ASSEMBLY);
267:     MatProductCreate(AT,A,NULL,&ATA);
268:     MatProductSetType(ATA,MATPRODUCT_AB);
269:     MatProductSetFromOptions(ATA);
270:     MatProductSymbolic(ATA);
271:     MatProductNumeric(ATA);
272:     MatCreateTile(1.0,Id,1.0,Zm,1.0,Zn,1.0,ATA,C);
273:     MatDestroy(&Id);
274:     MatDestroy(&Zm);
275:     MatDestroy(&Zn);
276:     MatDestroy(&ATA);
277:   } else {
278:     PetscNew(&ctx);
279:     ctx->A       = A;
280:     ctx->AT      = AT;
281:     ctx->swapped = svd->swapped;
282:     VecDuplicateEmpty(t,&ctx->x1);
283:     VecDuplicateEmpty(t,&ctx->y1);
284:     MatCreateVecsEmpty(A,&ctx->x2,NULL);
285:     MatCreateVecsEmpty(A,&ctx->y2,NULL);
286:     MatCreateVecs(A,NULL,&ctx->w);
287:     MatCreateShell(PetscObjectComm((PetscObject)svd),m+n,m+n,M+N,M+N,ctx,C);
288:     MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ECross);
289:     MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_ECross);
290: #if defined(PETSC_HAVE_CUDA)
291:     PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&cuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,"");
292:     if (cuda) MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross_CUDA);
293:     else
294: #endif
295:       MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross);
296:     MatGetVecType(A,&vtype);
297:     MatSetVecType(*C,vtype);
298:   }
299:   return 0;
300: }

302: /* Convergence test relative to the norm of R (used in GSVD only) */
303: static PetscErrorCode EPSConv_Cyclic(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
304: {
305:   SVD svd = (SVD)ctx;

307:   *errest = res/PetscMax(svd->nrma,svd->nrmb);
308:   return 0;
309: }

311: PetscErrorCode SVDSetUp_Cyclic(SVD svd)
312: {
313:   SVD_CYCLIC        *cyclic = (SVD_CYCLIC*)svd->data;
314:   PetscInt          M,N,m,n,p,k,i,isl,offset,nev,ncv,mpd,maxit;
315:   PetscReal         tol;
316:   const PetscScalar *isa,*oa;
317:   PetscScalar       *va;
318:   EPSProblemType    ptype;
319:   PetscBool         trackall,issinv;
320:   Vec               v,t;
321:   ST                st;
322:   Mat               Omega;
323:   MatType           Atype;

325:   MatGetSize(svd->A,&M,&N);
326:   MatGetLocalSize(svd->A,&m,&n);
327:   if (!cyclic->eps) SVDCyclicGetEPS(svd,&cyclic->eps);
328:   MatDestroy(&cyclic->C);
329:   MatDestroy(&cyclic->D);
330:   if (svd->isgeneralized) {
331:     if (svd->which==SVD_SMALLEST) {  /* alternative pencil */
332:       MatCreateVecs(svd->B,NULL,&t);
333:       SVDCyclicGetCyclicMat(svd,svd->B,svd->BT,&cyclic->C);
334:       SVDCyclicGetECrossMat(svd,svd->A,svd->AT,&cyclic->D,t);
335:     } else {
336:       MatCreateVecs(svd->A,NULL,&t);
337:       SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C);
338:       SVDCyclicGetECrossMat(svd,svd->B,svd->BT,&cyclic->D,t);
339:     }
340:     VecDestroy(&t);
341:     EPSSetOperators(cyclic->eps,cyclic->C,cyclic->D);
342:     EPSGetProblemType(cyclic->eps,&ptype);
343:     if (!ptype) EPSSetProblemType(cyclic->eps,EPS_GHEP);
344:   } else if (svd->ishyperbolic) {
345:     SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C);
346:     MatCreateVecs(cyclic->C,&v,NULL);
347:     VecSet(v,1.0);
348:     VecGetArrayRead(svd->omega,&oa);
349:     VecGetArray(v,&va);
350:     if (svd->swapped) PetscArraycpy(va+m,oa,n);
351:     else PetscArraycpy(va,oa,m);
352:     VecRestoreArrayRead(svd->omega,&oa);
353:     VecRestoreArray(v,&va);
354:     MatGetType(svd->OP,&Atype);
355:     MatCreate(PetscObjectComm((PetscObject)svd),&Omega);
356:     MatSetSizes(Omega,m+n,m+n,M+N,M+N);
357:     MatSetType(Omega,Atype);
358:     MatSetUp(Omega);
359:     MatDiagonalSet(Omega,v,INSERT_VALUES);
360:     EPSSetOperators(cyclic->eps,cyclic->C,Omega);
361:     EPSSetProblemType(cyclic->eps,EPS_GHIEP);
362:     MatDestroy(&Omega);
363:     VecDestroy(&v);
364:   } else {
365:     SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C);
366:     EPSSetOperators(cyclic->eps,cyclic->C,NULL);
367:     EPSSetProblemType(cyclic->eps,EPS_HEP);
368:   }
369:   if (!cyclic->usereps) {
370:     if (svd->which == SVD_LARGEST) {
371:       EPSGetST(cyclic->eps,&st);
372:       PetscObjectTypeCompare((PetscObject)st,STSINVERT,&issinv);
373:       if (issinv) EPSSetWhichEigenpairs(cyclic->eps,EPS_TARGET_MAGNITUDE);
374:       else if (svd->ishyperbolic) EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_MAGNITUDE);
375:       else EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL);
376:     } else {
377:       if (svd->isgeneralized) {  /* computes sigma^{-1} via alternative pencil */
378:         EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL);
379:       } else {
380:         if (svd->ishyperbolic) EPSSetWhichEigenpairs(cyclic->eps,EPS_TARGET_MAGNITUDE);
381:         else EPSSetEigenvalueComparison(cyclic->eps,SlepcCompareSmallestPosReal,NULL);
382:         EPSSetTarget(cyclic->eps,0.0);
383:       }
384:     }
385:     EPSGetDimensions(cyclic->eps,&nev,&ncv,&mpd);
387:     nev = PetscMax(nev,2*svd->nsv);
388:     if (ncv==PETSC_DEFAULT && svd->ncv!=PETSC_DEFAULT) ncv = PetscMax(3*svd->nsv,svd->ncv);
389:     if (mpd==PETSC_DEFAULT && svd->mpd!=PETSC_DEFAULT) mpd = svd->mpd;
390:     EPSSetDimensions(cyclic->eps,nev,ncv,mpd);
391:     EPSGetTolerances(cyclic->eps,&tol,&maxit);
392:     if (tol==PETSC_DEFAULT) tol = svd->tol==PETSC_DEFAULT? SLEPC_DEFAULT_TOL/10.0: svd->tol;
393:     if (maxit==PETSC_DEFAULT && svd->max_it!=PETSC_DEFAULT) maxit = svd->max_it;
394:     EPSSetTolerances(cyclic->eps,tol,maxit);
395:     switch (svd->conv) {
396:     case SVD_CONV_ABS:
397:       EPSSetConvergenceTest(cyclic->eps,EPS_CONV_ABS);break;
398:     case SVD_CONV_REL:
399:       EPSSetConvergenceTest(cyclic->eps,EPS_CONV_REL);break;
400:     case SVD_CONV_NORM:
401:       if (svd->isgeneralized) {
402:         if (!svd->nrma) MatNorm(svd->OP,NORM_INFINITY,&svd->nrma);
403:         if (!svd->nrmb) MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb);
404:         EPSSetConvergenceTestFunction(cyclic->eps,EPSConv_Cyclic,svd,NULL);
405:       } else {
406:         EPSSetConvergenceTest(cyclic->eps,EPS_CONV_NORM);break;
407:       }
408:       break;
409:     case SVD_CONV_MAXIT:
410:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Maxit convergence test not supported in this solver");
411:     case SVD_CONV_USER:
412:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
413:     }
414:   }
415:   SVDCheckUnsupported(svd,SVD_FEATURE_STOPPING);
416:   /* Transfer the trackall option from svd to eps */
417:   SVDGetTrackAll(svd,&trackall);
418:   EPSSetTrackAll(cyclic->eps,trackall);
419:   /* Transfer the initial subspace from svd to eps */
420:   if (svd->nini<0 || svd->ninil<0) {
421:     for (i=0;i<-PetscMin(svd->nini,svd->ninil);i++) {
422:       MatCreateVecs(cyclic->C,&v,NULL);
423:       VecGetArrayWrite(v,&va);
424:       if (svd->isgeneralized) MatGetLocalSize(svd->B,&p,NULL);
425:       k = (svd->isgeneralized && svd->which==SVD_SMALLEST)? p: m;  /* size of upper block row */
426:       if (i<-svd->ninil) {
427:         VecGetArrayRead(svd->ISL[i],&isa);
428:         if (svd->isgeneralized) {
429:           VecGetLocalSize(svd->ISL[i],&isl);
431:           offset = (svd->which==SVD_SMALLEST)? m: 0;
432:           PetscArraycpy(va,isa+offset,k);
433:         } else {
434:           VecGetLocalSize(svd->ISL[i],&isl);
436:           PetscArraycpy(va,isa,k);
437:         }
438:         VecRestoreArrayRead(svd->IS[i],&isa);
439:       } else PetscArrayzero(&va,k);
440:       if (i<-svd->nini) {
441:         VecGetLocalSize(svd->IS[i],&isl);
443:         VecGetArrayRead(svd->IS[i],&isa);
444:         PetscArraycpy(va+k,isa,n);
445:         VecRestoreArrayRead(svd->IS[i],&isa);
446:       } else PetscArrayzero(va+k,n);
447:       VecRestoreArrayWrite(v,&va);
448:       VecDestroy(&svd->IS[i]);
449:       svd->IS[i] = v;
450:     }
451:     svd->nini = PetscMin(svd->nini,svd->ninil);
452:     EPSSetInitialSpace(cyclic->eps,-svd->nini,svd->IS);
453:     SlepcBasisDestroy_Private(&svd->nini,&svd->IS);
454:     SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL);
455:   }
456:   EPSSetUp(cyclic->eps);
457:   EPSGetDimensions(cyclic->eps,NULL,&svd->ncv,&svd->mpd);
458:   svd->ncv = PetscMin(svd->ncv,PetscMin(M,N));
459:   EPSGetTolerances(cyclic->eps,NULL,&svd->max_it);
460:   if (svd->tol==PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;

462:   svd->leftbasis = PETSC_TRUE;
463:   SVDAllocateSolution(svd,0);
464:   return 0;
465: }

467: PetscErrorCode SVDCyclicCheckEigenvalue(SVD svd,PetscScalar er,PetscScalar ei,PetscReal *sigma,PetscBool *isreal)
468: {
469:   if (svd->ishyperbolic && PetscDefined(USE_COMPLEX) && PetscAbsReal(PetscImaginaryPart(er))>10*PetscAbsReal(PetscRealPart(er))) {
470:     *sigma = PetscImaginaryPart(er);
471:     if (isreal) *isreal = PETSC_FALSE;
472:   } else if (svd->ishyperbolic && !PetscDefined(USE_COMPLEX) && PetscAbsScalar(ei)>10*PetscAbsScalar(er)) {
473:     *sigma = PetscRealPart(ei);
474:     if (isreal) *isreal = PETSC_FALSE;
475:   } else {
476:     *sigma = PetscRealPart(er);
477:     if (isreal) *isreal = PETSC_TRUE;
478:   }
479:   return 0;
480: }

482: PetscErrorCode SVDSolve_Cyclic(SVD svd)
483: {
484:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;
485:   PetscInt       i,j,nconv;
486:   PetscScalar    er,ei;
487:   PetscReal      sigma;

489:   EPSSolve(cyclic->eps);
490:   EPSGetConverged(cyclic->eps,&nconv);
491:   EPSGetIterationNumber(cyclic->eps,&svd->its);
492:   EPSGetConvergedReason(cyclic->eps,(EPSConvergedReason*)&svd->reason);
493:   for (i=0,j=0;i<nconv;i++) {
494:     EPSGetEigenvalue(cyclic->eps,i,&er,&ei);
495:     SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL);
496:     if (sigma>0.0) {
497:       if (svd->isgeneralized && svd->which==SVD_SMALLEST) svd->sigma[j] = 1.0/sigma;
498:       else svd->sigma[j] = sigma;
499:       j++;
500:     }
501:   }
502:   svd->nconv = j;
503:   return 0;
504: }

506: static PetscErrorCode SVDComputeVectors_Cyclic_Standard(SVD svd)
507: {
508:   SVD_CYCLIC        *cyclic = (SVD_CYCLIC*)svd->data;
509:   PetscInt          i,j,m,nconv;
510:   PetscScalar       er,ei;
511:   PetscReal         sigma;
512:   const PetscScalar *px;
513:   Vec               x,x1,x2;

515:   MatCreateVecs(cyclic->C,&x,NULL);
516:   MatGetLocalSize(svd->A,&m,NULL);
517:   MatCreateVecsEmpty(svd->A,&x2,&x1);
518:   EPSGetConverged(cyclic->eps,&nconv);
519:   for (i=0,j=0;i<nconv;i++) {
520:     EPSGetEigenpair(cyclic->eps,i,&er,&ei,x,NULL);
521:     SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL);
522:     if (sigma<0.0) continue;
523:     VecGetArrayRead(x,&px);
524:     VecPlaceArray(x1,px);
525:     VecPlaceArray(x2,px+m);
526:     BVInsertVec(svd->U,j,x1);
527:     BVScaleColumn(svd->U,j,PETSC_SQRT2);
528:     BVInsertVec(svd->V,j,x2);
529:     BVScaleColumn(svd->V,j,PETSC_SQRT2);
530:     VecResetArray(x1);
531:     VecResetArray(x2);
532:     VecRestoreArrayRead(x,&px);
533:     j++;
534:   }
535:   VecDestroy(&x);
536:   VecDestroy(&x1);
537:   VecDestroy(&x2);
538:   return 0;
539: }

541: static PetscErrorCode SVDComputeVectors_Cyclic_Generalized(SVD svd)
542: {
543:   SVD_CYCLIC        *cyclic = (SVD_CYCLIC*)svd->data;
544:   PetscInt          i,j,m,p,nconv;
545:   PetscScalar       *dst,er,ei;
546:   PetscReal         sigma;
547:   const PetscScalar *src,*px;
548:   Vec               u,v,x,x1,x2,uv;

550:   MatGetLocalSize(svd->A,&m,NULL);
551:   MatGetLocalSize(svd->B,&p,NULL);
552:   MatCreateVecs(cyclic->C,&x,NULL);
553:   if (svd->which==SVD_SMALLEST) MatCreateVecsEmpty(svd->B,&x1,&x2);
554:   else MatCreateVecsEmpty(svd->A,&x2,&x1);
555:   MatCreateVecs(svd->A,NULL,&u);
556:   MatCreateVecs(svd->B,NULL,&v);
557:   EPSGetConverged(cyclic->eps,&nconv);
558:   for (i=0,j=0;i<nconv;i++) {
559:     EPSGetEigenpair(cyclic->eps,i,&er,&ei,x,NULL);
560:     SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL);
561:     if (sigma<0.0) continue;
562:     if (svd->which==SVD_SMALLEST) {
563:       /* evec_i = 1/sqrt(2)*[ v_i; w_i ],  w_i = x_i/c_i */
564:       VecGetArrayRead(x,&px);
565:       VecPlaceArray(x2,px);
566:       VecPlaceArray(x1,px+p);
567:       VecCopy(x2,v);
568:       VecScale(v,PETSC_SQRT2);  /* v_i = sqrt(2)*evec_i_1 */
569:       VecScale(x1,PETSC_SQRT2); /* w_i = sqrt(2)*evec_i_2 */
570:       MatMult(svd->A,x1,u);     /* A*w_i = u_i */
571:       VecScale(x1,1.0/PetscSqrtScalar(1.0+sigma*sigma));  /* x_i = w_i*c_i */
572:       BVInsertVec(svd->V,j,x1);
573:       VecResetArray(x2);
574:       VecResetArray(x1);
575:       VecRestoreArrayRead(x,&px);
576:     } else {
577:       /* evec_i = 1/sqrt(2)*[ u_i; w_i ],  w_i = x_i/s_i */
578:       VecGetArrayRead(x,&px);
579:       VecPlaceArray(x1,px);
580:       VecPlaceArray(x2,px+m);
581:       VecCopy(x1,u);
582:       VecScale(u,PETSC_SQRT2);  /* u_i = sqrt(2)*evec_i_1 */
583:       VecScale(x2,PETSC_SQRT2); /* w_i = sqrt(2)*evec_i_2 */
584:       MatMult(svd->B,x2,v);     /* B*w_i = v_i */
585:       VecScale(x2,1.0/PetscSqrtScalar(1.0+sigma*sigma));  /* x_i = w_i*s_i */
586:       BVInsertVec(svd->V,j,x2);
587:       VecResetArray(x1);
588:       VecResetArray(x2);
589:       VecRestoreArrayRead(x,&px);
590:     }
591:     /* copy [u;v] to U[j] */
592:     BVGetColumn(svd->U,j,&uv);
593:     VecGetArrayWrite(uv,&dst);
594:     VecGetArrayRead(u,&src);
595:     PetscArraycpy(dst,src,m);
596:     VecRestoreArrayRead(u,&src);
597:     VecGetArrayRead(v,&src);
598:     PetscArraycpy(dst+m,src,p);
599:     VecRestoreArrayRead(v,&src);
600:     VecRestoreArrayWrite(uv,&dst);
601:     BVRestoreColumn(svd->U,j,&uv);
602:     j++;
603:   }
604:   VecDestroy(&x);
605:   VecDestroy(&x1);
606:   VecDestroy(&x2);
607:   VecDestroy(&u);
608:   VecDestroy(&v);
609:   return 0;
610: }

612: #if defined(PETSC_USE_COMPLEX)
613: /* VecMaxAbs: returns the entry of x that has max(abs(x(i))), using w as a workspace vector */
614: static PetscErrorCode VecMaxAbs(Vec x,Vec w,PetscScalar *v)
615: {
616:   PetscMPIInt       size,rank,root;
617:   const PetscScalar *xx;
618:   const PetscInt    *ranges;
619:   PetscReal         val;
620:   PetscInt          p;

622:   VecCopy(x,w);
623:   VecAbs(w);
624:   VecMax(w,&p,&val);
625:   MPI_Comm_size(PetscObjectComm((PetscObject)x),&size);
626:   MPI_Comm_rank(PetscObjectComm((PetscObject)x),&rank);
627:   VecGetOwnershipRanges(x,&ranges);
628:   for (root=0;root<size;root++) if (p>=ranges[root] && p<ranges[root+1]) break;
629:   if (rank==root) {
630:     VecGetArrayRead(x,&xx);
631:     *v = xx[p-ranges[root]];
632:     VecRestoreArrayRead(x,&xx);
633:   }
634:   MPI_Bcast(v,1,MPIU_SCALAR,root,PetscObjectComm((PetscObject)x));
635:   return 0;
636: }
637: #endif

639: static PetscErrorCode SVDComputeVectors_Cyclic_Hyperbolic(SVD svd)
640: {
641:   SVD_CYCLIC        *cyclic = (SVD_CYCLIC*)svd->data;
642:   PetscInt          i,j,m,n,N,nconv;
643:   PetscScalar       er,ei;
644:   PetscReal         sigma,nrm;
645:   PetscBool         isreal;
646:   const PetscScalar *px;
647:   Vec               u,x,xi=NULL,x1,x2,x1i=NULL,x2i;
648:   Mat               Omega;
649:   MatType           Atype;
650:   BV                U=NULL,V=NULL;
651: #if !defined(PETSC_USE_COMPLEX)
652:   const PetscScalar *pxi;
653:   PetscReal         nrmr,nrmi;
654: #else
655:   PetscScalar       alpha;
656: #endif

658:   MatCreateVecs(cyclic->C,&x,svd->ishyperbolic?&xi:NULL);
659:   MatGetLocalSize(svd->A,&m,NULL);
660:   MatCreateVecsEmpty(svd->OP,&x2,&x1);
661: #if defined(PETSC_USE_COMPLEX)
662:   MatCreateVecs(svd->OP,&x2i,&x1i);
663: #else
664:   MatCreateVecsEmpty(svd->OP,&x2i,&x1i);
665: #endif
666:   /* set-up Omega-normalization of U */
667:   U = svd->swapped? svd->V: svd->U;
668:   V = svd->swapped? svd->U: svd->V;
669:   MatGetType(svd->A,&Atype);
670:   BVGetSizes(U,&n,&N,NULL);
671:   MatCreate(PetscObjectComm((PetscObject)svd),&Omega);
672:   MatSetSizes(Omega,n,n,N,N);
673:   MatSetType(Omega,Atype);
674:   MatSetUp(Omega);
675:   MatDiagonalSet(Omega,svd->omega,INSERT_VALUES);
676:   BVSetMatrix(U,Omega,PETSC_TRUE);
677:   MatDestroy(&Omega);

679:   EPSGetConverged(cyclic->eps,&nconv);
680:   for (i=0,j=0;i<nconv;i++) {
681:     EPSGetEigenpair(cyclic->eps,i,&er,&ei,x,xi);
682:     SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,&isreal);
683:     if (sigma<0.0) continue;
684:     VecGetArrayRead(x,&px);
685:     if (svd->swapped) {
686:       VecPlaceArray(x2,px);
687:       VecPlaceArray(x1,px+m);
688:     } else {
689:       VecPlaceArray(x1,px);
690:       VecPlaceArray(x2,px+n);
691:     }
692: #if defined(PETSC_USE_COMPLEX)
693:     BVInsertVec(U,j,x1);
694:     BVInsertVec(V,j,x2);
695:     if (!isreal) {
696:       VecMaxAbs(x1,x1i,&alpha);
697:       BVScaleColumn(U,j,PetscAbsScalar(alpha)/alpha);
698:       BVScaleColumn(V,j,PetscAbsScalar(alpha)/(alpha*PETSC_i));
699:     }
700: #else
701:     VecGetArrayRead(xi,&pxi);
702:     if (svd->swapped) {
703:       VecPlaceArray(x2i,pxi);
704:       VecPlaceArray(x1i,pxi+m);
705:     } else {
706:       VecPlaceArray(x1i,pxi);
707:       VecPlaceArray(x2i,pxi+n);
708:     }
709:     VecNorm(x2,NORM_2,&nrmr);
710:     VecNorm(x2i,NORM_2,&nrmi);
711:     if (nrmi>nrmr) {
712:       if (isreal) {
713:         BVInsertVec(U,j,x1i);
714:         BVInsertVec(V,j,x2i);
715:       } else {
716:         BVInsertVec(U,j,x1);
717:         BVInsertVec(V,j,x2i);
718:       }
719:     } else {
720:       if (isreal) {
721:         BVInsertVec(U,j,x1);
722:         BVInsertVec(V,j,x2);
723:       } else {
724:         BVInsertVec(U,j,x1i);
725:         BVScaleColumn(U,j,-1.0);
726:         BVInsertVec(V,j,x2);
727:       }
728:     }
729:     VecResetArray(x1i);
730:     VecResetArray(x2i);
731:     VecRestoreArrayRead(xi,&pxi);
732: #endif
733:     VecResetArray(x1);
734:     VecResetArray(x2);
735:     VecRestoreArrayRead(x,&px);
736:     BVGetColumn(U,j,&u);
737:     VecPointwiseMult(u,u,svd->omega);
738:     BVRestoreColumn(U,j,&u);
739:     BVNormColumn(U,j,NORM_2,&nrm);
740:     BVScaleColumn(U,j,1.0/PetscAbs(nrm));
741:     svd->sign[j] = PetscSign(nrm);
742:     BVNormColumn(V,j,NORM_2,&nrm);
743:     BVScaleColumn(V,j,1.0/nrm);
744:     j++;
745:   }
746:   VecDestroy(&x);
747:   VecDestroy(&x1);
748:   VecDestroy(&x2);
749:   VecDestroy(&xi);
750:   VecDestroy(&x1i);
751:   VecDestroy(&x2i);
752:   return 0;
753: }

755: PetscErrorCode SVDComputeVectors_Cyclic(SVD svd)
756: {
757:   switch (svd->problem_type) {
758:     case SVD_STANDARD:
759:       SVDComputeVectors_Cyclic_Standard(svd);
760:       break;
761:     case SVD_GENERALIZED:
762:       SVDComputeVectors_Cyclic_Generalized(svd);
763:       break;
764:     case SVD_HYPERBOLIC:
765:       SVDComputeVectors_Cyclic_Hyperbolic(svd);
766:       break;
767:     default:
768:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONG,"Unknown singular value problem type");
769:   }
770:   return 0;
771: }

773: static PetscErrorCode EPSMonitor_Cyclic(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
774: {
775:   PetscInt       i,j;
776:   SVD            svd = (SVD)ctx;
777:   PetscScalar    er,ei;
778:   PetscReal      sigma;

780:   nconv = 0;
781:   for (i=0,j=0;i<PetscMin(nest,svd->ncv);i++) {
782:     er = eigr[i]; ei = eigi[i];
783:     STBackTransform(eps->st,1,&er,&ei);
784:     SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL);
785:     if (sigma>0.0) {
786:       svd->sigma[j]  = sigma;
787:       svd->errest[j] = errest[i];
788:       if (errest[i] && errest[i] < svd->tol) nconv++;
789:       j++;
790:     }
791:   }
792:   nest = j;
793:   SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest);
794:   return 0;
795: }

797: PetscErrorCode SVDSetFromOptions_Cyclic(SVD svd,PetscOptionItems *PetscOptionsObject)
798: {
799:   PetscBool      set,val;
800:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;
801:   ST             st;

803:   PetscOptionsHeadBegin(PetscOptionsObject,"SVD Cyclic Options");

805:     PetscOptionsBool("-svd_cyclic_explicitmatrix","Use cyclic explicit matrix","SVDCyclicSetExplicitMatrix",cyclic->explicitmatrix,&val,&set);
806:     if (set) SVDCyclicSetExplicitMatrix(svd,val);

808:   PetscOptionsHeadEnd();

810:   if (!cyclic->eps) SVDCyclicGetEPS(svd,&cyclic->eps);
811:   if (!cyclic->explicitmatrix && !cyclic->usereps) {
812:     /* use as default an ST with shell matrix and Jacobi */
813:     EPSGetST(cyclic->eps,&st);
814:     STSetMatMode(st,ST_MATMODE_SHELL);
815:   }
816:   EPSSetFromOptions(cyclic->eps);
817:   return 0;
818: }

820: static PetscErrorCode SVDCyclicSetExplicitMatrix_Cyclic(SVD svd,PetscBool explicitmat)
821: {
822:   SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;

824:   if (cyclic->explicitmatrix != explicitmat) {
825:     cyclic->explicitmatrix = explicitmat;
826:     svd->state = SVD_STATE_INITIAL;
827:   }
828:   return 0;
829: }

831: /*@
832:    SVDCyclicSetExplicitMatrix - Indicate if the eigensolver operator
833:    H(A) = [ 0  A ; A^T 0 ] must be computed explicitly.

835:    Logically Collective on svd

837:    Input Parameters:
838: +  svd         - singular value solver
839: -  explicitmat - boolean flag indicating if H(A) is built explicitly

841:    Options Database Key:
842: .  -svd_cyclic_explicitmatrix <boolean> - Indicates the boolean flag

844:    Level: advanced

846: .seealso: SVDCyclicGetExplicitMatrix()
847: @*/
848: PetscErrorCode SVDCyclicSetExplicitMatrix(SVD svd,PetscBool explicitmat)
849: {
852:   PetscTryMethod(svd,"SVDCyclicSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
853:   return 0;
854: }

856: static PetscErrorCode SVDCyclicGetExplicitMatrix_Cyclic(SVD svd,PetscBool *explicitmat)
857: {
858:   SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;

860:   *explicitmat = cyclic->explicitmatrix;
861:   return 0;
862: }

864: /*@
865:    SVDCyclicGetExplicitMatrix - Returns the flag indicating if H(A) is built explicitly.

867:    Not Collective

869:    Input Parameter:
870: .  svd  - singular value solver

872:    Output Parameter:
873: .  explicitmat - the mode flag

875:    Level: advanced

877: .seealso: SVDCyclicSetExplicitMatrix()
878: @*/
879: PetscErrorCode SVDCyclicGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
880: {
883:   PetscUseMethod(svd,"SVDCyclicGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
884:   return 0;
885: }

887: static PetscErrorCode SVDCyclicSetEPS_Cyclic(SVD svd,EPS eps)
888: {
889:   SVD_CYCLIC      *cyclic = (SVD_CYCLIC*)svd->data;

891:   PetscObjectReference((PetscObject)eps);
892:   EPSDestroy(&cyclic->eps);
893:   cyclic->eps     = eps;
894:   cyclic->usereps = PETSC_TRUE;
895:   svd->state      = SVD_STATE_INITIAL;
896:   return 0;
897: }

899: /*@
900:    SVDCyclicSetEPS - Associate an eigensolver object (EPS) to the
901:    singular value solver.

903:    Collective on svd

905:    Input Parameters:
906: +  svd - singular value solver
907: -  eps - the eigensolver object

909:    Level: advanced

911: .seealso: SVDCyclicGetEPS()
912: @*/
913: PetscErrorCode SVDCyclicSetEPS(SVD svd,EPS eps)
914: {
918:   PetscTryMethod(svd,"SVDCyclicSetEPS_C",(SVD,EPS),(svd,eps));
919:   return 0;
920: }

922: static PetscErrorCode SVDCyclicGetEPS_Cyclic(SVD svd,EPS *eps)
923: {
924:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;

926:   if (!cyclic->eps) {
927:     EPSCreate(PetscObjectComm((PetscObject)svd),&cyclic->eps);
928:     PetscObjectIncrementTabLevel((PetscObject)cyclic->eps,(PetscObject)svd,1);
929:     EPSSetOptionsPrefix(cyclic->eps,((PetscObject)svd)->prefix);
930:     EPSAppendOptionsPrefix(cyclic->eps,"svd_cyclic_");
931:     PetscObjectSetOptions((PetscObject)cyclic->eps,((PetscObject)svd)->options);
932:     EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL);
933:     EPSMonitorSet(cyclic->eps,EPSMonitor_Cyclic,svd,NULL);
934:   }
935:   *eps = cyclic->eps;
936:   return 0;
937: }

939: /*@
940:    SVDCyclicGetEPS - Retrieve the eigensolver object (EPS) associated
941:    to the singular value solver.

943:    Not Collective

945:    Input Parameter:
946: .  svd - singular value solver

948:    Output Parameter:
949: .  eps - the eigensolver object

951:    Level: advanced

953: .seealso: SVDCyclicSetEPS()
954: @*/
955: PetscErrorCode SVDCyclicGetEPS(SVD svd,EPS *eps)
956: {
959:   PetscUseMethod(svd,"SVDCyclicGetEPS_C",(SVD,EPS*),(svd,eps));
960:   return 0;
961: }

963: PetscErrorCode SVDView_Cyclic(SVD svd,PetscViewer viewer)
964: {
965:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;
966:   PetscBool      isascii;

968:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
969:   if (isascii) {
970:     if (!cyclic->eps) SVDCyclicGetEPS(svd,&cyclic->eps);
971:     PetscViewerASCIIPrintf(viewer,"  %s matrix\n",cyclic->explicitmatrix?"explicit":"implicit");
972:     PetscViewerASCIIPushTab(viewer);
973:     EPSView(cyclic->eps,viewer);
974:     PetscViewerASCIIPopTab(viewer);
975:   }
976:   return 0;
977: }

979: PetscErrorCode SVDReset_Cyclic(SVD svd)
980: {
981:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;

983:   EPSReset(cyclic->eps);
984:   MatDestroy(&cyclic->C);
985:   MatDestroy(&cyclic->D);
986:   return 0;
987: }

989: PetscErrorCode SVDDestroy_Cyclic(SVD svd)
990: {
991:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;

993:   EPSDestroy(&cyclic->eps);
994:   PetscFree(svd->data);
995:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",NULL);
996:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",NULL);
997:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",NULL);
998:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",NULL);
999:   return 0;
1000: }

1002: SLEPC_EXTERN PetscErrorCode SVDCreate_Cyclic(SVD svd)
1003: {
1004:   SVD_CYCLIC     *cyclic;

1006:   PetscNew(&cyclic);
1007:   svd->data                = (void*)cyclic;
1008:   svd->ops->solve          = SVDSolve_Cyclic;
1009:   svd->ops->solveg         = SVDSolve_Cyclic;
1010:   svd->ops->solveh         = SVDSolve_Cyclic;
1011:   svd->ops->setup          = SVDSetUp_Cyclic;
1012:   svd->ops->setfromoptions = SVDSetFromOptions_Cyclic;
1013:   svd->ops->destroy        = SVDDestroy_Cyclic;
1014:   svd->ops->reset          = SVDReset_Cyclic;
1015:   svd->ops->view           = SVDView_Cyclic;
1016:   svd->ops->computevectors = SVDComputeVectors_Cyclic;
1017:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",SVDCyclicSetEPS_Cyclic);
1018:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",SVDCyclicGetEPS_Cyclic);
1019:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",SVDCyclicSetExplicitMatrix_Cyclic);
1020:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",SVDCyclicGetExplicitMatrix_Cyclic);
1021:   return 0;
1022: }