Actual source code: cyclic.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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>                /*I "slepcsvd.h" I*/
 17: #include <slepc/private/epsimpl.h>                /*I "slepceps.h" I*/
 18: #include "./cyclicimpl.h"

 20: static PetscErrorCode MatMult_Cyclic(Mat B,Vec x,Vec y)
 21: {
 22:   PetscErrorCode    ierr;
 23:   SVD               svd;
 24:   SVD_CYCLIC        *cyclic;
 25:   const PetscScalar *px;
 26:   PetscScalar       *py;
 27:   PetscInt          m;

 30:   MatShellGetContext(B,(void**)&svd);
 31:   cyclic = (SVD_CYCLIC*)svd->data;
 32:   SVDMatGetLocalSize(svd,&m,NULL);
 33:   VecGetArrayRead(x,&px);
 34:   VecGetArray(y,&py);
 35:   VecPlaceArray(cyclic->x1,px);
 36:   VecPlaceArray(cyclic->x2,px+m);
 37:   VecPlaceArray(cyclic->y1,py);
 38:   VecPlaceArray(cyclic->y2,py+m);
 39:   SVDMatMult(svd,PETSC_FALSE,cyclic->x2,cyclic->y1);
 40:   SVDMatMult(svd,PETSC_TRUE,cyclic->x1,cyclic->y2);
 41:   VecResetArray(cyclic->x1);
 42:   VecResetArray(cyclic->x2);
 43:   VecResetArray(cyclic->y1);
 44:   VecResetArray(cyclic->y2);
 45:   VecRestoreArrayRead(x,&px);
 46:   VecRestoreArray(y,&py);
 47:   return(0);
 48: }

 50: static PetscErrorCode MatGetDiagonal_Cyclic(Mat B,Vec diag)
 51: {

 55:   VecSet(diag,0.0);
 56:   return(0);
 57: }

 59: PetscErrorCode SVDSetUp_Cyclic(SVD svd)
 60: {
 61:   PetscErrorCode    ierr;
 62:   SVD_CYCLIC        *cyclic = (SVD_CYCLIC*)svd->data;
 63:   PetscInt          M,N,m,n,i,isl,Istart,Iend;
 64:   const PetscScalar *isa;
 65:   PetscScalar       *va;
 66:   PetscBool         trackall,issinv;
 67:   Vec               v;
 68:   Mat               Zm,Zn;
 69:   ST                st;
 70: #if defined(PETSC_HAVE_CUDA)
 71:   PetscBool         cuda;
 72: #endif

 75:   SVDMatGetSize(svd,&M,&N);
 76:   SVDMatGetLocalSize(svd,&m,&n);
 77:   if (!cyclic->mat) {
 78:     if (cyclic->explicitmatrix) {
 79:       if (!svd->A || !svd->AT) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Cannot use explicit cyclic matrix with implicit transpose");
 80:       MatCreate(PetscObjectComm((PetscObject)svd),&Zm);
 81:       MatSetSizes(Zm,m,m,M,M);
 82:       MatSetFromOptions(Zm);
 83:       MatSetUp(Zm);
 84:       MatGetOwnershipRange(Zm,&Istart,&Iend);
 85:       for (i=Istart;i<Iend;i++) {
 86:         MatSetValue(Zm,i,i,0.0,INSERT_VALUES);
 87:       }
 88:       MatAssemblyBegin(Zm,MAT_FINAL_ASSEMBLY);
 89:       MatAssemblyEnd(Zm,MAT_FINAL_ASSEMBLY);
 90:       MatCreate(PetscObjectComm((PetscObject)svd),&Zn);
 91:       MatSetSizes(Zn,n,n,N,N);
 92:       MatSetFromOptions(Zn);
 93:       MatSetUp(Zn);
 94:       MatGetOwnershipRange(Zn,&Istart,&Iend);
 95:       for (i=Istart;i<Iend;i++) {
 96:         MatSetValue(Zn,i,i,0.0,INSERT_VALUES);
 97:       }
 98:       MatAssemblyBegin(Zn,MAT_FINAL_ASSEMBLY);
 99:       MatAssemblyEnd(Zn,MAT_FINAL_ASSEMBLY);
100:       MatCreateTile(1.0,Zm,1.0,svd->A,1.0,svd->AT,1.0,Zn,&cyclic->mat);
101:       PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->mat);
102:       MatDestroy(&Zm);
103:       MatDestroy(&Zn);
104:     } else {
105:       SVDMatCreateVecsEmpty(svd,&cyclic->x2,&cyclic->x1);
106:       SVDMatCreateVecsEmpty(svd,&cyclic->y2,&cyclic->y1);
107:       PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->x1);
108:       PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->x2);
109:       PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->y1);
110:       PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->y2);
111:       MatCreateShell(PetscObjectComm((PetscObject)svd),m+n,m+n,M+N,M+N,svd,&cyclic->mat);
112:       MatShellSetOperation(cyclic->mat,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Cyclic);
113: #if defined(PETSC_HAVE_CUDA)
114:       PetscObjectTypeCompareAny((PetscObject)(svd->A?svd->A:svd->AT),&cuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,"");
115:       if (cuda) {
116:         MatShellSetOperation(cyclic->mat,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_Cyclic_CUDA);
117:         MatShellSetOperation(cyclic->mat,MATOP_MULT,(void(*)(void))MatMult_Cyclic_CUDA);
118:       } else
119: #endif
120:       {
121:         MatShellSetOperation(cyclic->mat,MATOP_MULT,(void(*)(void))MatMult_Cyclic);
122:       }
123:     }
124:     PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->mat);
125:   }

