Actual source code: blopex.c

slepc-3.16.2 2022-02-01
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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:    This file implements a wrapper to the BLOPEX package
 12: */

 14: #include <slepc/private/epsimpl.h>
 15: #include "blopex.h"
 16: #include <lobpcg.h>
 17: #include <interpreter.h>
 18: #include <multivector.h>
 19: #include <temp_multivector.h>

 21: PetscInt slepc_blopex_useconstr = -1;

 23: typedef struct {
 24:   lobpcg_Tolerance           tol;
 25:   lobpcg_BLASLAPACKFunctions blap_fn;
 26:   mv_InterfaceInterpreter    ii;
 27:   ST                         st;
 28:   Vec                        w;
 29:   PetscInt                   bs;     /* block size */
 30: } EPS_BLOPEX;

 32: static void Precond_FnSingleVector(void *data,void *x,void *y)
 33: {
 35:   EPS_BLOPEX     *blopex = (EPS_BLOPEX*)data;
 36:   MPI_Comm       comm = PetscObjectComm((PetscObject)blopex->st);
 37:   KSP            ksp;

 40:   STGetKSP(blopex->st,&ksp);CHKERRABORT(comm,ierr);
 41:   KSPSolve(ksp,(Vec)x,(Vec)y);CHKERRABORT(comm,ierr);
 42:   PetscFunctionReturnVoid();
 43: }

 45: static void Precond_FnMultiVector(void *data,void *x,void *y)
 46: {
 47:   EPS_BLOPEX *blopex = (EPS_BLOPEX*)data;

 50:   blopex->ii.Eval(Precond_FnSingleVector,data,x,y);
 51:   PetscFunctionReturnVoid();
 52: }

 54: static void OperatorASingleVector(void *data,void *x,void *y)
 55: {
 57:   EPS_BLOPEX     *blopex = (EPS_BLOPEX*)data;
 58:   MPI_Comm       comm = PetscObjectComm((PetscObject)blopex->st);
 59:   Mat            A,B;
 60:   PetscScalar    sigma;
 61:   PetscInt       nmat;

 64:   STGetNumMatrices(blopex->st,&nmat);CHKERRABORT(comm,ierr);
 65:   STGetMatrix(blopex->st,0,&A);CHKERRABORT(comm,ierr);
 66:   if (nmat>1) { STGetMatrix(blopex->st,1,&B);CHKERRABORT(comm,ierr); }
 67:   MatMult(A,(Vec)x,(Vec)y);CHKERRABORT(comm,ierr);
 68:   STGetShift(blopex->st,&sigma);CHKERRABORT(comm,ierr);
 69:   if (sigma != 0.0) {
 70:     if (nmat>1) {
 71:       MatMult(B,(Vec)x,blopex->w);CHKERRABORT(comm,ierr);
 72:     } else {
 73:       VecCopy((Vec)x,blopex->w);CHKERRABORT(comm,ierr);
 74:     }
 75:     VecAXPY((Vec)y,-sigma,blopex->w);CHKERRABORT(comm,ierr);
 76:   }
 77:   PetscFunctionReturnVoid();
 78: }

 80: static void OperatorAMultiVector(void *data,void *x,void *y)
 81: {
 82:   EPS_BLOPEX *blopex = (EPS_BLOPEX*)data;

 85:   blopex->ii.Eval(OperatorASingleVector,data,x,y);
 86:   PetscFunctionReturnVoid();
 87: }

 89: static void OperatorBSingleVector(void *data,void *x,void *y)
 90: {
 92:   EPS_BLOPEX     *blopex = (EPS_BLOPEX*)data;
 93:   MPI_Comm       comm = PetscObjectComm((PetscObject)blopex->st);
 94:   Mat            B;

 97:   STGetMatrix(blopex->st,1,&B);CHKERRABORT(comm,ierr);
 98:   MatMult(B,(Vec)x,(Vec)y);CHKERRABORT(comm,ierr);
 99:   PetscFunctionReturnVoid();
100: }

102: static void OperatorBMultiVector(void *data,void *x,void *y)
103: {
104:   EPS_BLOPEX *blopex = (EPS_BLOPEX*)data;

107:   blopex->ii.Eval(OperatorBSingleVector,data,x,y);
108:   PetscFunctionReturnVoid();
109: }

