Actual source code: cyclic.c

slepc-3.20.2 2024-03-15
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/bvimpl.h>
 18: #include "cyclic.h"

 20: static PetscErrorCode MatMult_Cyclic(Mat B,Vec x,Vec y)
 21: {
 22:   SVD_CYCLIC_SHELL  *ctx;
 23:   const PetscScalar *px;
 24:   PetscScalar       *py;
 25:   PetscInt          m;

 27:   PetscFunctionBegin;
 28:   PetscCall(MatShellGetContext(B,&ctx));
 29:   PetscCall(MatGetLocalSize(ctx->A,&m,NULL));
 30:   PetscCall(VecGetArrayRead(x,&px));
 31:   PetscCall(VecGetArrayWrite(y,&py));
 32:   PetscCall(VecPlaceArray(ctx->x1,px));
 33:   PetscCall(VecPlaceArray(ctx->x2,px+m));
 34:   PetscCall(VecPlaceArray(ctx->y1,py));
 35:   PetscCall(VecPlaceArray(ctx->y2,py+m));
 36:   PetscCall(MatMult(ctx->A,ctx->x2,ctx->y1));
 37:   PetscCall(MatMult(ctx->AT,ctx->x1,ctx->y2));
 38:   PetscCall(VecResetArray(ctx->x1));
 39:   PetscCall(VecResetArray(ctx->x2));
 40:   PetscCall(VecResetArray(ctx->y1));
 41:   PetscCall(VecResetArray(ctx->y2));
 42:   PetscCall(VecRestoreArrayRead(x,&px));
 43:   PetscCall(VecRestoreArrayWrite(y,&py));
 44:   PetscFunctionReturn(PETSC_SUCCESS);
 45: }

 47: static PetscErrorCode MatGetDiagonal_Cyclic(Mat B,Vec diag)
 48: {
 49:   PetscFunctionBegin;
 50:   PetscCall(VecSet(diag,0.0));
 51:   PetscFunctionReturn(PETSC_SUCCESS);
 52: }

 54: static PetscErrorCode MatDestroy_Cyclic(Mat B)
 55: {
 56:   SVD_CYCLIC_SHELL *ctx;

 58:   PetscFunctionBegin;
 59:   PetscCall(MatShellGetContext(B,&ctx));
 60:   PetscCall(VecDestroy(&ctx->x1));
 61:   PetscCall(VecDestroy(&ctx->x2));
 62:   PetscCall(VecDestroy(&ctx->y1));
 63:   PetscCall(VecDestroy(&ctx->y2));
 64:   if (ctx->misaligned) {
 65:     PetscCall(VecDestroy(&ctx->wx2));
 66:     PetscCall(VecDestroy(&ctx->wy2));
 67:   }
 68:   PetscCall(PetscFree(ctx));
 69:   PetscFunctionReturn(PETSC_SUCCESS);
 70: }

 72: /*
 73:    Builds cyclic matrix   C = | 0   A |
 74:                               | AT  0 |
 75: */
 76: static PetscErrorCode SVDCyclicGetCyclicMat(SVD svd,Mat A,Mat AT,Mat *C)
 77: {
 78:   SVD_CYCLIC       *cyclic = (SVD_CYCLIC*)svd->data;
 79:   SVD_CYCLIC_SHELL *ctx;
 80:   PetscInt         i,M,N,m,n,Istart,Iend;
 81:   VecType          vtype;
 82:   Mat              Zm,Zn;
 83: #if defined(PETSC_HAVE_CUDA)
 84:   PetscBool        cuda;
 85:   const PetscInt   *ranges;
 86:   PetscMPIInt      size;
 87: #endif

 89:   PetscFunctionBegin;
 90:   PetscCall(MatGetSize(A,&M,&N));
 91:   PetscCall(MatGetLocalSize(A,&m,&n));

 93:   if (cyclic->explicitmatrix) {
 94:     PetscCheck(svd->expltrans,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Cannot use explicit cyclic matrix with implicit transpose");
 95:     PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Zm));
 96:     PetscCall(MatSetSizes(Zm,m,m,M,M));
 97:     PetscCall(MatSetFromOptions(Zm));
 98:     PetscCall(MatSetUp(Zm));
 99:     PetscCall(MatGetOwnershipRange(Zm,&Istart,&Iend));
100:     for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(Zm,i,i,0.0,INSERT_VALUES));
101:     PetscCall(MatAssemblyBegin(Zm,MAT_FINAL_ASSEMBLY));
102:     PetscCall(MatAssemblyEnd(Zm,MAT_FINAL_ASSEMBLY));
103:     PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Zn));
104:     PetscCall(MatSetSizes(Zn,n,n,N,N));
105:     PetscCall(MatSetFromOptions(Zn));
106:     PetscCall(MatSetUp(Zn));
107:     PetscCall(MatGetOwnershipRange(Zn,&Istart,&Iend));
108:     for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(Zn,i,i,0.0,INSERT_VALUES));
109:     PetscCall(MatAssemblyBegin(Zn,MAT_FINAL_ASSEMBLY));
110:     PetscCall(MatAssemblyEnd(Zn,MAT_FINAL_ASSEMBLY));
111:     PetscCall(MatCreateTile(1.0,Zm,1.0,A,1.0,AT,1.0,Zn,C));
112:     PetscCall(MatDestroy(&Zm));
113:     PetscCall(MatDestroy(&Zn));
114:   } else {
115:     PetscCall(PetscNew(&ctx));
116:     ctx->A       = A;
117:     ctx->AT      = AT;
118:     ctx->swapped = svd->swapped;
119:     PetscCall(MatCreateVecsEmpty(A,&ctx->x2,&ctx->x1));
120:     PetscCall(MatCreateVecsEmpty(A,&ctx->y2,&ctx->y1));
121:     PetscCall(MatCreateShell(PetscObjectComm((PetscObject)svd),m+n,m+n,M+N,M+N,ctx,C));
122:     PetscCall(MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Cyclic));
123:     PetscCall(MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_Cyclic));
124: #if defined(PETSC_HAVE_CUDA)
125:     PetscCall(PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&cuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,""));
126:     if (cuda) PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic_CUDA));
127:     else
128: #endif
129:       PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic));
130:     PetscCall(MatGetVecType(A,&vtype));
131:     PetscCall(MatSetVecType(*C,vtype));
132: #if defined(PETSC_HAVE_CUDA)
133:     if (cuda) {
134:       /* check alignment of bottom block */
135:       PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ctx->x1),&size));
136:       PetscCall(VecGetOwnershipRanges(ctx->x1,&ranges));
137:       for (i=0;i<size;i++) {
138:         ctx->misaligned = (((ranges[i+1]-ranges[i])*sizeof(PetscScalar))%16)? PETSC_TRUE: PETSC_FALSE;
139:         if (ctx->misaligned) break;
140:       }
141:       if (ctx->misaligned) {  /* create work vectors for MatMult */
142:         PetscCall(VecDuplicate(ctx->x2,&ctx->wx2));
143:         PetscCall(VecDuplicate(ctx->y2,&ctx->wy2));
144:       }
145:     }
146: #endif
147:   }
148:   PetscFunctionReturn(PETSC_SUCCESS);
149: }