127:   if (!cyclic->eps) { SVDCyclicGetEPS(svd,&cyclic->eps); }
128:   EPSSetOperators(cyclic->eps,cyclic->mat,NULL);
129:   EPSSetProblemType(cyclic->eps,EPS_HEP);
130:   if (!cyclic->usereps) {
131:     if (svd->which == SVD_LARGEST) {
132:       EPSGetST(cyclic->eps,&st);
133:       PetscObjectTypeCompare((PetscObject)st,STSINVERT,&issinv);
134:       if (issinv) {
135:         EPSSetWhichEigenpairs(cyclic->eps,EPS_TARGET_MAGNITUDE);
136:       } else {
137:         EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL);
138:       }
139:     } else {
140:       EPSSetEigenvalueComparison(cyclic->eps,SlepcCompareSmallestPosReal,NULL);
141:       EPSSetTarget(cyclic->eps,0.0);
142:     }
143:     EPSSetDimensions(cyclic->eps,svd->nsv,svd->ncv?svd->ncv:PETSC_DEFAULT,svd->mpd?svd->mpd:PETSC_DEFAULT);
144:     EPSSetTolerances(cyclic->eps,svd->tol==PETSC_DEFAULT?SLEPC_DEFAULT_TOL/10.0:svd->tol,svd->max_it?svd->max_it:PETSC_DEFAULT);
145:     switch (svd->conv) {
146:     case SVD_CONV_ABS:
147:       EPSSetConvergenceTest(cyclic->eps,EPS_CONV_ABS);break;
148:     case SVD_CONV_REL:
149:       EPSSetConvergenceTest(cyclic->eps,EPS_CONV_REL);break;
150:     case SVD_CONV_USER:
151:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
152:     }
153:   }
154:   if (svd->stop!=SVD_STOP_BASIC) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined stopping test not supported in this solver");
155:   /* Transfer the trackall option from svd to eps */
156:   SVDGetTrackAll(svd,&trackall);
157:   EPSSetTrackAll(cyclic->eps,trackall);
158:   /* Transfer the initial subspace from svd to eps */
159:   if (svd->nini<0 || svd->ninil<0) {
160:     for (i=0;i<-PetscMin(svd->nini,svd->ninil);i++) {
161:       MatCreateVecs(cyclic->mat,&v,NULL);
162:       VecGetArray(v,&va);
163:       if (i<-svd->ninil) {
164:         VecGetSize(svd->ISL[i],&isl);
165:         if (isl!=m) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Size mismatch for left initial vector");
166:         VecGetArrayRead(svd->ISL[i],&isa);
167:         PetscMemcpy(va,isa,sizeof(PetscScalar)*m);
168:         VecRestoreArrayRead(svd->IS[i],&isa);
169:       } else {
170:         PetscMemzero(&va,sizeof(PetscScalar)*m);
171:       }
172:       if (i<-svd->nini) {
173:         VecGetSize(svd->IS[i],&isl);
174:         if (isl!=n) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Size mismatch for right initial vector");
175:         VecGetArrayRead(svd->IS[i],&isa);
176:         PetscMemcpy(va+m,isa,sizeof(PetscScalar)*n);
177:         VecRestoreArrayRead(svd->IS[i],&isa);
178:       } else {
179:         PetscMemzero(va+m,sizeof(PetscScalar)*n);
180:       }
181:       VecRestoreArray(v,&va);
182:       VecDestroy(&svd->IS[i]);
183:       svd->IS[i] = v;
184:     }
185:     svd->nini = PetscMin(svd->nini,svd->ninil);
186:     EPSSetInitialSpace(cyclic->eps,-svd->nini,svd->IS);
187:     SlepcBasisDestroy_Private(&svd->nini,&svd->IS);
188:     SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL);
189:   }
190:   EPSSetUp(cyclic->eps);
191:   EPSGetDimensions(cyclic->eps,NULL,&svd->ncv,&svd->mpd);
192:   svd->ncv = PetscMin(svd->ncv,PetscMin(M,N));
193:   EPSGetTolerances(cyclic->eps,NULL,&svd->max_it);
194:   if (svd->tol==PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;

196:   svd->leftbasis = PETSC_TRUE;
197:   SVDAllocateSolution(svd,0);
198:   return(0);
199: }