111: PetscErrorCode EPSSetDimensions_BLOPEX(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
112: {
114:   EPS_BLOPEX     *ctx = (EPS_BLOPEX*)eps->data;
115:   PetscInt       k;

118:   k = ((eps->nev-1)/ctx->bs+1)*ctx->bs;
119:   if (*ncv!=PETSC_DEFAULT) { /* ncv set */
120:     if (*ncv<k) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv is not sufficiently large");
121:   } else *ncv = k;
122:   if (*mpd==PETSC_DEFAULT) *mpd = *ncv;
123:   else { PetscInfo(eps,"Warning: given value of mpd ignored\n"); }
124:   return(0);
125: }

127: PetscErrorCode EPSSetUp_BLOPEX(EPS eps)
128: {
130:   EPS_BLOPEX     *blopex = (EPS_BLOPEX*)eps->data;
131:   PetscBool      flg;
132:   KSP            ksp;

135:   EPSCheckHermitianDefinite(eps);
136:   if (!blopex->bs) blopex->bs = PetscMin(16,eps->nev);
137:   EPSSetDimensions_BLOPEX(eps,eps->nev,&eps->ncv,&eps->mpd);
138:   if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
139:   if (!eps->which) eps->which = EPS_SMALLEST_REAL;
140:   if (eps->which!=EPS_SMALLEST_REAL) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only smallest real eigenvalues");
141:   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_STOPPING);
142:   EPSCheckIgnored(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_EXTRACTION);

144:   blopex->st = eps->st;

146:   if (eps->converged == EPSConvergedRelative) {
147:     blopex->tol.absolute = 0.0;
148:     blopex->tol.relative = SlepcDefaultTol(eps->tol);
149:   } else if (eps->converged == EPSConvergedAbsolute) {
150:     blopex->tol.absolute = SlepcDefaultTol(eps->tol);
151:     blopex->tol.relative = 0.0;
152:   } else SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Convergence test not supported in this solver");

154:   SLEPCSetupInterpreter(&blopex->ii);

156:   STGetKSP(eps->st,&ksp);
157:   PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);
158:   if (!flg) { PetscInfo(eps,"Warning: ignoring KSP, should use KSPPREONLY\n"); }

160:   /* allocate memory */
161:   if (!eps->V) { EPSGetBV(eps,&eps->V); }
162:   PetscObjectTypeCompareAny((PetscObject)eps->V,&flg,BVVECS,BVCONTIGUOUS,"");
163:   if (!flg) {  /* blopex only works with BVVECS or BVCONTIGUOUS */
164:     BVSetType(eps->V,BVCONTIGUOUS);
165:   }
166:   EPSAllocateSolution(eps,0);
167:   if (!blopex->w) {
168:     BVCreateVec(eps->V,&blopex->w);
169:     PetscLogObjectParent((PetscObject)eps,(PetscObject)blopex->w);
170:   }

172: #if defined(PETSC_USE_COMPLEX)
173:   blopex->blap_fn.zpotrf = PETSC_zpotrf_interface;
174:   blopex->blap_fn.zhegv = PETSC_zsygv_interface;
175: #else
176:   blopex->blap_fn.dpotrf = PETSC_dpotrf_interface;
177:   blopex->blap_fn.dsygv = PETSC_dsygv_interface;
178: #endif
179:   return(0);
180: }

