Actual source code: blopex.c
slepc-3.18.2 2023-01-26
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: return 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) BVCreateVec(eps->V,&blopex->w);
154: #if defined(PETSC_USE_COMPLEX)
155: blopex->blap_fn.zpotrf = PETSC_zpotrf_interface;
156: blopex->blap_fn.zhegv = PETSC_zsygv_interface;
157: #else
158: blopex->blap_fn.dpotrf = PETSC_dpotrf_interface;
159: blopex->blap_fn.dsygv = PETSC_dsygv_interface;
160: #endif
161: return 0;
162: }
164: PetscErrorCode EPSSolve_BLOPEX(EPS eps)
165: {
166: EPS_BLOPEX *blopex = (EPS_BLOPEX*)eps->data;
167: PetscScalar sigma,*eigr=NULL;
168: PetscReal *errest=NULL;
169: int i,j,info,its,nconv;
170: double *residhist=NULL;
171: mv_MultiVectorPtr eigenvectors,constraints;
172: #if defined(PETSC_USE_COMPLEX)
173: komplex *lambda=NULL,*lambdahist=NULL;
174: #else
175: double *lambda=NULL,*lambdahist=NULL;
176: #endif
178: STGetShift(eps->st,&sigma);
179: PetscMalloc1(blopex->bs,&lambda);
180: 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);
182: /* Complete the initial basis with random vectors */
183: for (i=0;i<eps->nini;i++) { /* in case the initial vectors were also set with VecSetRandom */
184: BVSetRandomColumn(eps->V,eps->nini);
185: }
186: for (i=eps->nini;i<eps->ncv;i++) BVSetRandomColumn(eps->V,i);
188: while (eps->reason == EPS_CONVERGED_ITERATING) {
190: /* Create multivector of constraints from leading columns of V */
191: PetscObjectComposedDataSetInt((PetscObject)eps->V,slepc_blopex_useconstr,1);
192: BVSetActiveColumns(eps->V,0,eps->nconv);
193: constraints = mv_MultiVectorCreateFromSampleVector(&blopex->ii,eps->nds+eps->nconv,eps->V);
195: /* Create multivector where eigenvectors of this run will be stored */
196: PetscObjectComposedDataSetInt((PetscObject)eps->V,slepc_blopex_useconstr,0);
197: BVSetActiveColumns(eps->V,eps->nconv,eps->nconv+blopex->bs);
198: eigenvectors = mv_MultiVectorCreateFromSampleVector(&blopex->ii,blopex->bs,eps->V);
200: #if defined(PETSC_USE_COMPLEX)
201: info = lobpcg_solve_complex(eigenvectors,blopex,OperatorAMultiVector,
202: eps->isgeneralized?blopex:NULL,eps->isgeneralized?OperatorBMultiVector:NULL,
203: blopex,Precond_FnMultiVector,constraints,
204: blopex->blap_fn,blopex->tol,eps->max_it,0,&its,
205: lambda,lambdahist,blopex->bs,eps->errest+eps->nconv,residhist,blopex->bs);
206: #else
207: info = lobpcg_solve_double(eigenvectors,blopex,OperatorAMultiVector,
208: eps->isgeneralized?blopex:NULL,eps->isgeneralized?OperatorBMultiVector:NULL,
209: blopex,Precond_FnMultiVector,constraints,
210: blopex->blap_fn,blopex->tol,eps->max_it,0,&its,
211: lambda,lambdahist,blopex->bs,eps->errest+eps->nconv,residhist,blopex->bs);
212: #endif
214: mv_MultiVectorDestroy(constraints);
215: mv_MultiVectorDestroy(eigenvectors);
217: for (j=0;j<blopex->bs;j++) {
218: #if defined(PETSC_USE_COMPLEX)
219: eps->eigr[eps->nconv+j] = PetscCMPLX(lambda[j].real,lambda[j].imag);
220: #else
221: eps->eigr[eps->nconv+j] = lambda[j];
222: #endif
223: }
225: if (eps->numbermonitors>0) {
226: for (i=0;i<its;i++) {
227: nconv = 0;
228: for (j=0;j<blopex->bs;j++) {
229: #if defined(PETSC_USE_COMPLEX)
230: eigr[eps->nconv+j] = PetscCMPLX(lambdahist[j+i*blopex->bs].real,lambdahist[j+i*blopex->bs].imag);
231: #else
232: eigr[eps->nconv+j] = lambdahist[j+i*blopex->bs];
233: #endif
234: errest[eps->nconv+j] = residhist[j+i*blopex->bs];
235: if (residhist[j+i*blopex->bs]<=eps->tol) nconv++;
236: }
237: EPSMonitor(eps,eps->its+i,eps->nconv+nconv,eigr,eps->eigi,errest,eps->nconv+blopex->bs);
238: }
239: }
241: eps->its += its;
242: if (info==-1) {
243: eps->reason = EPS_DIVERGED_ITS;
244: break;
245: } else {
246: for (i=0;i<blopex->bs;i++) {
247: if (sigma != 0.0) eps->eigr[eps->nconv+i] += sigma;
248: }
249: eps->nconv += blopex->bs;
250: if (eps->nconv>=eps->nev) eps->reason = EPS_CONVERGED_TOL;
251: }
252: }
254: PetscFree(lambda);
255: if (eps->numbermonitors>0) PetscFree4(lambdahist,eigr,residhist,errest);
256: return 0;
257: }
259: static PetscErrorCode EPSBLOPEXSetBlockSize_BLOPEX(EPS eps,PetscInt bs)
260: {
261: EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;
263: if (bs==PETSC_DEFAULT) {
264: ctx->bs = 0;
265: eps->state = EPS_STATE_INITIAL;
266: } else {
268: ctx->bs = bs;
269: }
270: return 0;
271: }
273: /*@
274: EPSBLOPEXSetBlockSize - Sets the block size of the BLOPEX solver.
276: Logically Collective on eps
278: Input Parameters:
279: + eps - the eigenproblem solver context
280: - bs - the block size
282: Options Database Key:
283: . -eps_blopex_blocksize - Sets the block size
285: Level: advanced
287: .seealso: EPSBLOPEXGetBlockSize()
288: @*/
289: PetscErrorCode EPSBLOPEXSetBlockSize(EPS eps,PetscInt bs)
290: {
293: PetscTryMethod(eps,"EPSBLOPEXSetBlockSize_C",(EPS,PetscInt),(eps,bs));
294: return 0;
295: }
297: static PetscErrorCode EPSBLOPEXGetBlockSize_BLOPEX(EPS eps,PetscInt *bs)
298: {
299: EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;
301: *bs = ctx->bs;
302: return 0;
303: }
305: /*@
306: EPSBLOPEXGetBlockSize - Gets the block size used in the BLOPEX solver.
308: Not Collective
310: Input Parameter:
311: . eps - the eigenproblem solver context
313: Output Parameter:
314: . bs - the block size
316: Level: advanced
318: .seealso: EPSBLOPEXSetBlockSize()
319: @*/
320: PetscErrorCode EPSBLOPEXGetBlockSize(EPS eps,PetscInt *bs)
321: {
324: PetscUseMethod(eps,"EPSBLOPEXGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
325: return 0;
326: }
328: PetscErrorCode EPSReset_BLOPEX(EPS eps)
329: {
330: EPS_BLOPEX *blopex = (EPS_BLOPEX*)eps->data;
332: VecDestroy(&blopex->w);
333: return 0;
334: }
336: PetscErrorCode EPSDestroy_BLOPEX(EPS eps)
337: {
338: LOBPCG_DestroyRandomContext();
339: PetscFree(eps->data);
340: PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXSetBlockSize_C",NULL);
341: PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXGetBlockSize_C",NULL);
342: return 0;
343: }
345: PetscErrorCode EPSView_BLOPEX(EPS eps,PetscViewer viewer)
346: {
347: EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;
348: PetscBool isascii;
350: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
351: if (isascii) PetscViewerASCIIPrintf(viewer," block size %" PetscInt_FMT "\n",ctx->bs);
352: return 0;
353: }
355: PetscErrorCode EPSSetFromOptions_BLOPEX(EPS eps,PetscOptionItems *PetscOptionsObject)
356: {
357: PetscBool flg;
358: PetscInt bs;
360: PetscOptionsHeadBegin(PetscOptionsObject,"EPS BLOPEX Options");
362: PetscOptionsInt("-eps_blopex_blocksize","Block size","EPSBLOPEXSetBlockSize",20,&bs,&flg);
363: if (flg) EPSBLOPEXSetBlockSize(eps,bs);
365: PetscOptionsHeadEnd();
367: LOBPCG_SetFromOptionsRandomContext();
368: return 0;
369: }
371: SLEPC_EXTERN PetscErrorCode EPSCreate_BLOPEX(EPS eps)
372: {
373: EPS_BLOPEX *ctx;
375: PetscNew(&ctx);
376: eps->data = (void*)ctx;
378: eps->categ = EPS_CATEGORY_PRECOND;
380: eps->ops->solve = EPSSolve_BLOPEX;
381: eps->ops->setup = EPSSetUp_BLOPEX;
382: eps->ops->setupsort = EPSSetUpSort_Basic;
383: eps->ops->setfromoptions = EPSSetFromOptions_BLOPEX;
384: eps->ops->destroy = EPSDestroy_BLOPEX;
385: eps->ops->reset = EPSReset_BLOPEX;
386: eps->ops->view = EPSView_BLOPEX;
387: eps->ops->backtransform = EPSBackTransform_Default;
388: eps->ops->setdefaultst = EPSSetDefaultST_GMRES;
390: LOBPCG_InitRandomContext(PetscObjectComm((PetscObject)eps),NULL);
391: PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXSetBlockSize_C",EPSBLOPEXSetBlockSize_BLOPEX);
392: PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXGetBlockSize_C",EPSBLOPEXGetBlockSize_BLOPEX);
393: if (slepc_blopex_useconstr < 0) PetscObjectComposedDataRegister(&slepc_blopex_useconstr);
394: return 0;
395: }