Actual source code: blopex.c

slepc-3.18.3 2023-03-24
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:    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: {
 34:   EPS_BLOPEX     *blopex = (EPS_BLOPEX*)data;
 35:   MPI_Comm       comm = PetscObjectComm((PetscObject)blopex->st);
 36:   KSP            ksp;

 38:   comm,STGetKSP(blopex->st,&ksp);
 39:   comm,KSPSolve(ksp,(Vec)x,(Vec)y);
 40:   return;
 41: }

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

 47:   blopex->ii.Eval(Precond_FnSingleVector,data,x,y);
 48:   return;
 49: }

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

 59:   comm,STGetNumMatrices(blopex->st,&nmat);
 60:   comm,STGetMatrix(blopex->st,0,&A);
 61:   if (nmat>1) comm,STGetMatrix(blopex->st,1,&B);
 62:   comm,MatMult(A,(Vec)x,(Vec)y);
 63:   comm,STGetShift(blopex->st,&sigma);
 64:   if (sigma != 0.0) {
 65:     if (nmat>1) comm,MatMult(B,(Vec)x,blopex->w);
 66:     else comm,VecCopy((Vec)x,blopex->w);
 67:     comm,VecAXPY((Vec)y,-sigma,blopex->w);
 68:   }
 69:   return;
 70: }

 72: static void OperatorAMultiVector(void *data,void *x,void *y)
 73: {
 74:   EPS_BLOPEX *blopex = (EPS_BLOPEX*)data;

 76:   blopex->ii.Eval(OperatorASingleVector,data,x,y);
 77:   return;
 78: }

 80: static void OperatorBSingleVector(void *data,void *x,void *y)
 81: {
 82:   EPS_BLOPEX     *blopex = (EPS_BLOPEX*)data;
 83:   MPI_Comm       comm = PetscObjectComm((PetscObject)blopex->st);
 84:   Mat            B;

 86:   comm,STGetMatrix(blopex->st,1,&B);
 87:   comm,MatMult(B,(Vec)x,(Vec)y);
 88:   return;
 89: }

 91: static void OperatorBMultiVector(void *data,void *x,void *y)
 92: {
 93:   EPS_BLOPEX *blopex = (EPS_BLOPEX*)data;

 95:   blopex->ii.Eval(OperatorBSingleVector,data,x,y);
 96:   return;
 97: }

 99: PetscErrorCode EPSSetDimensions_BLOPEX(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
100: {
101:   EPS_BLOPEX     *ctx = (EPS_BLOPEX*)eps->data;
102:   PetscInt       k;

104:   k = ((eps->nev-1)/ctx->bs+1)*ctx->bs;
105:   if (*ncv!=PETSC_DEFAULT) { /* ncv set */
107:   } else *ncv = k;
108:   if (*mpd==PETSC_DEFAULT) *mpd = *ncv;
109:   else PetscInfo(eps,"Warning: given value of mpd ignored\n");
110:   return 0;
111: }

113: PetscErrorCode EPSSetUp_BLOPEX(EPS eps)
114: {
115:   EPS_BLOPEX     *blopex = (EPS_BLOPEX*)eps->data;
116:   PetscBool      flg;
117:   KSP            ksp;

119:   EPSCheckHermitianDefinite(eps);
120:   if (!blopex->bs) blopex->bs = PetscMin(16,eps->nev);
121:   EPSSetDimensions_BLOPEX(eps,eps->nev,&eps->ncv,&eps->mpd);
122:   if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
123:   if (!eps->which) eps->which = EPS_SMALLEST_REAL;
125:   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_STOPPING);
126:   EPSCheckIgnored(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_EXTRACTION);

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

131:   if (eps->converged == EPSConvergedRelative) {
132:     blopex->tol.absolute = 0.0;
133:     blopex->tol.relative = SlepcDefaultTol(eps->tol);
134:   } else {  /* EPSConvergedAbsolute */
135:     blopex->tol.absolute = SlepcDefaultTol(eps->tol);
136:     blopex->tol.relative = 0.0;
137:   }

139:   SLEPCSetupInterpreter(&blopex->ii);

141:   STGetKSP(eps->st,&ksp);
142:   PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);
143:   if (!flg) PetscInfo(eps,"Warning: ignoring KSP, should use KSPPREONLY\n");

145:   /* allocate memory */
146:   if (!eps->V) EPSGetBV(eps,&eps->V);
147:   PetscObjectTypeCompareAny((PetscObject)eps->V,&flg,BVVECS,BVCONTIGUOUS,"");
148:   if (!flg) {  /* blopex only works with BVVECS or BVCONTIGUOUS */
149:     BVSetType(eps->V,BVCONTIGUOUS);
150:   }
151:   EPSAllocateSolution(eps,0);
152:   if (!blopex->w) BVCreateVec(eps->V,&blopex->w);