151: static PetscErrorCode MatMult_ECross(Mat B,Vec x,Vec y)
152: {
153:   SVD_CYCLIC_SHELL  *ctx;
154:   const PetscScalar *px;
155:   PetscScalar       *py;
156:   PetscInt          mn,m,n;

158:   PetscFunctionBegin;
159:   PetscCall(MatShellGetContext(B,&ctx));
160:   PetscCall(MatGetLocalSize(ctx->A,NULL,&n));
161:   PetscCall(VecGetLocalSize(y,&mn));
162:   m = mn-n;
163:   PetscCall(VecGetArrayRead(x,&px));
164:   PetscCall(VecGetArrayWrite(y,&py));
165:   PetscCall(VecPlaceArray(ctx->x1,px));
166:   PetscCall(VecPlaceArray(ctx->x2,px+m));
167:   PetscCall(VecPlaceArray(ctx->y1,py));
168:   PetscCall(VecPlaceArray(ctx->y2,py+m));
169:   PetscCall(VecCopy(ctx->x1,ctx->y1));
170:   PetscCall(MatMult(ctx->A,ctx->x2,ctx->w));
171:   PetscCall(MatMult(ctx->AT,ctx->w,ctx->y2));
172:   PetscCall(VecResetArray(ctx->x1));
173:   PetscCall(VecResetArray(ctx->x2));
174:   PetscCall(VecResetArray(ctx->y1));
175:   PetscCall(VecResetArray(ctx->y2));
176:   PetscCall(VecRestoreArrayRead(x,&px));
177:   PetscCall(VecRestoreArrayWrite(y,&py));
178:   PetscFunctionReturn(PETSC_SUCCESS);
179: }

181: static PetscErrorCode MatGetDiagonal_ECross(Mat B,Vec d)
182: {
183:   SVD_CYCLIC_SHELL  *ctx;
184:   PetscScalar       *pd;
185:   PetscMPIInt       len;
186:   PetscInt          mn,m,n,N,i,j,start,end,ncols;
187:   PetscScalar       *work1,*work2,*diag;
188:   const PetscInt    *cols;
189:   const PetscScalar *vals;

191:   PetscFunctionBegin;
192:   PetscCall(MatShellGetContext(B,&ctx));
193:   PetscCall(MatGetLocalSize(ctx->A,NULL,&n));
194:   PetscCall(VecGetLocalSize(d,&mn));
195:   m = mn-n;
196:   PetscCall(VecGetArrayWrite(d,&pd));
197:   PetscCall(VecPlaceArray(ctx->y1,pd));
198:   PetscCall(VecSet(ctx->y1,1.0));
199:   PetscCall(VecResetArray(ctx->y1));
200:   PetscCall(VecPlaceArray(ctx->y2,pd+m));
201:   if (!ctx->diag) {
202:     /* compute diagonal from rows and store in ctx->diag */
203:     PetscCall(VecDuplicate(ctx->y2,&ctx->diag));
204:     PetscCall(MatGetSize(ctx->A,NULL,&N));
205:     PetscCall(PetscCalloc2(N,&work1,N,&work2));
206:     if (ctx->swapped) {
207:       PetscCall(MatGetOwnershipRange(ctx->AT,&start,&end));
208:       for (i=start;i<end;i++) {
209:         PetscCall(MatGetRow(ctx->AT,i,&ncols,NULL,&vals));
210:         for (j=0;j<ncols;j++) work1[i] += vals[j]*vals[j];
211:         PetscCall(MatRestoreRow(ctx->AT,i,&ncols,NULL,&vals));
212:       }
213:     } else {
214:       PetscCall(MatGetOwnershipRange(ctx->A,&start,&end));
215:       for (i=start;i<end;i++) {
216:         PetscCall(MatGetRow(ctx->A,i,&ncols,&cols,&vals));
217:         for (j=0;j<ncols;j++) work1[cols[j]] += vals[j]*vals[j];
218:         PetscCall(MatRestoreRow(ctx->A,i,&ncols,&cols,&vals));
219:       }
220:     }
221:     PetscCall(PetscMPIIntCast(N,&len));
222:     PetscCall(MPIU_Allreduce(work1,work2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)B)));
223:     PetscCall(VecGetOwnershipRange(ctx->diag,&start,&end));
224:     PetscCall(VecGetArrayWrite(ctx->diag,&diag));
225:     for (i=start;i<end;i++) diag[i-start] = work2[i];
226:     PetscCall(VecRestoreArrayWrite(ctx->diag,&diag));
227:     PetscCall(PetscFree2(work1,work2));
228:   }
229:   PetscCall(VecCopy(ctx->diag,ctx->y2));
230:   PetscCall(VecResetArray(ctx->y2));
231:   PetscCall(VecRestoreArrayWrite(d,&pd));
232:   PetscFunctionReturn(PETSC_SUCCESS);
233: }

235: static PetscErrorCode MatDestroy_ECross(Mat B)
236: {
237:   SVD_CYCLIC_SHELL *ctx;

239:   PetscFunctionBegin;
240:   PetscCall(MatShellGetContext(B,&ctx));
241:   PetscCall(VecDestroy(&ctx->x1));
242:   PetscCall(VecDestroy(&ctx->x2));
243:   PetscCall(VecDestroy(&ctx->y1));
244:   PetscCall(VecDestroy(&ctx->y2));
245:   PetscCall(VecDestroy(&ctx->diag));
246:   PetscCall(VecDestroy(&ctx->w));
247:   if (ctx->misaligned) {
248:     PetscCall(VecDestroy(&ctx->wx2));
249:     PetscCall(VecDestroy(&ctx->wy2));
250:   }
251:   PetscCall(PetscFree(ctx));
252:   PetscFunctionReturn(PETSC_SUCCESS);
253: }

