Actual source code: cyclic.c

slepc-3.22.1 2024-10-28
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) || defined(PETSC_HAVE_HIP)
 84:   PetscBool        gpu;
 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(MatGetOwnershipRange(Zm,&Istart,&Iend));
 99:     for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(Zm,i,i,0.0,INSERT_VALUES));
100:     PetscCall(MatAssemblyBegin(Zm,MAT_FINAL_ASSEMBLY));
101:     PetscCall(MatAssemblyEnd(Zm,MAT_FINAL_ASSEMBLY));
102:     PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Zn));
103:     PetscCall(MatSetSizes(Zn,n,n,N,N));
104:     PetscCall(MatSetFromOptions(Zn));
105:     PetscCall(MatGetOwnershipRange(Zn,&Istart,&Iend));
106:     for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(Zn,i,i,0.0,INSERT_VALUES));
107:     PetscCall(MatAssemblyBegin(Zn,MAT_FINAL_ASSEMBLY));
108:     PetscCall(MatAssemblyEnd(Zn,MAT_FINAL_ASSEMBLY));
109:     PetscCall(MatCreateTile(1.0,Zm,1.0,A,1.0,AT,1.0,Zn,C));
110:     PetscCall(MatDestroy(&Zm));
111:     PetscCall(MatDestroy(&Zn));
112:   } else {
113:     PetscCall(PetscNew(&ctx));
114:     ctx->A       = A;
115:     ctx->AT      = AT;
116:     ctx->swapped = svd->swapped;
117:     PetscCall(MatCreateVecsEmpty(A,&ctx->x2,&ctx->x1));
118:     PetscCall(MatCreateVecsEmpty(A,&ctx->y2,&ctx->y1));
119:     PetscCall(MatCreateShell(PetscObjectComm((PetscObject)svd),m+n,m+n,M+N,M+N,ctx,C));
120:     PetscCall(MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Cyclic));
121:     PetscCall(MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_Cyclic));
122: #if defined(PETSC_HAVE_CUDA)
123:     PetscCall(PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&gpu,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,""));
124:     if (gpu) PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic_CUDA));
125:     else
126: #elif defined(PETSC_HAVE_HIP)
127:     PetscCall(PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&gpu,MATSEQAIJHIPSPARSE,MATMPIAIJHIPSPARSE,""));
128:     if (gpu) PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic_HIP));
129:     else
130: #endif
131:       PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cyclic));
132:     PetscCall(MatGetVecType(A,&vtype));
133:     PetscCall(MatSetVecType(*C,vtype));
134: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
135:     if (gpu) {
136:       /* check alignment of bottom block */
137:       PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ctx->x1),&size));
138:       PetscCall(VecGetOwnershipRanges(ctx->x1,&ranges));
139:       for (i=0;i<size;i++) {
140:         ctx->misaligned = (((ranges[i+1]-ranges[i])*sizeof(PetscScalar))%16)? PETSC_TRUE: PETSC_FALSE;
141:         if (ctx->misaligned) break;
142:       }
143:       if (ctx->misaligned) {  /* create work vectors for MatMult */
144:         PetscCall(VecDuplicate(ctx->x2,&ctx->wx2));
145:         PetscCall(VecDuplicate(ctx->y2,&ctx->wy2));
146:       }
147:     }
148: #endif
149:   }
150:   PetscFunctionReturn(PETSC_SUCCESS);
151: }

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

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

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

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

237: static PetscErrorCode MatDestroy_ECross(Mat B)
238: {
239:   SVD_CYCLIC_SHELL *ctx;

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

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

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

281:   if (cyclic->explicitmatrix) {
282:     PetscCheck(svd->expltrans,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Cannot use explicit cyclic matrix with implicit transpose");
283:     PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)svd),m,m,M,M,1.0,&Id));
284:     PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Zm));
285:     PetscCall(MatSetSizes(Zm,m,n,M,N));
286:     PetscCall(MatSetFromOptions(Zm));
287:     PetscCall(MatGetOwnershipRange(Zm,&Istart,&Iend));
288:     for (i=Istart;i<Iend;i++) {
289:       if (i<N) PetscCall(MatSetValue(Zm,i,i,0.0,INSERT_VALUES));
290:     }
291:     PetscCall(MatAssemblyBegin(Zm,MAT_FINAL_ASSEMBLY));
292:     PetscCall(MatAssemblyEnd(Zm,MAT_FINAL_ASSEMBLY));
293:     PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Zn));
294:     PetscCall(MatSetSizes(Zn,n,m,N,M));
295:     PetscCall(MatSetFromOptions(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),&gpu,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,""));
327:     if (gpu) PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross_CUDA));
328:     else
329: #elif defined(PETSC_HAVE_HIP)
330:     PetscCall(PetscObjectTypeCompareAny((PetscObject)(svd->swapped?AT:A),&gpu,MATSEQAIJHIPSPARSE,MATMPIAIJHIPSPARSE,""));
331:     if (gpu) PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross_HIP));
332:     else
333: #endif
334:       PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_ECross));
335:     PetscCall(MatGetVecType(A,&vtype));
336:     PetscCall(MatSetVecType(*C,vtype));
337: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
338:     if (gpu) {
339:       /* check alignment of bottom block */
340:       PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ctx->x1),&size));
341:       PetscCall(VecGetOwnershipRanges(ctx->x1,&ranges));
342:       for (i=0;i<size;i++) {
343:         ctx->misaligned = (((ranges[i+1]-ranges[i])*sizeof(PetscScalar))%16)? PETSC_TRUE: PETSC_FALSE;
344:         if (ctx->misaligned) break;
345:       }
346:       if (ctx->misaligned) {  /* create work vectors for MatMult */
347:         PetscCall(VecDuplicate(ctx->x2,&ctx->wx2));
348:         PetscCall(VecDuplicate(ctx->y2,&ctx->wy2));
349:       }
350:     }
351: #endif
352:   }
353:   PetscFunctionReturn(PETSC_SUCCESS);
354: }

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

