Actual source code: blopex.c
slepc-main 2024-11-09
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: }