255: /*
256:    Builds extended cross product matrix   C = | I_m   0  |
257:                                               |  0  AT*A |
258:    t is an auxiliary Vec used to take the dimensions of the upper block
259: */
260: static PetscErrorCode SVDCyclicGetECrossMat(SVD svd,Mat A,Mat AT,Mat *C,Vec t)
261: {
262:   SVD_CYCLIC       *cyclic = (SVD_CYCLIC*)svd->data;
263:   SVD_CYCLIC_SHELL *ctx;
264:   PetscInt         i,M,N,m,n,Istart,Iend;
265:   VecType          vtype;
266:   Mat              Id,Zm,Zn,ATA;
267: #if defined(PETSC_HAVE_CUDA)
268:   PetscBool        cuda;
269:   const PetscInt   *ranges;
270:   PetscMPIInt      size;
271: #endif

273:   PetscFunctionBegin;
274:   PetscCall(MatGetSize(A,NULL,&N));
275:   PetscCall(MatGetLocalSize(A,NULL,&n));
276:   PetscCall(VecGetSize(t,&M));
277:   PetscCall(VecGetLocalSize(t,&m));

279:   if (cyclic->explicitmatrix) {
280:     PetscCheck(svd->expltrans,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Cannot use explicit cyclic matrix with implicit transpose");
281:     PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)svd),m,m,M,M,1.0,&Id));
282:     PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Zm));
283:     PetscCall(MatSetSizes(Zm,m,n,M,N));
284:     PetscCall(MatSetFromOptions(Zm));
285:     PetscCall(MatSetUp(Zm));
286:     PetscCall(MatGetOwnershipRange(Zm,&Istart,&Iend));
287:     for (i=Istart;i<Iend;i++) {
288:       if (i<N) PetscCall(MatSetValue(Zm,i,i,0.0,INSERT_VALUES));
289:     }
290:     PetscCall(MatAssemblyBegin(Zm,MAT_FINAL_ASSEMBLY));
291:     PetscCall(MatAssemblyEnd(Zm,MAT_FINAL_ASSEMBLY));
292:     PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Zn));
293:     PetscCall(MatSetSizes(Zn,n,m,N,M));
294:     PetscCall(MatSetFromOptions(Zn));
295:     PetscCall(MatSetUp(Zn));
296:     PetscCall(MatGetOwnershipRange(Zn,&Istart,&Iend));
297:     for (i=Istart;i<Iend;i++) {
298:       if (i<m) PetscCall(MatSetValue(Zn,i,i,0.0,INSERT_VALUES));
299:     }
300:     PetscCall(MatAssemblyBegin(Zn,MAT_FINAL_ASSEMBLY));
301:     PetscCall(MatAssemblyEnd(Zn,MAT_FINAL_ASSEMBLY));
302:     PetscCall(MatProductCreate(AT,A,NULL,&ATA));
303:     PetscCall(MatProductSetType(ATA,MATPRODUCT_AB));
304:     PetscCall(MatProductSetFromOptions(ATA));
305:     PetscCall(MatProductSymbolic(ATA));
306:     PetscCall(MatProductNumeric(ATA));
307:     PetscCall(MatCreateTile(1.0,Id,1.0,Zm,1.0,Zn,1.0,ATA,C));
308:     PetscCall(MatDestroy(&Id));
309:     PetscCall(MatDestroy(&Zm));
310:     PetscCall(MatDestroy(&Zn));
311:     PetscCall(MatDestroy(&ATA));
312:   } else {
313:     PetscCall(PetscNew(&ctx));
314:     ctx->A       = A;
315:     ctx->AT      = AT;
316:     ctx->swapped = svd->swapped;
317:     PetscCall(VecDuplicateEmpty(t,&ctx->x1));
318:     PetscCall(VecDuplicateEmpty(t,&ctx->y1));
319:     PetscCall(MatCreateVecsEmpty(A,&ctx->x2,NULL));
320:     PetscCall(MatCreateVecsEmpty(A,&ctx->y2,NULL));
321:     PetscCall(MatCreateVecs(A,NULL,&ctx->w));
322:     PetscCall(MatCreateShell(PetscObjectComm((PetscObject)svd),m+n,m+n,M+N,M+N,ctx,C));
323:     PetscCall(MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ECross));
324:     PetscCall(MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_ECross));
325: #if defined(PETSC_HAVE_CUDA)
326:     PetscCall(PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&cuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,""));
327:     if (cuda) PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross_CUDA));
328:     else
329: #endif
330:       PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross));
331:     PetscCall(MatGetVecType(A,&vtype));
332:     PetscCall(MatSetVecType(*C,vtype));
333: #if defined(PETSC_HAVE_CUDA)
334:     if (cuda) {
335:       /* check alignment of bottom block */
336:       PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ctx->x1),&size));
337:       PetscCall(VecGetOwnershipRanges(ctx->x1,&ranges));
338:       for (i=0;i<size;i++) {
339:         ctx->misaligned = (((ranges[i+1]-ranges[i])*sizeof(PetscScalar))%16)? PETSC_TRUE: PETSC_FALSE;
340:         if (ctx->misaligned) break;
341:       }
342:       if (ctx->misaligned) {  /* create work vectors for MatMult */
343:         PetscCall(VecDuplicate(ctx->x2,&ctx->wx2));
344:         PetscCall(VecDuplicate(ctx->y2,&ctx->wy2));
345:       }
346:     }
347: #endif
348:   }
349:   PetscFunctionReturn(PETSC_SUCCESS);
350: }

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

357:   PetscFunctionBegin;
358:   *errest = res/PetscMax(svd->nrma,svd->nrmb);
359:   PetscFunctionReturn(PETSC_SUCCESS);
360: }