154: #if defined(PETSC_USE_COMPLEX)
155:   blopex->blap_fn.zpotrf = PETSC_zpotrf_interface;
156:   blopex->blap_fn.zhegv = PETSC_zsygv_interface;
157: #else
158:   blopex->blap_fn.dpotrf = PETSC_dpotrf_interface;
159:   blopex->blap_fn.dsygv = PETSC_dsygv_interface;
160: #endif
161:   return 0;
162: }

164: PetscErrorCode EPSSolve_BLOPEX(EPS eps)
165: {
166:   EPS_BLOPEX        *blopex = (EPS_BLOPEX*)eps->data;
167:   PetscScalar       sigma,*eigr=NULL;
168:   PetscReal         *errest=NULL;
169:   int               i,j,info,its,nconv;
170:   double            *residhist=NULL;
171:   mv_MultiVectorPtr eigenvectors,constraints;
172: #if defined(PETSC_USE_COMPLEX)
173:   komplex           *lambda=NULL,*lambdahist=NULL;
174: #else
175:   double            *lambda=NULL,*lambdahist=NULL;
176: #endif

178:   STGetShift(eps->st,&sigma);
179:   PetscMalloc1(blopex->bs,&lambda);
180:   if (eps->numbermonitors>0) PetscMalloc4(blopex->bs*(eps->max_it+1),&lambdahist,eps->ncv,&eigr,blopex->bs*(eps->max_it+1),&residhist,eps->ncv,&errest);

182:   /* Complete the initial basis with random vectors */
183:   for (i=0;i<eps->nini;i++) {  /* in case the initial vectors were also set with VecSetRandom */
184:     BVSetRandomColumn(eps->V,eps->nini);
185:   }
186:   for (i=eps->nini;i<eps->ncv;i++) BVSetRandomColumn(eps->V,i);

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

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

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

200: #if defined(PETSC_USE_COMPLEX)
201:     info = lobpcg_solve_complex(eigenvectors,blopex,OperatorAMultiVector,
202:           eps->isgeneralized?blopex:NULL,eps->isgeneralized?OperatorBMultiVector:NULL,
203:           blopex,Precond_FnMultiVector,constraints,
204:           blopex->blap_fn,blopex->tol,eps->max_it,0,&its,
205:           lambda,lambdahist,blopex->bs,eps->errest+eps->nconv,residhist,blopex->bs);
206: #else
207:     info = lobpcg_solve_double(eigenvectors,blopex,OperatorAMultiVector,
208:           eps->isgeneralized?blopex:NULL,eps->isgeneralized?OperatorBMultiVector:NULL,
209:           blopex,Precond_FnMultiVector,constraints,
210:           blopex->blap_fn,blopex->tol,eps->max_it,0,&its,
211:           lambda,lambdahist,blopex->bs,eps->errest+eps->nconv,residhist,blopex->bs);
212: #endif
214:     mv_MultiVectorDestroy(constraints);
215:     mv_MultiVectorDestroy(eigenvectors);

217:     for (j=0;j<blopex->bs;j++) {
218: #if defined(PETSC_USE_COMPLEX)
219:       eps->eigr[eps->nconv+j] = PetscCMPLX(lambda[j].real,lambda[j].imag);
220: #else
221:       eps->eigr[eps->nconv+j] = lambda[j];
222: #endif
223:     }

225:     if (eps->numbermonitors>0) {
226:       for (i=0;i<its;i++) {
227:         nconv = 0;
228:         for (j=0;j<blopex->bs;j++) {
229: #if defined(PETSC_USE_COMPLEX)
230:           eigr[eps->nconv+j] = PetscCMPLX(lambdahist[j+i*blopex->bs].real,lambdahist[j+i*blopex->bs].imag);
231: #else
232:           eigr[eps->nconv+j] = lambdahist[j+i*blopex->bs];
233: #endif
234:           errest[eps->nconv+j] = residhist[j+i*blopex->bs];
235:           if (residhist[j+i*blopex->bs]<=eps->tol) nconv++;
236:         }
237:         EPSMonitor(eps,eps->its+i,eps->nconv+nconv,eigr,eps->eigi,errest,eps->nconv+blopex->bs);
238:       }
239:     }

241:     eps->its += its;
242:     if (info==-1) {
243:       eps->reason = EPS_DIVERGED_ITS;
244:       break;
245:     } else {
246:       for (i=0;i<blopex->bs;i++) {
247:         if (sigma != 0.0) eps->eigr[eps->nconv+i] += sigma;
248:       }
249:       eps->nconv += blopex->bs;
250:       if (eps->nconv>=eps->nev) eps->reason = EPS_CONVERGED_TOL;
251:     }
252:   }

254:   PetscFree(lambda);
255:   if (eps->numbermonitors>0) PetscFree4(lambdahist,eigr,residhist,errest);
256:   return 0;
257: }