201: PetscErrorCode SVDSolve_Cyclic(SVD svd)
202: {
203:   PetscErrorCode    ierr;
204:   SVD_CYCLIC        *cyclic = (SVD_CYCLIC*)svd->data;
205:   PetscInt          i,j,M,N,m,n;
206:   PetscScalar       sigma;
207:   const PetscScalar *px;
208:   Vec               x,x1,x2;

211:   EPSSolve(cyclic->eps);
212:   EPSGetConverged(cyclic->eps,&svd->nconv);
213:   EPSGetIterationNumber(cyclic->eps,&svd->its);
214:   EPSGetConvergedReason(cyclic->eps,(EPSConvergedReason*)&svd->reason);

216:   MatCreateVecs(cyclic->mat,&x,NULL);
217:   SVDMatGetSize(svd,&M,&N);
218:   SVDMatGetLocalSize(svd,&m,&n);
219:   SVDMatCreateVecsEmpty(svd,&x2,&x1);
220:   for (i=0,j=0;i<svd->nconv;i++) {
221:     EPSGetEigenpair(cyclic->eps,i,&sigma,NULL,x,NULL);
222:     if (PetscRealPart(sigma) > 0.0) {
223:       svd->sigma[j] = PetscRealPart(sigma);
224:       VecGetArrayRead(x,&px);
225:       VecPlaceArray(x1,px);
226:       VecPlaceArray(x2,px+m);
227:       BVInsertVec(svd->U,j,x1);
228:       BVScaleColumn(svd->U,j,1.0/PetscSqrtReal(2.0));
229:       BVInsertVec(svd->V,j,x2);
230:       BVScaleColumn(svd->V,j,1.0/PetscSqrtReal(2.0));
231:       VecResetArray(x1);
232:       VecResetArray(x2);
233:       VecRestoreArrayRead(x,&px);
234:       j++;
235:     }
236:   }
237:   svd->nconv = j;

239:   VecDestroy(&x);
240:   VecDestroy(&x1);
241:   VecDestroy(&x2);
242:   return(0);
243: }

245: static PetscErrorCode EPSMonitor_Cyclic(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
246: {
247:   PetscInt       i,j;
248:   SVD            svd = (SVD)ctx;
249:   PetscScalar    er,ei;

253:   nconv = 0;
254:   for (i=0,j=0;i<PetscMin(nest,svd->ncv);i++) {
255:     er = eigr[i]; ei = eigi[i];
256:     STBackTransform(eps->st,1,&er,&ei);
257:     if (PetscRealPart(er) > 0.0) {
258:       svd->sigma[j] = PetscRealPart(er);
259:       svd->errest[j] = errest[i];
260:       if (errest[i] && errest[i] < svd->tol) nconv++;
261:       j++;
262:     }
263:   }
264:   nest = j;
265:   SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest);
266:   return(0);
267: }