362: static PetscErrorCode SVDSetUp_Cyclic(SVD svd)
363: {
364:   SVD_CYCLIC        *cyclic = (SVD_CYCLIC*)svd->data;
365:   PetscInt          M,N,m,n,p,k,i,isl,offset,nev,ncv,mpd,maxit;
366:   PetscReal         tol;
367:   const PetscScalar *isa,*oa;
368:   PetscScalar       *va;
369:   EPSProblemType    ptype;
370:   PetscBool         trackall,issinv;
371:   Vec               v,t;
372:   ST                st;
373:   Mat               Omega;
374:   MatType           Atype;

376:   PetscFunctionBegin;
377:   PetscCall(MatGetSize(svd->A,&M,&N));
378:   PetscCall(MatGetLocalSize(svd->A,&m,&n));
379:   if (!cyclic->eps) PetscCall(SVDCyclicGetEPS(svd,&cyclic->eps));
380:   PetscCall(MatDestroy(&cyclic->C));
381:   PetscCall(MatDestroy(&cyclic->D));
382:   if (svd->isgeneralized) {
383:     if (svd->which==SVD_SMALLEST) {  /* alternative pencil */
384:       PetscCall(MatCreateVecs(svd->B,NULL,&t));
385:       PetscCall(SVDCyclicGetCyclicMat(svd,svd->B,svd->BT,&cyclic->C));
386:       PetscCall(SVDCyclicGetECrossMat(svd,svd->A,svd->AT,&cyclic->D,t));
387:     } else {
388:       PetscCall(MatCreateVecs(svd->A,NULL,&t));
389:       PetscCall(SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C));
390:       PetscCall(SVDCyclicGetECrossMat(svd,svd->B,svd->BT,&cyclic->D,t));
391:     }
392:     PetscCall(VecDestroy(&t));
393:     PetscCall(EPSSetOperators(cyclic->eps,cyclic->C,cyclic->D));
394:     PetscCall(EPSGetProblemType(cyclic->eps,&ptype));
395:     if (!ptype) PetscCall(EPSSetProblemType(cyclic->eps,EPS_GHEP));
396:   } else if (svd->ishyperbolic) {
397:     PetscCall(SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C));
398:     PetscCall(MatCreateVecs(cyclic->C,&v,NULL));
399:     PetscCall(VecSet(v,1.0));
400:     PetscCall(VecGetArrayRead(svd->omega,&oa));
401:     PetscCall(VecGetArray(v,&va));
402:     if (svd->swapped) PetscCall(PetscArraycpy(va+m,oa,n));
403:     else PetscCall(PetscArraycpy(va,oa,m));
404:     PetscCall(VecRestoreArrayRead(svd->omega,&oa));
405:     PetscCall(VecRestoreArray(v,&va));
406:     PetscCall(MatGetType(svd->OP,&Atype));
407:     PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Omega));
408:     PetscCall(MatSetSizes(Omega,m+n,m+n,M+N,M+N));
409:     PetscCall(MatSetType(Omega,Atype));
410:     PetscCall(MatSetUp(Omega));
411:     PetscCall(MatDiagonalSet(Omega,v,INSERT_VALUES));
412:     PetscCall(EPSSetOperators(cyclic->eps,cyclic->C,Omega));
413:     PetscCall(EPSSetProblemType(cyclic->eps,EPS_GHIEP));
414:     PetscCall(MatDestroy(&Omega));
415:     PetscCall(VecDestroy(&v));
416:   } else {
417:     PetscCall(SVDCyclicGetCyclicMat(svd,svd->A,svd->AT,&cyclic->C));
418:     PetscCall(EPSSetOperators(cyclic->eps,cyclic->C,NULL));
419:     PetscCall(EPSSetProblemType(cyclic->eps,EPS_HEP));
420:   }
421:   if (!cyclic->usereps) {
422:     if (svd->which == SVD_LARGEST) {
423:       PetscCall(EPSGetST(cyclic->eps,&st));
424:       PetscCall(PetscObjectTypeCompare((PetscObject)st,STSINVERT,&issinv));
425:       if (issinv) PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_TARGET_MAGNITUDE));
426:       else if (svd->ishyperbolic) PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_MAGNITUDE));
427:       else PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL));
428:     } else {
429:       if (svd->isgeneralized) {  /* computes sigma^{-1} via alternative pencil */
430:         PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL));
431:       } else {
432:         if (svd->ishyperbolic) PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_TARGET_MAGNITUDE));
433:         else PetscCall(EPSSetEigenvalueComparison(cyclic->eps,SlepcCompareSmallestPosReal,NULL));
434:         PetscCall(EPSSetTarget(cyclic->eps,0.0));
435:       }
436:     }
437:     PetscCall(EPSGetDimensions(cyclic->eps,&nev,&ncv,&mpd));
438:     PetscCheck(nev==1 || nev>=2*svd->nsv,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONG,"The number of requested eigenvalues %" PetscInt_FMT " must be at least 2*%" PetscInt_FMT,nev,svd->nsv);
439:     nev = PetscMax(nev,2*svd->nsv);
440:     if (ncv==PETSC_DEFAULT && svd->ncv!=PETSC_DEFAULT) ncv = PetscMax(3*svd->nsv,svd->ncv);
441:     if (mpd==PETSC_DEFAULT && svd->mpd!=PETSC_DEFAULT) mpd = svd->mpd;
442:     PetscCall(EPSSetDimensions(cyclic->eps,nev,ncv,mpd));
443:     PetscCall(EPSGetTolerances(cyclic->eps,&tol,&maxit));
444:     if (tol==(PetscReal)PETSC_DEFAULT) tol = svd->tol==(PetscReal)PETSC_DEFAULT? SLEPC_DEFAULT_TOL/10.0: svd->tol;
445:     if (maxit==PETSC_DEFAULT && svd->max_it!=PETSC_DEFAULT) maxit = svd->max_it;
446:     PetscCall(EPSSetTolerances(cyclic->eps,tol,maxit));
447:     switch (svd->conv) {
448:     case SVD_CONV_ABS:
449:       PetscCall(EPSSetConvergenceTest(cyclic->eps,EPS_CONV_ABS));break;
450:     case SVD_CONV_REL:
451:       PetscCall(EPSSetConvergenceTest(cyclic->eps,EPS_CONV_REL));break;
452:     case SVD_CONV_NORM:
453:       if (svd->isgeneralized) {
454:         if (!svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
455:         if (!svd->nrmb) PetscCall(MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb));
456:         PetscCall(EPSSetConvergenceTestFunction(cyclic->eps,EPSConv_Cyclic,svd,NULL));
457:       } else {
458:         PetscCall(EPSSetConvergenceTest(cyclic->eps,EPS_CONV_NORM));break;
459:       }
460:       break;
461:     case SVD_CONV_MAXIT:
462:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Maxit convergence test not supported in this solver");
463:     case SVD_CONV_USER:
464:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
465:     }
466:   }
467:   SVDCheckUnsupported(svd,SVD_FEATURE_STOPPING);
468:   /* Transfer the trackall option from svd to eps */
469:   PetscCall(SVDGetTrackAll(svd,&trackall));
470:   PetscCall(EPSSetTrackAll(cyclic->eps,trackall));
471:   /* Transfer the initial subspace from svd to eps */
472:   if (svd->nini<0 || svd->ninil<0) {
473:     for (i=0;i<-PetscMin(svd->nini,svd->ninil);i++) {
474:       PetscCall(MatCreateVecs(cyclic->C,&v,NULL));
475:       PetscCall(VecGetArrayWrite(v,&va));
476:       if (svd->isgeneralized) PetscCall(MatGetLocalSize(svd->B,&p,NULL));
477:       k = (svd->isgeneralized && svd->which==SVD_SMALLEST)? p: m;  /* size of upper block row */
478:       if (i<-svd->ninil) {
479:         PetscCall(VecGetArrayRead(svd->ISL[i],&isa));
480:         if (svd->isgeneralized) {
481:           PetscCall(VecGetLocalSize(svd->ISL[i],&isl));
482:           PetscCheck(isl==m+p,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Size mismatch for left initial vector");
483:           offset = (svd->which==SVD_SMALLEST)? m: 0;
484:           PetscCall(PetscArraycpy(va,isa+offset,k));
485:         } else {
486:           PetscCall(VecGetLocalSize(svd->ISL[i],&isl));
487:           PetscCheck(isl==k,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Size mismatch for left initial vector");
488:           PetscCall(PetscArraycpy(va,isa,k));
489:         }
490:         PetscCall(VecRestoreArrayRead(svd->IS[i],&isa));
491:       } else PetscCall(PetscArrayzero(&va,k));
492:       if (i<-svd->nini) {
493:         PetscCall(VecGetLocalSize(svd->IS[i],&isl));
494:         PetscCheck(isl==n,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Size mismatch for right initial vector");
495:         PetscCall(VecGetArrayRead(svd->IS[i],&isa));
496:         PetscCall(PetscArraycpy(va+k,isa,n));
497:         PetscCall(VecRestoreArrayRead(svd->IS[i],&isa));
498:       } else PetscCall(PetscArrayzero(va+k,n));
499:       PetscCall(VecRestoreArrayWrite(v,&va));
500:       PetscCall(VecDestroy(&svd->IS[i]));
501:       svd->IS[i] = v;
502:     }
503:     svd->nini = PetscMin(svd->nini,svd->ninil);
504:     PetscCall(EPSSetInitialSpace(cyclic->eps,-svd->nini,svd->IS));
505:     PetscCall(SlepcBasisDestroy_Private(&svd->nini,&svd->IS));
506:     PetscCall(SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL));
507:   }
508:   PetscCall(EPSSetUp(cyclic->eps));
509:   PetscCall(EPSGetDimensions(cyclic->eps,NULL,&svd->ncv,&svd->mpd));
510:   svd->ncv = PetscMin(svd->ncv,PetscMin(M,N));
511:   PetscCall(EPSGetTolerances(cyclic->eps,NULL,&svd->max_it));
512:   if (svd->tol==(PetscReal)PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;

514:   svd->leftbasis = PETSC_TRUE;
515:   PetscCall(SVDAllocateSolution(svd,0));
516:   PetscFunctionReturn(PETSC_SUCCESS);
517: }

519: static PetscErrorCode SVDCyclicCheckEigenvalue(SVD svd,PetscScalar er,PetscScalar ei,PetscReal *sigma,PetscBool *isreal)
520: {
521:   PetscFunctionBegin;
522:   if (svd->ishyperbolic && PetscDefined(USE_COMPLEX) && PetscAbsReal(PetscImaginaryPart(er))>10*PetscAbsReal(PetscRealPart(er))) {
523:     *sigma = PetscImaginaryPart(er);
524:     if (isreal) *isreal = PETSC_FALSE;
525:   } else if (svd->ishyperbolic && !PetscDefined(USE_COMPLEX) && PetscAbsScalar(ei)>10*PetscAbsScalar(er)) {
526:     *sigma = PetscRealPart(ei);
527:     if (isreal) *isreal = PETSC_FALSE;
528:   } else {
529:     *sigma = PetscRealPart(er);
530:     if (isreal) *isreal = PETSC_TRUE;
531:   }
532:   PetscFunctionReturn(PETSC_SUCCESS);
533: }

535: static PetscErrorCode SVDSolve_Cyclic(SVD svd)
536: {
537:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;
538:   PetscInt       i,j,nconv;
539:   PetscScalar    er,ei;
540:   PetscReal      sigma;

542:   PetscFunctionBegin;
543:   PetscCall(EPSSolve(cyclic->eps));
544:   PetscCall(EPSGetConverged(cyclic->eps,&nconv));
545:   PetscCall(EPSGetIterationNumber(cyclic->eps,&svd->its));
546:   PetscCall(EPSGetConvergedReason(cyclic->eps,(EPSConvergedReason*)&svd->reason));
547:   for (i=0,j=0;i<nconv;i++) {
548:     PetscCall(EPSGetEigenvalue(cyclic->eps,i,&er,&ei));
549:     PetscCall(SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL));
550:     if (sigma>0.0) {
551:       if (svd->isgeneralized && svd->which==SVD_SMALLEST) svd->sigma[j] = 1.0/sigma;
552:       else svd->sigma[j] = sigma;
553:       j++;
554:     }
555:   }
556:   svd->nconv = j;
557:   PetscFunctionReturn(PETSC_SUCCESS);
558: }

560: static PetscErrorCode SVDComputeVectors_Cyclic_Standard(SVD svd)
561: {
562:   SVD_CYCLIC        *cyclic = (SVD_CYCLIC*)svd->data;
563:   PetscInt          i,j,m,nconv;
564:   PetscScalar       er,ei;
565:   PetscReal         sigma;
566:   const PetscScalar *px;
567:   Vec               x,x1,x2;

569:   PetscFunctionBegin;
570:   PetscCall(MatCreateVecs(cyclic->C,&x,NULL));
571:   PetscCall(MatGetLocalSize(svd->A,&m,NULL));
572:   PetscCall(MatCreateVecsEmpty(svd->A,&x2,&x1));
573:   PetscCall(EPSGetConverged(cyclic->eps,&nconv));
574:   for (i=0,j=0;i<nconv;i++) {
575:     PetscCall(EPSGetEigenpair(cyclic->eps,i,&er,&ei,x,NULL));
576:     PetscCall(SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL));
577:     if (sigma<0.0) continue;
578:     PetscCall(VecGetArrayRead(x,&px));
579:     PetscCall(VecPlaceArray(x1,px));
580:     PetscCall(VecPlaceArray(x2,px+m));
581:     PetscCall(BVInsertVec(svd->U,j,x1));
582:     PetscCall(BVScaleColumn(svd->U,j,PETSC_SQRT2));
583:     PetscCall(BVInsertVec(svd->V,j,x2));
584:     PetscCall(BVScaleColumn(svd->V,j,PETSC_SQRT2));
585:     PetscCall(VecResetArray(x1));
586:     PetscCall(VecResetArray(x2));
587:     PetscCall(VecRestoreArrayRead(x,&px));
588:     j++;
589:   }
590:   PetscCall(VecDestroy(&x));
591:   PetscCall(VecDestroy(&x1));
592:   PetscCall(VecDestroy(&x2));
593:   PetscFunctionReturn(PETSC_SUCCESS);
594: }

