Actual source code: blopex.c

slepc-3.20.2 2024-03-15
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:   PetscFunctionBegin;
 39:   PetscCallAbort(comm,STGetKSP(blopex->st,&ksp));
 40:   PetscCallAbort(comm,KSPSolve(ksp,(Vec)x,(Vec)y));
 41:   PetscFunctionReturnVoid();
 42: }

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

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

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

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

 75: static void OperatorAMultiVector(void *data,void *x,void *y)
 76: {
 77:   EPS_BLOPEX *blopex = (EPS_BLOPEX*)data;

 79:   PetscFunctionBegin;
 80:   blopex->ii.Eval(OperatorASingleVector,data,x,y);
 81:   PetscFunctionReturnVoid();
 82: }

 84: static void OperatorBSingleVector(void *data,void *x,void *y)
 85: {
 86:   EPS_BLOPEX     *blopex = (EPS_BLOPEX*)data;
 87:   MPI_Comm       comm = PetscObjectComm((PetscObject)blopex->st);
 88:   Mat            B;

 90:   PetscFunctionBegin;
 91:   PetscCallAbort(comm,STGetMatrix(blopex->st,1,&B));
 92:   PetscCallAbort(comm,MatMult(B,(Vec)x,(Vec)y));
 93:   PetscFunctionReturnVoid();
 94: }

 96: static void OperatorBMultiVector(void *data,void *x,void *y)
 97: {
 98:   EPS_BLOPEX *blopex = (EPS_BLOPEX*)data;

100:   PetscFunctionBegin;
101:   blopex->ii.Eval(OperatorBSingleVector,data,x,y);
102:   PetscFunctionReturnVoid();
103: }

105: static PetscErrorCode EPSSetDimensions_BLOPEX(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
106: {
107:   EPS_BLOPEX     *ctx = (EPS_BLOPEX*)eps->data;
108:   PetscInt       k;

110:   PetscFunctionBegin;
111:   k = ((eps->nev-1)/ctx->bs+1)*ctx->bs;
112:   if (*ncv!=PETSC_DEFAULT) { /* ncv set */
113:     PetscCheck(*ncv>=k,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv is not sufficiently large");
114:   } else *ncv = k;
115:   if (*mpd==PETSC_DEFAULT) *mpd = *ncv;
116:   else PetscCall(PetscInfo(eps,"Warning: given value of mpd ignored\n"));
117:   PetscFunctionReturn(PETSC_SUCCESS);
118: }

120: static PetscErrorCode EPSSetUp_BLOPEX(EPS eps)
121: {
122:   EPS_BLOPEX     *blopex = (EPS_BLOPEX*)eps->data;
123:   PetscBool      flg;
124:   KSP            ksp;

126:   PetscFunctionBegin;
127:   EPSCheckHermitianDefinite(eps);
128:   if (!blopex->bs) blopex->bs = PetscMin(16,eps->nev);
129:   PetscCall(EPSSetDimensions_BLOPEX(eps,eps->nev,&eps->ncv,&eps->mpd));
130:   if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
131:   if (!eps->which) eps->which = EPS_SMALLEST_REAL;
132:   PetscCheck(eps->which==EPS_SMALLEST_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only smallest real eigenvalues");
133:   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_STOPPING);
134:   EPSCheckIgnored(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_EXTRACTION);

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

138:   PetscCheck(eps->converged==EPSConvergedRelative || eps->converged==EPSConvergedAbsolute,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Convergence test not supported in this solver");
139:   if (eps->converged == EPSConvergedRelative) {
140:     blopex->tol.absolute = 0.0;
141:     blopex->tol.relative = SlepcDefaultTol(eps->tol);
142:   } else {  /* EPSConvergedAbsolute */
143:     blopex->tol.absolute = SlepcDefaultTol(eps->tol);
144:     blopex->tol.relative = 0.0;
145:   }

147:   SLEPCSetupInterpreter(&blopex->ii);

149:   PetscCall(STGetKSP(eps->st,&ksp));
150:   PetscCall(PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg));
151:   if (!flg) PetscCall(PetscInfo(eps,"Warning: ignoring KSP, should use KSPPREONLY\n"));

153:   /* allocate memory */
154:   if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
155:   PetscCall(PetscObjectTypeCompareAny((PetscObject)eps->V,&flg,BVVECS,BVCONTIGUOUS,""));
156:   if (!flg) {  /* blopex only works with BVVECS or BVCONTIGUOUS */
157:     PetscCall(BVSetType(eps->V,BVCONTIGUOUS));
158:   }
159:   PetscCall(EPSAllocateSolution(eps,0));
160:   if (!blopex->w) PetscCall(BVCreateVec(eps->V,&blopex->w));

162: #if defined(PETSC_USE_COMPLEX)
163:   blopex->blap_fn.zpotrf = PETSC_zpotrf_interface;
164:   blopex->blap_fn.zhegv = PETSC_zsygv_interface;
165: #else
166:   blopex->blap_fn.dpotrf = PETSC_dpotrf_interface;
167:   blopex->blap_fn.dsygv = PETSC_dsygv_interface;
168: #endif
169:   PetscFunctionReturn(PETSC_SUCCESS);
170: }

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

186:   PetscFunctionBegin;
187:   PetscCall(STGetShift(eps->st,&sigma));
188:   PetscCall(PetscMalloc1(blopex->bs,&lambda));
189:   if (eps->numbermonitors>0) PetscCall(PetscMalloc4(blopex->bs*(eps->max_it+1),&lambdahist,eps->ncv,&eigr,blopex->bs*(eps->max_it+1),&residhist,eps->ncv,&errest));

191:   /* Complete the initial basis with random vectors */
192:   for (i=0;i<eps->nini;i++) {  /* in case the initial vectors were also set with VecSetRandom */
193:     PetscCall(BVSetRandomColumn(eps->V,eps->nini));
194:   }
195:   for (i=eps->nini;i<eps->ncv;i++) PetscCall(BVSetRandomColumn(eps->V,i));

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

199:     /* Create multivector of constraints from leading columns of V */
200:     PetscCall(PetscObjectComposedDataSetInt((PetscObject)eps->V,slepc_blopex_useconstr,1));
201:     PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
202:     constraints = mv_MultiVectorCreateFromSampleVector(&blopex->ii,eps->nds+eps->nconv,eps->V);

204:     /* Create multivector where eigenvectors of this run will be stored */
205:     PetscCall(PetscObjectComposedDataSetInt((PetscObject)eps->V,slepc_blopex_useconstr,0));
206:     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,eps->nconv+blopex->bs));
207:     eigenvectors = mv_MultiVectorCreateFromSampleVector(&blopex->ii,blopex->bs,eps->V);