182: PetscErrorCode EPSSolve_BLOPEX(EPS eps)
183: {
184:   EPS_BLOPEX        *blopex = (EPS_BLOPEX*)eps->data;
185:   PetscScalar       sigma,*eigr=NULL;
186:   PetscReal         *errest=NULL;
187:   int               i,j,info,its,nconv;
188:   double            *residhist=NULL;
189:   PetscErrorCode    ierr;
190:   mv_MultiVectorPtr eigenvectors,constraints;
191: #if defined(PETSC_USE_COMPLEX)
192:   komplex           *lambda=NULL,*lambdahist=NULL;
193: #else
194:   double            *lambda=NULL,*lambdahist=NULL;
195: #endif

198:   STGetShift(eps->st,&sigma);
199:   PetscMalloc1(blopex->bs,&lambda);
200:   if (eps->numbermonitors>0) {
201:     PetscMalloc4(blopex->bs*(eps->max_it+1),&lambdahist,eps->ncv,&eigr,blopex->bs*(eps->max_it+1),&residhist,eps->ncv,&errest);
202:   }

204:   /* Complete the initial basis with random vectors */
205:   for (i=0;i<eps->nini;i++) {  /* in case the initial vectors were also set with VecSetRandom */
206:     BVSetRandomColumn(eps->V,eps->nini);
207:   }
208:   for (i=eps->nini;i<eps->ncv;i++) {
209:     BVSetRandomColumn(eps->V,i);
210:   }

212:   while (eps->reason == EPS_CONVERGED_ITERATING) {

214:     /* Create multivector of constraints from leading columns of V */
215:     PetscObjectComposedDataSetInt((PetscObject)eps->V,slepc_blopex_useconstr,1);
216:     BVSetActiveColumns(eps->V,0,eps->nconv);
217:     constraints = mv_MultiVectorCreateFromSampleVector(&blopex->ii,eps->nds+eps->nconv,eps->V);

219:     /* Create multivector where eigenvectors of this run will be stored */
220:     PetscObjectComposedDataSetInt((PetscObject)eps->V,slepc_blopex_useconstr,0);
221:     BVSetActiveColumns(eps->V,eps->nconv,eps->nconv+blopex->bs);
222:     eigenvectors = mv_MultiVectorCreateFromSampleVector(&blopex->ii,blopex->bs,eps->V);

224: #if defined(PETSC_USE_COMPLEX)
225:     info = lobpcg_solve_complex(eigenvectors,blopex,OperatorAMultiVector,
226:           eps->isgeneralized?blopex:NULL,eps->isgeneralized?OperatorBMultiVector:NULL,
227:           blopex,Precond_FnMultiVector,constraints,
228:           blopex->blap_fn,blopex->tol,eps->max_it,0,&its,
229:           lambda,lambdahist,blopex->bs,eps->errest+eps->nconv,residhist,blopex->bs);
230: #else
231:     info = lobpcg_solve_double(eigenvectors,blopex,OperatorAMultiVector,
232:           eps->isgeneralized?blopex:NULL,eps->isgeneralized?OperatorBMultiVector:NULL,
233:           blopex,Precond_FnMultiVector,constraints,
234:           blopex->blap_fn,blopex->tol,eps->max_it,0,&its,
235:           lambda,lambdahist,blopex->bs,eps->errest+eps->nconv,residhist,blopex->bs);
236: #endif
237:     if (info>0) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"BLOPEX failed with exit code=%d",info);
238:     mv_MultiVectorDestroy(constraints);
239:     mv_MultiVectorDestroy(eigenvectors);

241:     for (j=0;j<blopex->bs;j++) {
242: #if defined(PETSC_USE_COMPLEX)
243:       eps->eigr[eps->nconv+j] = PetscCMPLX(lambda[j].real,lambda[j].imag);
244: #else
245:       eps->eigr[eps->nconv+j] = lambda[j];
246: #endif
247:     }

249:     if (eps->numbermonitors>0) {
250:       for (i=0;i<its;i++) {
251:         nconv = 0;
252:         for (j=0;j<blopex->bs;j++) {
253: #if defined(PETSC_USE_COMPLEX)
254:           eigr[eps->nconv+j] = PetscCMPLX(lambdahist[j+i*blopex->bs].real,lambdahist[j+i*blopex->bs].imag);
255: #else
256:           eigr[eps->nconv+j] = lambdahist[j+i*blopex->bs];
257: #endif
258:           errest[eps->nconv+j] = residhist[j+i*blopex->bs];
259:           if (residhist[j+i*blopex->bs]<=eps->tol) nconv++;
260:         }
261:         EPSMonitor(eps,eps->its+i,eps->nconv+nconv,eigr,eps->eigi,errest,eps->nconv+blopex->bs);
262:       }
263:     }

265:     eps->its += its;
266:     if (info==-1) {
267:       eps->reason = EPS_DIVERGED_ITS;
268:       break;
269:     } else {
270:       for (i=0;i<blopex->bs;i++) {
271:         if (sigma != 0.0) eps->eigr[eps->nconv+i] += sigma;
272:       }
273:       eps->nconv += blopex->bs;
274:       if (eps->nconv>=eps->nev) eps->reason = EPS_CONVERGED_TOL;
275:     }
276:   }

278:   PetscFree(lambda);
279:   if (eps->numbermonitors>0) {
280:     PetscFree4(lambdahist,eigr,residhist,errest);
281:   }
282:   return(0);
283: }