596: static PetscErrorCode SVDComputeVectors_Cyclic_Generalized(SVD svd)
597: {
598:   SVD_CYCLIC        *cyclic = (SVD_CYCLIC*)svd->data;
599:   PetscInt          i,j,m,p,nconv;
600:   PetscScalar       *dst,er,ei;
601:   PetscReal         sigma;
602:   const PetscScalar *src,*px;
603:   Vec               u,v,x,x1,x2,uv;

605:   PetscFunctionBegin;
606:   PetscCall(MatGetLocalSize(svd->A,&m,NULL));
607:   PetscCall(MatGetLocalSize(svd->B,&p,NULL));
608:   PetscCall(MatCreateVecs(cyclic->C,&x,NULL));
609:   if (svd->which==SVD_SMALLEST) PetscCall(MatCreateVecsEmpty(svd->B,&x1,&x2));
610:   else PetscCall(MatCreateVecsEmpty(svd->A,&x2,&x1));
611:   PetscCall(MatCreateVecs(svd->A,NULL,&u));
612:   PetscCall(MatCreateVecs(svd->B,NULL,&v));
613:   PetscCall(EPSGetConverged(cyclic->eps,&nconv));
614:   for (i=0,j=0;i<nconv;i++) {
615:     PetscCall(EPSGetEigenpair(cyclic->eps,i,&er,&ei,x,NULL));
616:     PetscCall(SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL));
617:     if (sigma<0.0) continue;
618:     if (svd->which==SVD_SMALLEST) {
619:       /* evec_i = 1/sqrt(2)*[ v_i; w_i ],  w_i = x_i/c_i */
620:       PetscCall(VecGetArrayRead(x,&px));
621:       PetscCall(VecPlaceArray(x2,px));
622:       PetscCall(VecPlaceArray(x1,px+p));
623:       PetscCall(VecCopy(x2,v));
624:       PetscCall(VecScale(v,PETSC_SQRT2));  /* v_i = sqrt(2)*evec_i_1 */
625:       PetscCall(VecScale(x1,PETSC_SQRT2)); /* w_i = sqrt(2)*evec_i_2 */
626:       PetscCall(MatMult(svd->A,x1,u));     /* A*w_i = u_i */
627:       PetscCall(VecScale(x1,1.0/PetscSqrtScalar(1.0+sigma*sigma)));  /* x_i = w_i*c_i */
628:       PetscCall(BVInsertVec(svd->V,j,x1));
629:       PetscCall(VecResetArray(x2));
630:       PetscCall(VecResetArray(x1));
631:       PetscCall(VecRestoreArrayRead(x,&px));
632:     } else {
633:       /* evec_i = 1/sqrt(2)*[ u_i; w_i ],  w_i = x_i/s_i */
634:       PetscCall(VecGetArrayRead(x,&px));
635:       PetscCall(VecPlaceArray(x1,px));
636:       PetscCall(VecPlaceArray(x2,px+m));
637:       PetscCall(VecCopy(x1,u));
638:       PetscCall(VecScale(u,PETSC_SQRT2));  /* u_i = sqrt(2)*evec_i_1 */
639:       PetscCall(VecScale(x2,PETSC_SQRT2)); /* w_i = sqrt(2)*evec_i_2 */
640:       PetscCall(MatMult(svd->B,x2,v));     /* B*w_i = v_i */
641:       PetscCall(VecScale(x2,1.0/PetscSqrtScalar(1.0+sigma*sigma)));  /* x_i = w_i*s_i */
642:       PetscCall(BVInsertVec(svd->V,j,x2));
643:       PetscCall(VecResetArray(x1));
644:       PetscCall(VecResetArray(x2));
645:       PetscCall(VecRestoreArrayRead(x,&px));
646:     }
647:     /* copy [u;v] to U[j] */
648:     PetscCall(BVGetColumn(svd->U,j,&uv));
649:     PetscCall(VecGetArrayWrite(uv,&dst));
650:     PetscCall(VecGetArrayRead(u,&src));
651:     PetscCall(PetscArraycpy(dst,src,m));
652:     PetscCall(VecRestoreArrayRead(u,&src));
653:     PetscCall(VecGetArrayRead(v,&src));
654:     PetscCall(PetscArraycpy(dst+m,src,p));
655:     PetscCall(VecRestoreArrayRead(v,&src));
656:     PetscCall(VecRestoreArrayWrite(uv,&dst));
657:     PetscCall(BVRestoreColumn(svd->U,j,&uv));
658:     j++;
659:   }
660:   PetscCall(VecDestroy(&x));
661:   PetscCall(VecDestroy(&x1));
662:   PetscCall(VecDestroy(&x2));
663:   PetscCall(VecDestroy(&u));
664:   PetscCall(VecDestroy(&v));
665:   PetscFunctionReturn(PETSC_SUCCESS);
666: }