269: PetscErrorCode SVDSetFromOptions_Cyclic(PetscOptionItems *PetscOptionsObject,SVD svd)
270: {
272:   PetscBool      set,val;
273:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;
274:   ST             st;

277:   PetscOptionsHead(PetscOptionsObject,"SVD Cyclic Options");

279:     PetscOptionsBool("-svd_cyclic_explicitmatrix","Use cyclic explicit matrix","SVDCyclicSetExplicitMatrix",cyclic->explicitmatrix,&val,&set);
280:     if (set) { SVDCyclicSetExplicitMatrix(svd,val); }

282:   PetscOptionsTail();

284:   if (!cyclic->eps) { SVDCyclicGetEPS(svd,&cyclic->eps); }
285:   if (!cyclic->explicitmatrix && !cyclic->usereps) {
286:     /* use as default an ST with shell matrix and Jacobi */
287:     EPSGetST(cyclic->eps,&st);
288:     STSetMatMode(st,ST_MATMODE_SHELL);
289:   }
290:   EPSSetFromOptions(cyclic->eps);
291:   return(0);
292: }

294: static PetscErrorCode SVDCyclicSetExplicitMatrix_Cyclic(SVD svd,PetscBool explicitmatrix)
295: {
296:   SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;

299:   if (cyclic->explicitmatrix != explicitmatrix) {
300:     cyclic->explicitmatrix = explicitmatrix;
301:     svd->state = SVD_STATE_INITIAL;
302:   }
303:   return(0);
304: }

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

310:    Logically Collective on SVD

312:    Input Parameters:
313: +  svd      - singular value solver
314: -  explicit - boolean flag indicating if H(A) is built explicitly

316:    Options Database Key:
317: .  -svd_cyclic_explicitmatrix <boolean> - Indicates the boolean flag

319:    Level: advanced

321: .seealso: SVDCyclicGetExplicitMatrix()
322: @*/
323: PetscErrorCode SVDCyclicSetExplicitMatrix(SVD svd,PetscBool explicitmatrix)
324: {

330:   PetscTryMethod(svd,"SVDCyclicSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmatrix));
331:   return(0);
332: }

334: static PetscErrorCode SVDCyclicGetExplicitMatrix_Cyclic(SVD svd,PetscBool *explicitmatrix)
335: {
336:   SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;

339:   *explicitmatrix = cyclic->explicitmatrix;
340:   return(0);
341: }

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

346:    Not Collective

348:    Input Parameter:
349: .  svd  - singular value solver

351:    Output Parameter:
352: .  explicit - the mode flag

354:    Level: advanced

356: .seealso: SVDCyclicSetExplicitMatrix()
357: @*/
358: PetscErrorCode SVDCyclicGetExplicitMatrix(SVD svd,PetscBool *explicitmatrix)
359: {

365:   PetscUseMethod(svd,"SVDCyclicGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmatrix));
366:   return(0);
367: }

369: static PetscErrorCode SVDCyclicSetEPS_Cyclic(SVD svd,EPS eps)
370: {
371:   PetscErrorCode  ierr;
372:   SVD_CYCLIC      *cyclic = (SVD_CYCLIC*)svd->data;

375:   PetscObjectReference((PetscObject)eps);
376:   EPSDestroy(&cyclic->eps);
377:   cyclic->eps = eps;
378:   cyclic->usereps = PETSC_TRUE;
379:   PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->eps);
380:   svd->state = SVD_STATE_INITIAL;
381:   return(0);
382: }

384: /*@
385:    SVDCyclicSetEPS - Associate an eigensolver object (EPS) to the
386:    singular value solver.

388:    Collective on SVD

390:    Input Parameters:
391: +  svd - singular value solver
392: -  eps - the eigensolver object

394:    Level: advanced

396: .seealso: SVDCyclicGetEPS()
397: @*/
398: PetscErrorCode SVDCyclicSetEPS(SVD svd,EPS eps)
399: {

406:   PetscTryMethod(svd,"SVDCyclicSetEPS_C",(SVD,EPS),(svd,eps));
407:   return(0);
408: }

410: static PetscErrorCode SVDCyclicGetEPS_Cyclic(SVD svd,EPS *eps)
411: {
413:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;

416:   if (!cyclic->eps) {
417:     EPSCreate(PetscObjectComm((PetscObject)svd),&cyclic->eps);
418:     PetscObjectIncrementTabLevel((PetscObject)cyclic->eps,(PetscObject)svd,1);
419:     EPSSetOptionsPrefix(cyclic->eps,((PetscObject)svd)->prefix);
420:     EPSAppendOptionsPrefix(cyclic->eps,"svd_cyclic_");
421:     PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->eps);
422:     PetscObjectSetOptions((PetscObject)cyclic->eps,((PetscObject)svd)->options);
423:     EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL);
424:     EPSMonitorSet(cyclic->eps,EPSMonitor_Cyclic,svd,NULL);
425:   }
426:   *eps = cyclic->eps;
427:   return(0);
428: }