259: static PetscErrorCode EPSBLOPEXSetBlockSize_BLOPEX(EPS eps,PetscInt bs)
260: {
261:   EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;

263:   if (bs==PETSC_DEFAULT) {
264:     ctx->bs    = 0;
265:     eps->state = EPS_STATE_INITIAL;
266:   } else {
268:     ctx->bs = bs;
269:   }
270:   return 0;
271: }

273: /*@
274:    EPSBLOPEXSetBlockSize - Sets the block size of the BLOPEX solver.

276:    Logically Collective on eps

278:    Input Parameters:
279: +  eps - the eigenproblem solver context
280: -  bs  - the block size

282:    Options Database Key:
283: .  -eps_blopex_blocksize - Sets the block size

285:    Level: advanced

287: .seealso: EPSBLOPEXGetBlockSize()
288: @*/
289: PetscErrorCode EPSBLOPEXSetBlockSize(EPS eps,PetscInt bs)
290: {
293:   PetscTryMethod(eps,"EPSBLOPEXSetBlockSize_C",(EPS,PetscInt),(eps,bs));
294:   return 0;
295: }

297: static PetscErrorCode EPSBLOPEXGetBlockSize_BLOPEX(EPS eps,PetscInt *bs)
298: {
299:   EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;

301:   *bs = ctx->bs;
302:   return 0;
303: }

305: /*@
306:    EPSBLOPEXGetBlockSize - Gets the block size used in the BLOPEX solver.

308:    Not Collective

310:    Input Parameter:
311: .  eps - the eigenproblem solver context

313:    Output Parameter:
314: .  bs - the block size

316:    Level: advanced

318: .seealso: EPSBLOPEXSetBlockSize()
319: @*/
320: PetscErrorCode EPSBLOPEXGetBlockSize(EPS eps,PetscInt *bs)
321: {
324:   PetscUseMethod(eps,"EPSBLOPEXGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
325:   return 0;
326: }

328: PetscErrorCode EPSReset_BLOPEX(EPS eps)
329: {
330:   EPS_BLOPEX     *blopex = (EPS_BLOPEX*)eps->data;

332:   VecDestroy(&blopex->w);
333:   return 0;
334: }

336: PetscErrorCode EPSDestroy_BLOPEX(EPS eps)
337: {
338:   LOBPCG_DestroyRandomContext();
339:   PetscFree(eps->data);
340:   PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXSetBlockSize_C",NULL);
341:   PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXGetBlockSize_C",NULL);
342:   return 0;
343: }

345: PetscErrorCode EPSView_BLOPEX(EPS eps,PetscViewer viewer)
346: {
347:   EPS_BLOPEX     *ctx = (EPS_BLOPEX*)eps->data;
348:   PetscBool      isascii;

350:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
351:   if (isascii) PetscViewerASCIIPrintf(viewer,"  block size %" PetscInt_FMT "\n",ctx->bs);
352:   return 0;
353: }

355: PetscErrorCode EPSSetFromOptions_BLOPEX(EPS eps,PetscOptionItems *PetscOptionsObject)
356: {
357:   PetscBool      flg;
358:   PetscInt       bs;

360:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS BLOPEX Options");

362:     PetscOptionsInt("-eps_blopex_blocksize","Block size","EPSBLOPEXSetBlockSize",20,&bs,&flg);
363:     if (flg) EPSBLOPEXSetBlockSize(eps,bs);

365:   PetscOptionsHeadEnd();

367:   LOBPCG_SetFromOptionsRandomContext();
368:   return 0;
369: }

371: SLEPC_EXTERN PetscErrorCode EPSCreate_BLOPEX(EPS eps)
372: {
373:   EPS_BLOPEX     *ctx;

375:   PetscNew(&ctx);
376:   eps->data = (void*)ctx;

378:   eps->categ = EPS_CATEGORY_PRECOND;

380:   eps->ops->solve          = EPSSolve_BLOPEX;
381:   eps->ops->setup          = EPSSetUp_BLOPEX;
382:   eps->ops->setupsort      = EPSSetUpSort_Basic;
383:   eps->ops->setfromoptions = EPSSetFromOptions_BLOPEX;
384:   eps->ops->destroy        = EPSDestroy_BLOPEX;
385:   eps->ops->reset          = EPSReset_BLOPEX;
386:   eps->ops->view           = EPSView_BLOPEX;
387:   eps->ops->backtransform  = EPSBackTransform_Default;
388:   eps->ops->setdefaultst   = EPSSetDefaultST_GMRES;

390:   LOBPCG_InitRandomContext(PetscObjectComm((PetscObject)eps),NULL);
391:   PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXSetBlockSize_C",EPSBLOPEXSetBlockSize_BLOPEX);
392:   PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXGetBlockSize_C",EPSBLOPEXGetBlockSize_BLOPEX);
393:   if (slepc_blopex_useconstr < 0) PetscObjectComposedDataRegister(&slepc_blopex_useconstr);
394:   return 0;
395: }