285: static PetscErrorCode EPSBLOPEXSetBlockSize_BLOPEX(EPS eps,PetscInt bs)
286: {
287:   EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;

290:   if (bs==PETSC_DEFAULT) {
291:     ctx->bs    = 0;
292:     eps->state = EPS_STATE_INITIAL;
293:   } else {
294:     if (bs<=0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block size must be >0");
295:     ctx->bs = bs;
296:   }
297:   return(0);
298: }

300: /*@
301:    EPSBLOPEXSetBlockSize - Sets the block size of the BLOPEX solver.

303:    Logically Collective on eps

305:    Input Parameters:
306: +  eps - the eigenproblem solver context
307: -  bs  - the block size

309:    Options Database Key:
310: .  -eps_blopex_blocksize - Sets the block size

312:    Level: advanced

314: .seealso: EPSBLOPEXGetBlockSize()
315: @*/
316: PetscErrorCode EPSBLOPEXSetBlockSize(EPS eps,PetscInt bs)
317: {

323:   PetscTryMethod(eps,"EPSBLOPEXSetBlockSize_C",(EPS,PetscInt),(eps,bs));
324:   return(0);
325: }

327: static PetscErrorCode EPSBLOPEXGetBlockSize_BLOPEX(EPS eps,PetscInt *bs)
328: {
329:   EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;

332:   *bs = ctx->bs;
333:   return(0);
334: }

336: /*@
337:    EPSBLOPEXGetBlockSize - Gets the block size used in the BLOPEX solver.

339:    Not Collective

341:    Input Parameter:
342: .  eps - the eigenproblem solver context

344:    Output Parameter:
345: .  bs - the block size

347:    Level: advanced

349: .seealso: EPSBLOPEXSetBlockSize()
350: @*/
351: PetscErrorCode EPSBLOPEXGetBlockSize(EPS eps,PetscInt *bs)
352: {

358:   PetscUseMethod(eps,"EPSBLOPEXGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
359:   return(0);
360: }

362: PetscErrorCode EPSReset_BLOPEX(EPS eps)
363: {
365:   EPS_BLOPEX     *blopex = (EPS_BLOPEX*)eps->data;

368:   VecDestroy(&blopex->w);
369:   return(0);
370: }

372: PetscErrorCode EPSDestroy_BLOPEX(EPS eps)
373: {

377:   LOBPCG_DestroyRandomContext();
378:   PetscFree(eps->data);
379:   PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXSetBlockSize_C",NULL);
380:   PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXGetBlockSize_C",NULL);
381:   return(0);
382: }

384: PetscErrorCode EPSView_BLOPEX(EPS eps,PetscViewer viewer)
385: {
387:   EPS_BLOPEX     *ctx = (EPS_BLOPEX*)eps->data;
388:   PetscBool      isascii;

391:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
392:   if (isascii) {
393:     PetscViewerASCIIPrintf(viewer,"  block size %D\n",ctx->bs);
394:   }
395:   return(0);
396: }

398: PetscErrorCode EPSSetFromOptions_BLOPEX(PetscOptionItems *PetscOptionsObject,EPS eps)
399: {
401:   PetscBool      flg;
402:   PetscInt       bs;

405:   PetscOptionsHead(PetscOptionsObject,"EPS BLOPEX Options");

407:     PetscOptionsInt("-eps_blopex_blocksize","Block size","EPSBLOPEXSetBlockSize",20,&bs,&flg);
408:     if (flg) { EPSBLOPEXSetBlockSize(eps,bs); }

410:   PetscOptionsTail();

412:   LOBPCG_SetFromOptionsRandomContext();
413:   return(0);
414: }

416: SLEPC_EXTERN PetscErrorCode EPSCreate_BLOPEX(EPS eps)
417: {
418:   EPS_BLOPEX     *ctx;

422:   PetscNewLog(eps,&ctx);
423:   eps->data = (void*)ctx;

425:   eps->categ = EPS_CATEGORY_PRECOND;

427:   eps->ops->solve          = EPSSolve_BLOPEX;
428:   eps->ops->setup          = EPSSetUp_BLOPEX;
429:   eps->ops->setupsort      = EPSSetUpSort_Basic;
430:   eps->ops->setfromoptions = EPSSetFromOptions_BLOPEX;
431:   eps->ops->destroy        = EPSDestroy_BLOPEX;
432:   eps->ops->reset          = EPSReset_BLOPEX;
433:   eps->ops->view           = EPSView_BLOPEX;
434:   eps->ops->backtransform  = EPSBackTransform_Default;
435:   eps->ops->setdefaultst   = EPSSetDefaultST_GMRES;

437:   LOBPCG_InitRandomContext(PetscObjectComm((PetscObject)eps),NULL);
438:   PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXSetBlockSize_C",EPSBLOPEXSetBlockSize_BLOPEX);
439:   PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXGetBlockSize_C",EPSBLOPEXGetBlockSize_BLOPEX);
440:   if (slepc_blopex_useconstr < 0) { PetscObjectComposedDataRegister(&slepc_blopex_useconstr); }
441:   return(0);
442: }