209: #if defined(PETSC_USE_COMPLEX)
210:     info = lobpcg_solve_complex(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: #else
216:     info = lobpcg_solve_double(eigenvectors,blopex,OperatorAMultiVector,
217:           eps->isgeneralized?blopex:NULL,eps->isgeneralized?OperatorBMultiVector:NULL,
218:           blopex,Precond_FnMultiVector,constraints,
219:           blopex->blap_fn,blopex->tol,eps->max_it,0,&its,
220:           lambda,lambdahist,blopex->bs,eps->errest+eps->nconv,residhist,blopex->bs);
221: #endif
222:     PetscCheck(info==0,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"BLOPEX failed with exit code=%d",info);
223:     mv_MultiVectorDestroy(constraints);
224:     mv_MultiVectorDestroy(eigenvectors);

226:     for (j=0;j<blopex->bs;j++) {
227: #if defined(PETSC_USE_COMPLEX)
228:       eps->eigr[eps->nconv+j] = PetscCMPLX(lambda[j].real,lambda[j].imag);
229: #else
230:       eps->eigr[eps->nconv+j] = lambda[j];
231: #endif
232:     }

234:     if (eps->numbermonitors>0) {
235:       for (i=0;i<its;i++) {
236:         nconv = 0;
237:         for (j=0;j<blopex->bs;j++) {
238: #if defined(PETSC_USE_COMPLEX)
239:           eigr[eps->nconv+j] = PetscCMPLX(lambdahist[j+i*blopex->bs].real,lambdahist[j+i*blopex->bs].imag);
240: #else
241:           eigr[eps->nconv+j] = lambdahist[j+i*blopex->bs];
242: #endif
243:           errest[eps->nconv+j] = residhist[j+i*blopex->bs];
244:           if (residhist[j+i*blopex->bs]<=eps->tol) nconv++;
245:         }
246:         PetscCall(EPSMonitor(eps,eps->its+i,eps->nconv+nconv,eigr,eps->eigi,errest,eps->nconv+blopex->bs));
247:       }
248:     }

250:     eps->its += its;
251:     if (info==-1) {
252:       eps->reason = EPS_DIVERGED_ITS;
253:       break;
254:     } else {
255:       for (i=0;i<blopex->bs;i++) {
256:         if (sigma != 0.0) eps->eigr[eps->nconv+i] += sigma;
257:       }
258:       eps->nconv += blopex->bs;
259:       if (eps->nconv>=eps->nev) eps->reason = EPS_CONVERGED_TOL;
260:     }
261:   }

263:   PetscCall(PetscFree(lambda));
264:   if (eps->numbermonitors>0) PetscCall(PetscFree4(lambdahist,eigr,residhist,errest));
265:   PetscFunctionReturn(PETSC_SUCCESS);
266: }

