Actual source code: blopex.c

slepc-3.22.1 2024-10-28
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_DETERMINE) { /* 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_DETERMINE) *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:   EPSCheckNotStructured(eps);
129:   if (!blopex->bs) blopex->bs = PetscMin(16,eps->nev);
130:   PetscCall(EPSSetDimensions_BLOPEX(eps,eps->nev,&eps->ncv,&eps->mpd));
131:   if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
132:   if (!eps->which) eps->which = EPS_SMALLEST_REAL;
133:   PetscCheck(eps->which==EPS_SMALLEST_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only smallest real eigenvalues");
134:   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_STOPPING);
135:   EPSCheckIgnored(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_EXTRACTION);

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

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

148:   SLEPCSetupInterpreter(&blopex->ii);

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

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

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

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

187:   PetscFunctionBegin;
188:   PetscCall(STGetShift(eps->st,&sigma));
189:   PetscCall(PetscMalloc1(blopex->bs,&lambda));
190:   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));

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

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

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

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

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

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

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

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

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

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

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

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

287:    Logically Collective

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

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

296:    Level: advanced

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

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

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

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

321:    Not Collective

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

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

329:    Level: advanced

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

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

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

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

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

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

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

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

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

383:   PetscOptionsHeadEnd();

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

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

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

397:   eps->categ = EPS_CATEGORY_PRECOND;

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

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