668: #if defined(PETSC_USE_COMPLEX)
669: /* VecMaxAbs: returns the entry of x that has max(abs(x(i))), using w as a workspace vector */
670: static PetscErrorCode VecMaxAbs(Vec x,Vec w,PetscScalar *v)
671: {
672:   PetscMPIInt       size,rank,root;
673:   const PetscScalar *xx;
674:   const PetscInt    *ranges;
675:   PetscReal         val;
676:   PetscInt          p;

678:   PetscFunctionBegin;
679:   PetscCall(VecCopy(x,w));
680:   PetscCall(VecAbs(w));
681:   PetscCall(VecMax(w,&p,&val));
682:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)x),&size));
683:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)x),&rank));
684:   PetscCall(VecGetOwnershipRanges(x,&ranges));
685:   for (root=0;root<size;root++) if (p>=ranges[root] && p<ranges[root+1]) break;
686:   if (rank==root) {
687:     PetscCall(VecGetArrayRead(x,&xx));
688:     *v = xx[p-ranges[root]];
689:     PetscCall(VecRestoreArrayRead(x,&xx));
690:   }
691:   PetscCallMPI(MPI_Bcast(v,1,MPIU_SCALAR,root,PetscObjectComm((PetscObject)x)));
692:   PetscFunctionReturn(PETSC_SUCCESS);
693: }
694: #endif