268: static PetscErrorCode EPSBLOPEXSetBlockSize_BLOPEX(EPS eps,PetscInt bs)
269: {
270:   EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;

272:   PetscFunctionBegin;
273:   if (bs==PETSC_DEFAULT) {
274:     ctx->bs    = 0;
275:     eps->state = EPS_STATE_INITIAL;
276:   } else {
277:     PetscCheck(bs>0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block size must be >0");
278:     ctx->bs = bs;
279:   }
280:   PetscFunctionReturn(PETSC_SUCCESS);
281: }

283: /*@
284:    EPSBLOPEXSetBlockSize - Sets the block size of the BLOPEX solver.

286:    Logically Collective

288:    Input Parameters:
289: +  eps - the eigenproblem solver context
290: -  bs  - the block size

292:    Options Database Key:
293: .  -eps_blopex_blocksize - Sets the block size

295:    Level: advanced

297: .seealso: EPSBLOPEXGetBlockSize()
298: @*/
299: PetscErrorCode EPSBLOPEXSetBlockSize(EPS eps,PetscInt bs)
300: {
301:   PetscFunctionBegin;
304:   PetscTryMethod(eps,"EPSBLOPEXSetBlockSize_C",(EPS,PetscInt),(eps,bs));
305:   PetscFunctionReturn(PETSC_SUCCESS);
306: }

308: static PetscErrorCode EPSBLOPEXGetBlockSize_BLOPEX(EPS eps,PetscInt *bs)
309: {
310:   EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;

312:   PetscFunctionBegin;
313:   *bs = ctx->bs;
314:   PetscFunctionReturn(PETSC_SUCCESS);
315: }

317: /*@
318:    EPSBLOPEXGetBlockSize - Gets the block size used in the BLOPEX solver.

320:    Not Collective

322:    Input Parameter:
323: .  eps - the eigenproblem solver context

325:    Output Parameter:
326: .  bs - the block size

328:    Level: advanced

330: .seealso: EPSBLOPEXSetBlockSize()
331: @*/
332: PetscErrorCode EPSBLOPEXGetBlockSize(EPS eps,PetscInt *bs)
333: {
334:   PetscFunctionBegin;
336:   PetscAssertPointer(bs,2);
337:   PetscUseMethod(eps,"EPSBLOPEXGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
338:   PetscFunctionReturn(PETSC_SUCCESS);
339: }

341: static PetscErrorCode EPSReset_BLOPEX(EPS eps)
342: {
343:   EPS_BLOPEX     *blopex = (EPS_BLOPEX*)eps->data;

345:   PetscFunctionBegin;
346:   PetscCall(VecDestroy(&blopex->w));
347:   PetscFunctionReturn(PETSC_SUCCESS);
348: }

350: static PetscErrorCode EPSDestroy_BLOPEX(EPS eps)
351: {
352:   PetscFunctionBegin;
353:   LOBPCG_DestroyRandomContext();
354:   PetscCall(PetscFree(eps->data));
355:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXSetBlockSize_C",NULL));
356:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXGetBlockSize_C",NULL));
357:   PetscFunctionReturn(PETSC_SUCCESS);
358: }

360: static PetscErrorCode EPSView_BLOPEX(EPS eps,PetscViewer viewer)
361: {
362:   EPS_BLOPEX     *ctx = (EPS_BLOPEX*)eps->data;
363:   PetscBool      isascii;

365:   PetscFunctionBegin;
366:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
367:   if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer,"  block size %" PetscInt_FMT "\n",ctx->bs));
368:   PetscFunctionReturn(PETSC_SUCCESS);
369: }

371: static PetscErrorCode EPSSetFromOptions_BLOPEX(EPS eps,PetscOptionItems *PetscOptionsObject)
372: {
373:   PetscBool      flg;
374:   PetscInt       bs;

376:   PetscFunctionBegin;
377:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS BLOPEX Options");

379:     PetscCall(PetscOptionsInt("-eps_blopex_blocksize","Block size","EPSBLOPEXSetBlockSize",20,&bs,&flg));
380:     if (flg) PetscCall(EPSBLOPEXSetBlockSize(eps,bs));

382:   PetscOptionsHeadEnd();

384:   LOBPCG_SetFromOptionsRandomContext();
385:   PetscFunctionReturn(PETSC_SUCCESS);
386: }

388: SLEPC_EXTERN PetscErrorCode EPSCreate_BLOPEX(EPS eps)
389: {
390:   EPS_BLOPEX     *ctx;

392:   PetscFunctionBegin;
393:   PetscCall(PetscNew(&ctx));
394:   eps->data = (void*)ctx;

396:   eps->categ = EPS_CATEGORY_PRECOND;

398:   eps->ops->solve          = EPSSolve_BLOPEX;
399:   eps->ops->setup          = EPSSetUp_BLOPEX;
400:   eps->ops->setupsort      = EPSSetUpSort_Basic;
401:   eps->ops->setfromoptions = EPSSetFromOptions_BLOPEX;
402:   eps->ops->destroy        = EPSDestroy_BLOPEX;
403:   eps->ops->reset          = EPSReset_BLOPEX;
404:   eps->ops->view           = EPSView_BLOPEX;
405:   eps->ops->backtransform  = EPSBackTransform_Default;
406:   eps->ops->setdefaultst   = EPSSetDefaultST_GMRES;

408:   LOBPCG_InitRandomContext(PetscObjectComm((PetscObject)eps),NULL);
409:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXSetBlockSize_C",EPSBLOPEXSetBlockSize_BLOPEX));
410:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXGetBlockSize_C",EPSBLOPEXGetBlockSize_BLOPEX));
411:   if (slepc_blopex_useconstr < 0) PetscCall(PetscObjectComposedDataRegister(&slepc_blopex_useconstr));
412:   PetscFunctionReturn(PETSC_SUCCESS);
413: }