Actual source code: blopex.c
slepc-3.20.2 2024-03-15
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: }