Actual source code: blopex.c

slepc-3.17.2 2022-08-09
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:   PetscFunctionReturn(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) {
153:     BVCreateVec(eps->V,&blopex->w);
154:     PetscLogObjectParent((PetscObject)eps,(PetscObject)blopex->w);
155:   }

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

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

181:   STGetShift(eps->st,&sigma);
182:   PetscMalloc1(blopex->bs,&lambda);
183:   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);

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

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

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

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

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

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

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

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

257:   PetscFree(lambda);
258:   if (eps->numbermonitors>0) PetscFree4(lambdahist,eigr,residhist,errest);
259:   PetscFunctionReturn(0);
260: }

262: static PetscErrorCode EPSBLOPEXSetBlockSize_BLOPEX(EPS eps,PetscInt bs)
263: {
264:   EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;

266:   if (bs==PETSC_DEFAULT) {
267:     ctx->bs    = 0;
268:     eps->state = EPS_STATE_INITIAL;
269:   } else {
271:     ctx->bs = bs;
272:   }
273:   PetscFunctionReturn(0);
274: }

276: /*@
277:    EPSBLOPEXSetBlockSize - Sets the block size of the BLOPEX solver.

279:    Logically Collective on eps

281:    Input Parameters:
282: +  eps - the eigenproblem solver context
283: -  bs  - the block size

285:    Options Database Key:
286: .  -eps_blopex_blocksize - Sets the block size

288:    Level: advanced

290: .seealso: EPSBLOPEXGetBlockSize()
291: @*/
292: PetscErrorCode EPSBLOPEXSetBlockSize(EPS eps,PetscInt bs)
293: {
296:   PetscTryMethod(eps,"EPSBLOPEXSetBlockSize_C",(EPS,PetscInt),(eps,bs));
297:   PetscFunctionReturn(0);
298: }

300: static PetscErrorCode EPSBLOPEXGetBlockSize_BLOPEX(EPS eps,PetscInt *bs)
301: {
302:   EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;

304:   *bs = ctx->bs;
305:   PetscFunctionReturn(0);
306: }

308: /*@
309:    EPSBLOPEXGetBlockSize - Gets the block size used in the BLOPEX solver.

311:    Not Collective

313:    Input Parameter:
314: .  eps - the eigenproblem solver context

316:    Output Parameter:
317: .  bs - the block size

319:    Level: advanced

321: .seealso: EPSBLOPEXSetBlockSize()
322: @*/
323: PetscErrorCode EPSBLOPEXGetBlockSize(EPS eps,PetscInt *bs)
324: {
327:   PetscUseMethod(eps,"EPSBLOPEXGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
328:   PetscFunctionReturn(0);
329: }

331: PetscErrorCode EPSReset_BLOPEX(EPS eps)
332: {
333:   EPS_BLOPEX     *blopex = (EPS_BLOPEX*)eps->data;

335:   VecDestroy(&blopex->w);
336:   PetscFunctionReturn(0);
337: }

339: PetscErrorCode EPSDestroy_BLOPEX(EPS eps)
340: {
341:   LOBPCG_DestroyRandomContext();
342:   PetscFree(eps->data);
343:   PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXSetBlockSize_C",NULL);
344:   PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXGetBlockSize_C",NULL);
345:   PetscFunctionReturn(0);
346: }

348: PetscErrorCode EPSView_BLOPEX(EPS eps,PetscViewer viewer)
349: {
350:   EPS_BLOPEX     *ctx = (EPS_BLOPEX*)eps->data;
351:   PetscBool      isascii;

353:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
354:   if (isascii) PetscViewerASCIIPrintf(viewer,"  block size %" PetscInt_FMT "\n",ctx->bs);
355:   PetscFunctionReturn(0);
356: }

358: PetscErrorCode EPSSetFromOptions_BLOPEX(PetscOptionItems *PetscOptionsObject,EPS eps)
359: {
360:   PetscBool      flg;
361:   PetscInt       bs;

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

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

368:   PetscOptionsTail();

370:   LOBPCG_SetFromOptionsRandomContext();
371:   PetscFunctionReturn(0);
372: }

374: SLEPC_EXTERN PetscErrorCode EPSCreate_BLOPEX(EPS eps)
375: {
376:   EPS_BLOPEX     *ctx;

378:   PetscNewLog(eps,&ctx);
379:   eps->data = (void*)ctx;

381:   eps->categ = EPS_CATEGORY_PRECOND;

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

393:   LOBPCG_InitRandomContext(PetscObjectComm((PetscObject)eps),NULL);
394:   PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXSetBlockSize_C",EPSBLOPEXSetBlockSize_BLOPEX);
395:   PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXGetBlockSize_C",EPSBLOPEXGetBlockSize_BLOPEX);
396:   if (slepc_blopex_useconstr < 0) PetscObjectComposedDataRegister(&slepc_blopex_useconstr);
397:   PetscFunctionReturn(0);
398: }