361:   PetscFunctionBegin;
362:   *errest = res/PetscMax(svd->nrma,svd->nrmb);
363:   PetscFunctionReturn(PETSC_SUCCESS);
364: }

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

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

517:   svd->leftbasis = PETSC_TRUE;
518:   PetscCall(SVDAllocateSolution(svd,0));
519:   PetscFunctionReturn(PETSC_SUCCESS);
520: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

858:   PetscFunctionBegin;
859:   PetscOptionsHeadBegin(PetscOptionsObject,"SVD Cyclic Options");

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

864:   PetscOptionsHeadEnd();

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

876: static PetscErrorCode SVDCyclicSetExplicitMatrix_Cyclic(SVD svd,PetscBool explicitmat)
877: {
878:   SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;

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

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

892:    Logically Collective

894:    Input Parameters:
895: +  svd         - singular value solver
896: -  explicitmat - boolean flag indicating if H(A) is built explicitly

898:    Options Database Key:
899: .  -svd_cyclic_explicitmatrix <boolean> - Indicates the boolean flag

901:    Level: advanced

903: .seealso: SVDCyclicGetExplicitMatrix()
904: @*/
905: PetscErrorCode SVDCyclicSetExplicitMatrix(SVD svd,PetscBool explicitmat)
906: {
907:   PetscFunctionBegin;
910:   PetscTryMethod(svd,"SVDCyclicSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
911:   PetscFunctionReturn(PETSC_SUCCESS);
912: }

914: static PetscErrorCode SVDCyclicGetExplicitMatrix_Cyclic(SVD svd,PetscBool *explicitmat)
915: {
916:   SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;

918:   PetscFunctionBegin;
919:   *explicitmat = cyclic->explicitmatrix;
920:   PetscFunctionReturn(PETSC_SUCCESS);
921: }

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

926:    Not Collective

928:    Input Parameter:
929: .  svd  - singular value solver

931:    Output Parameter:
932: .  explicitmat - the mode flag

934:    Level: advanced

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

947: static PetscErrorCode SVDCyclicSetEPS_Cyclic(SVD svd,EPS eps)
948: {
949:   SVD_CYCLIC      *cyclic = (SVD_CYCLIC*)svd->data;

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

960: /*@
961:    SVDCyclicSetEPS - Associate an eigensolver object (EPS) to the
962:    singular value solver.

964:    Collective

966:    Input Parameters:
967: +  svd - singular value solver
968: -  eps - the eigensolver object

970:    Level: advanced

972: .seealso: SVDCyclicGetEPS()
973: @*/
974: PetscErrorCode SVDCyclicSetEPS(SVD svd,EPS eps)
975: {
976:   PetscFunctionBegin;
979:   PetscCheckSameComm(svd,1,eps,2);
980:   PetscTryMethod(svd,"SVDCyclicSetEPS_C",(SVD,EPS),(svd,eps));
981:   PetscFunctionReturn(PETSC_SUCCESS);
982: }

984: static PetscErrorCode SVDCyclicGetEPS_Cyclic(SVD svd,EPS *eps)
985: {
986:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;

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

1002: /*@
1003:    SVDCyclicGetEPS - Retrieve the eigensolver object (EPS) associated
1004:    to the singular value solver.

1006:    Collective

1008:    Input Parameter:
1009: .  svd - singular value solver

1011:    Output Parameter:
1012: .  eps - the eigensolver object

1014:    Level: advanced

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

1027: static PetscErrorCode SVDView_Cyclic(SVD svd,PetscViewer viewer)
1028: {
1029:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;
1030:   PetscBool      isascii;

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

1044: static PetscErrorCode SVDReset_Cyclic(SVD svd)
1045: {
1046:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;

1048:   PetscFunctionBegin;
1049:   PetscCall(EPSReset(cyclic->eps));
1050:   PetscCall(MatDestroy(&cyclic->C));
1051:   PetscCall(MatDestroy(&cyclic->D));
1052:   PetscFunctionReturn(PETSC_SUCCESS);
1053: }

1055: static PetscErrorCode SVDDestroy_Cyclic(SVD svd)
1056: {
1057:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;

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

1069: SLEPC_EXTERN PetscErrorCode SVDCreate_Cyclic(SVD svd)
1070: {
1071:   SVD_CYCLIC     *cyclic;

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