696: static PetscErrorCode SVDComputeVectors_Cyclic_Hyperbolic(SVD svd)
697: {
698:   SVD_CYCLIC        *cyclic = (SVD_CYCLIC*)svd->data;
699:   PetscInt          i,j,m,n,nconv;
700:   PetscScalar       er,ei;
701:   PetscReal         sigma,nrm;
702:   PetscBool         isreal;
703:   const PetscScalar *px;
704:   Vec               u,x,xi=NULL,x1,x2,x1i=NULL,x2i;
705:   BV                U=NULL,V=NULL;
706: #if !defined(PETSC_USE_COMPLEX)
707:   const PetscScalar *pxi;
708:   PetscReal         nrmr,nrmi;
709: #else
710:   PetscScalar       alpha;
711: #endif

713:   PetscFunctionBegin;
714:   PetscCall(MatCreateVecs(cyclic->C,&x,svd->ishyperbolic?&xi:NULL));
715:   PetscCall(MatGetLocalSize(svd->A,&m,NULL));
716:   PetscCall(MatCreateVecsEmpty(svd->OP,&x2,&x1));
717: #if defined(PETSC_USE_COMPLEX)
718:   PetscCall(MatCreateVecs(svd->OP,&x2i,&x1i));
719: #else
720:   PetscCall(MatCreateVecsEmpty(svd->OP,&x2i,&x1i));
721: #endif
722:   /* set-up Omega-normalization of U */
723:   U = svd->swapped? svd->V: svd->U;
724:   V = svd->swapped? svd->U: svd->V;
725:   PetscCall(BVGetSizes(U,&n,NULL,NULL));
726:   PetscCall(BV_SetMatrixDiagonal(U,svd->omega,svd->A));
727:   PetscCall(EPSGetConverged(cyclic->eps,&nconv));
728:   for (i=0,j=0;i<nconv;i++) {
729:     PetscCall(EPSGetEigenpair(cyclic->eps,i,&er,&ei,x,xi));
730:     PetscCall(SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,&isreal));
731:     if (sigma<0.0) continue;
732:     PetscCall(VecGetArrayRead(x,&px));
733:     if (svd->swapped) {
734:       PetscCall(VecPlaceArray(x2,px));
735:       PetscCall(VecPlaceArray(x1,px+m));
736:     } else {
737:       PetscCall(VecPlaceArray(x1,px));
738:       PetscCall(VecPlaceArray(x2,px+n));
739:     }
740: #if defined(PETSC_USE_COMPLEX)
741:     PetscCall(BVInsertVec(U,j,x1));
742:     PetscCall(BVInsertVec(V,j,x2));
743:     if (!isreal) {
744:       PetscCall(VecMaxAbs(x1,x1i,&alpha));
745:       PetscCall(BVScaleColumn(U,j,PetscAbsScalar(alpha)/alpha));
746:       PetscCall(BVScaleColumn(V,j,PetscAbsScalar(alpha)/(alpha*PETSC_i)));
747:     }
748: #else
749:     PetscCall(VecGetArrayRead(xi,&pxi));
750:     if (svd->swapped) {
751:       PetscCall(VecPlaceArray(x2i,pxi));
752:       PetscCall(VecPlaceArray(x1i,pxi+m));
753:     } else {
754:       PetscCall(VecPlaceArray(x1i,pxi));
755:       PetscCall(VecPlaceArray(x2i,pxi+n));
756:     }
757:     PetscCall(VecNorm(x2,NORM_2,&nrmr));
758:     PetscCall(VecNorm(x2i,NORM_2,&nrmi));
759:     if (nrmi>nrmr) {
760:       if (isreal) {
761:         PetscCall(BVInsertVec(U,j,x1i));
762:         PetscCall(BVInsertVec(V,j,x2i));
763:       } else {
764:         PetscCall(BVInsertVec(U,j,x1));
765:         PetscCall(BVInsertVec(V,j,x2i));
766:       }
767:     } else {
768:       if (isreal) {
769:         PetscCall(BVInsertVec(U,j,x1));
770:         PetscCall(BVInsertVec(V,j,x2));
771:       } else {
772:         PetscCall(BVInsertVec(U,j,x1i));
773:         PetscCall(BVScaleColumn(U,j,-1.0));
774:         PetscCall(BVInsertVec(V,j,x2));
775:       }
776:     }
777:     PetscCall(VecResetArray(x1i));
778:     PetscCall(VecResetArray(x2i));
779:     PetscCall(VecRestoreArrayRead(xi,&pxi));
780: #endif
781:     PetscCall(VecResetArray(x1));
782:     PetscCall(VecResetArray(x2));
783:     PetscCall(VecRestoreArrayRead(x,&px));
784:     PetscCall(BVGetColumn(U,j,&u));
785:     PetscCall(VecPointwiseMult(u,u,svd->omega));
786:     PetscCall(BVRestoreColumn(U,j,&u));
787:     PetscCall(BVNormColumn(U,j,NORM_2,&nrm));
788:     PetscCall(BVScaleColumn(U,j,1.0/PetscAbs(nrm)));
789:     svd->sign[j] = PetscSign(nrm);
790:     PetscCall(BVNormColumn(V,j,NORM_2,&nrm));
791:     PetscCall(BVScaleColumn(V,j,1.0/nrm));
792:     j++;
793:   }
794:   PetscCall(VecDestroy(&x));
795:   PetscCall(VecDestroy(&x1));
796:   PetscCall(VecDestroy(&x2));
797:   PetscCall(VecDestroy(&xi));
798:   PetscCall(VecDestroy(&x1i));
799:   PetscCall(VecDestroy(&x2i));
800:   PetscFunctionReturn(PETSC_SUCCESS);
801: }

803: static PetscErrorCode SVDComputeVectors_Cyclic(SVD svd)
804: {
805:   PetscFunctionBegin;
806:   switch (svd->problem_type) {
807:     case SVD_STANDARD:
808:       PetscCall(SVDComputeVectors_Cyclic_Standard(svd));
809:       break;
810:     case SVD_GENERALIZED:
811:       PetscCall(SVDComputeVectors_Cyclic_Generalized(svd));
812:       break;
813:     case SVD_HYPERBOLIC:
814:       PetscCall(SVDComputeVectors_Cyclic_Hyperbolic(svd));
815:       break;
816:     default:
817:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_WRONG,"Unknown singular value problem type");
818:   }
819:   PetscFunctionReturn(PETSC_SUCCESS);
820: }

822: static PetscErrorCode EPSMonitor_Cyclic(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
823: {
824:   PetscInt       i,j;
825:   SVD            svd = (SVD)ctx;
826:   PetscScalar    er,ei;
827:   PetscReal      sigma;
828:   ST             st;

830:   PetscFunctionBegin;
831:   nconv = 0;
832:   PetscCall(EPSGetST(eps,&st));
833:   for (i=0,j=0;i<PetscMin(nest,svd->ncv);i++) {
834:     er = eigr[i]; ei = eigi[i];
835:     PetscCall(STBackTransform(st,1,&er,&ei));
836:     PetscCall(SVDCyclicCheckEigenvalue(svd,er,ei,&sigma,NULL));
837:     if (sigma>0.0) {
838:       svd->sigma[j]  = sigma;
839:       svd->errest[j] = errest[i];
840:       if (errest[i] && errest[i] < svd->tol) nconv++;
841:       j++;
842:     }
843:   }
844:   nest = j;
845:   PetscCall(SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest));
846:   PetscFunctionReturn(PETSC_SUCCESS);
847: }

849: static PetscErrorCode SVDSetFromOptions_Cyclic(SVD svd,PetscOptionItems *PetscOptionsObject)
850: {
851:   PetscBool      set,val;
852:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;
853:   ST             st;

855:   PetscFunctionBegin;
856:   PetscOptionsHeadBegin(PetscOptionsObject,"SVD Cyclic Options");

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

861:   PetscOptionsHeadEnd();

863:   if (!cyclic->eps) PetscCall(SVDCyclicGetEPS(svd,&cyclic->eps));
864:   if (!cyclic->explicitmatrix && !cyclic->usereps) {
865:     /* use as default an ST with shell matrix and Jacobi */
866:     PetscCall(EPSGetST(cyclic->eps,&st));
867:     PetscCall(STSetMatMode(st,ST_MATMODE_SHELL));
868:   }
869:   PetscCall(EPSSetFromOptions(cyclic->eps));
870:   PetscFunctionReturn(PETSC_SUCCESS);
871: }

873: static PetscErrorCode SVDCyclicSetExplicitMatrix_Cyclic(SVD svd,PetscBool explicitmat)
874: {
875:   SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;

877:   PetscFunctionBegin;
878:   if (cyclic->explicitmatrix != explicitmat) {
879:     cyclic->explicitmatrix = explicitmat;
880:     svd->state = SVD_STATE_INITIAL;
881:   }
882:   PetscFunctionReturn(PETSC_SUCCESS);
883: }

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

889:    Logically Collective

891:    Input Parameters:
892: +  svd         - singular value solver
893: -  explicitmat - boolean flag indicating if H(A) is built explicitly

895:    Options Database Key:
896: .  -svd_cyclic_explicitmatrix <boolean> - Indicates the boolean flag

898:    Level: advanced