430: /*@
431:    SVDCyclicGetEPS - Retrieve the eigensolver object (EPS) associated
432:    to the singular value solver.

434:    Not Collective

436:    Input Parameter:
437: .  svd - singular value solver

439:    Output Parameter:
440: .  eps - the eigensolver object

442:    Level: advanced

444: .seealso: SVDCyclicSetEPS()
445: @*/
446: PetscErrorCode SVDCyclicGetEPS(SVD svd,EPS *eps)
447: {

453:   PetscUseMethod(svd,"SVDCyclicGetEPS_C",(SVD,EPS*),(svd,eps));
454:   return(0);
455: }

457: PetscErrorCode SVDView_Cyclic(SVD svd,PetscViewer viewer)
458: {
460:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;
461:   PetscBool      isascii;

464:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
465:   if (isascii) {
466:     if (!cyclic->eps) { SVDCyclicGetEPS(svd,&cyclic->eps); }
467:     PetscViewerASCIIPrintf(viewer,"  %s matrix\n",cyclic->explicitmatrix?"explicit":"implicit");
468:     PetscViewerASCIIPushTab(viewer);
469:     EPSView(cyclic->eps,viewer);
470:     PetscViewerASCIIPopTab(viewer);
471:   }
472:   return(0);
473: }

475: PetscErrorCode SVDReset_Cyclic(SVD svd)
476: {
478:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;

481:   EPSReset(cyclic->eps);
482:   MatDestroy(&cyclic->mat);
483:   VecDestroy(&cyclic->x1);
484:   VecDestroy(&cyclic->x2);
485:   VecDestroy(&cyclic->y1);
486:   VecDestroy(&cyclic->y2);
487:   return(0);
488: }

490: PetscErrorCode SVDDestroy_Cyclic(SVD svd)
491: {
493:   SVD_CYCLIC     *cyclic = (SVD_CYCLIC*)svd->data;

496:   EPSDestroy(&cyclic->eps);
497:   PetscFree(svd->data);
498:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",NULL);
499:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",NULL);
500:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",NULL);
501:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",NULL);
502:   return(0);
503: }

505: SLEPC_EXTERN PetscErrorCode SVDCreate_Cyclic(SVD svd)
506: {
508:   SVD_CYCLIC     *cyclic;

511:   PetscNewLog(svd,&cyclic);
512:   svd->data                      = (void*)cyclic;
513:   svd->ops->solve                = SVDSolve_Cyclic;
514:   svd->ops->setup                = SVDSetUp_Cyclic;
515:   svd->ops->setfromoptions       = SVDSetFromOptions_Cyclic;
516:   svd->ops->destroy              = SVDDestroy_Cyclic;
517:   svd->ops->reset                = SVDReset_Cyclic;
518:   svd->ops->view                 = SVDView_Cyclic;
519:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",SVDCyclicSetEPS_Cyclic);
520:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",SVDCyclicGetEPS_Cyclic);
521:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",SVDCyclicSetExplicitMatrix_Cyclic);
522:   PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",SVDCyclicGetExplicitMatrix_Cyclic);
523:   return(0);
524: }