900: .seealso: SVDCyclicGetExplicitMatrix()
901: @*/
902: PetscErrorCode SVDCyclicSetExplicitMatrix(SVD svd,PetscBool explicitmat)
903: {
904:   PetscFunctionBegin;
907:   PetscTryMethod(svd,"SVDCyclicSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
908:   PetscFunctionReturn(PETSC_SUCCESS);
909: }

911: static PetscErrorCode SVDCyclicGetExplicitMatrix_Cyclic(SVD svd,PetscBool *explicitmat)
912: {
913:   SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;

915:   PetscFunctionBegin;
916:   *explicitmat = cyclic->explicitmatrix;
917:   PetscFunctionReturn(PETSC_SUCCESS);
918: }

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

923:    Not Collective

925:    Input Parameter:
926: .  svd  - singular value solver

928:    Output Parameter:
929: .  explicitmat - the mode flag

931:    Level: advanced

933: .seealso: SVDCyclicSetExplicitMatrix()
934: @*/
935: PetscErrorCode SVDCyclicGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
936: {
937:   PetscFunctionBegin;
939:   PetscAssertPointer(explicitmat,2);
940:   PetscUseMethod(svd,"SVDCyclicGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
941:   PetscFunctionReturn(PETSC_SUCCESS);
942: }

944: static PetscErrorCode SVDCyclicSetEPS_Cyclic(SVD svd,EPS eps)
945: {
946:   SVD_CYCLIC      *cyclic = (SVD_CYCLIC*)svd->data;

948:   PetscFunctionBegin;
949:   PetscCall(PetscObjectReference((PetscObject)eps));
950:   PetscCall(EPSDestroy(&cyclic->eps));
951:   cyclic->eps     = eps;
952:   cyclic->usereps = PETSC_TRUE;
953:   svd->state      = SVD_STATE_INITIAL;
954:   PetscFunctionReturn(PETSC_SUCCESS);
955: }

957: /*@
958:    SVDCyclicSetEPS - Associate an eigensolver object (EPS) to the
959:    singular value solver.

961:    Collective

963:    Input Parameters:
964: +  svd - singular value solver
965: -  eps - the eigensolver object

967:    Level: advanced

969: .seealso: SVDCyclicGetEPS()
970: @*/
971: PetscErrorCode SVDCyclicSetEPS(SVD svd,EPS eps)
972: {
973:   PetscFunctionBegin;
976:   PetscCheckSameComm(svd,1,eps,2);
977:   PetscTryMethod(svd,"SVDCyclicSetEPS_C",(SVD,EPS),(svd,eps));
978:   PetscFunctionReturn(PETSC_SUCCESS);
979: }

981: static PetscErrorCode SVDCyclicGetEPS_Cyclic(SVD svd,EPS *eps)
982: {
983:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;

985:   PetscFunctionBegin;
986:   if (!cyclic->eps) {
987:     PetscCall(EPSCreate(PetscObjectComm((PetscObject)svd),&cyclic->eps));
988:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)cyclic->eps,(PetscObject)svd,1));
989:     PetscCall(EPSSetOptionsPrefix(cyclic->eps,((PetscObject)svd)->prefix));
990:     PetscCall(EPSAppendOptionsPrefix(cyclic->eps,"svd_cyclic_"));
991:     PetscCall(PetscObjectSetOptions((PetscObject)cyclic->eps,((PetscObject)svd)->options));
992:     PetscCall(EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL));
993:     PetscCall(EPSMonitorSet(cyclic->eps,EPSMonitor_Cyclic,svd,NULL));
994:   }
995:   *eps = cyclic->eps;
996:   PetscFunctionReturn(PETSC_SUCCESS);
997: }

999: /*@
1000:    SVDCyclicGetEPS - Retrieve the eigensolver object (EPS) associated
1001:    to the singular value solver.

1003:    Collective

1005:    Input Parameter:
1006: .  svd - singular value solver

1008:    Output Parameter:
1009: .  eps - the eigensolver object

1011:    Level: advanced

1013: .seealso: SVDCyclicSetEPS()
1014: @*/
1015: PetscErrorCode SVDCyclicGetEPS(SVD svd,EPS *eps)
1016: {
1017:   PetscFunctionBegin;
1019:   PetscAssertPointer(eps,2);
1020:   PetscUseMethod(svd,"SVDCyclicGetEPS_C",(SVD,EPS*),(svd,eps));
1021:   PetscFunctionReturn(PETSC_SUCCESS);
1022: }

1024: static PetscErrorCode SVDView_Cyclic(SVD svd,PetscViewer viewer)
1025: {
1026:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;
1027:   PetscBool      isascii;

1029:   PetscFunctionBegin;
1030:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1031:   if (isascii) {
1032:     if (!cyclic->eps) PetscCall(SVDCyclicGetEPS(svd,&cyclic->eps));
1033:     PetscCall(PetscViewerASCIIPrintf(viewer,"  %s matrix\n",cyclic->explicitmatrix?"explicit":"implicit"));
1034:     PetscCall(PetscViewerASCIIPushTab(viewer));
1035:     PetscCall(EPSView(cyclic->eps,viewer));
1036:     PetscCall(PetscViewerASCIIPopTab(viewer));
1037:   }
1038:   PetscFunctionReturn(PETSC_SUCCESS);
1039: }

1041: static PetscErrorCode SVDReset_Cyclic(SVD svd)
1042: {
1043:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;

1045:   PetscFunctionBegin;
1046:   PetscCall(EPSReset(cyclic->eps));
1047:   PetscCall(MatDestroy(&cyclic->C));
1048:   PetscCall(MatDestroy(&cyclic->D));
1049:   PetscFunctionReturn(PETSC_SUCCESS);
1050: }

1052: static PetscErrorCode SVDDestroy_Cyclic(SVD svd)
1053: {
1054:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;

1056:   PetscFunctionBegin;
1057:   PetscCall(EPSDestroy(&cyclic->eps));
1058:   PetscCall(PetscFree(svd->data));
1059:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",NULL));
1060:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",NULL));
1061:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",NULL));
1062:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",NULL));
1063:   PetscFunctionReturn(PETSC_SUCCESS);
1064: }

1066: SLEPC_EXTERN PetscErrorCode SVDCreate_Cyclic(SVD svd)
1067: {
1068:   SVD_CYCLIC     *cyclic;

1070:   PetscFunctionBegin;
1071:   PetscCall(PetscNew(&cyclic));
1072:   svd->data                = (void*)cyclic;
1073:   svd->ops->solve          = SVDSolve_Cyclic;
1074:   svd->ops->solveg         = SVDSolve_Cyclic;
1075:   svd->ops->solveh         = SVDSolve_Cyclic;
1076:   svd->ops->setup          = SVDSetUp_Cyclic;
1077:   svd->ops->setfromoptions = SVDSetFromOptions_Cyclic;
1078:   svd->ops->destroy        = SVDDestroy_Cyclic;
1079:   svd->ops->reset          = SVDReset_Cyclic;
1080:   svd->ops->view           = SVDView_Cyclic;
1081:   svd->ops->computevectors = SVDComputeVectors_Cyclic;
1082:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",SVDCyclicSetEPS_Cyclic));
1083:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",SVDCyclicGetEPS_Cyclic));
1084:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",SVDCyclicSetExplicitMatrix_Cyclic));
1085:   PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",SVDCyclicGetExplicitMatrix_Cyclic));
1086:   PetscFunctionReturn(PETSC_SUCCESS);
1087: }