Actual source code: blopex.c
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: if (*nev==0) *nev = 1;
112: k = ((*nev-1)/ctx->bs+1)*ctx->bs;
113: if (*ncv!=PETSC_DETERMINE) { /* ncv set */
114: PetscCheck(*ncv>=k,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv is not sufficiently large");
115: } else *ncv = k;
116: if (*mpd==PETSC_DETERMINE) *mpd = *ncv;
117: else PetscCall(PetscInfo(eps,"Warning: given value of mpd ignored\n"));
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
121: static PetscErrorCode EPSSetUp_BLOPEX(EPS eps)
122: {
123: EPS_BLOPEX *blopex = (EPS_BLOPEX*)eps->data;
124: PetscBool flg;
125: KSP ksp;
127: PetscFunctionBegin;
128: EPSCheckHermitianDefinite(eps);
129: EPSCheckNotStructured(eps);
130: if (!blopex->bs) blopex->bs = PetscMin(16,eps->nev?eps->nev:1);
131: PetscCall(EPSSetDimensions_BLOPEX(eps,&eps->nev,&eps->ncv,&eps->mpd));
132: if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
133: if (!eps->which) eps->which = EPS_SMALLEST_REAL;
134: PetscCheck(eps->which==EPS_SMALLEST_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only smallest real eigenvalues");
135: EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_STOPPING);
136: EPSCheckIgnored(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_EXTRACTION);
138: blopex->st = eps->st;
140: PetscCheck(eps->converged==EPSConvergedRelative || eps->converged==EPSConvergedAbsolute,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Convergence test not supported in this solver");
141: if (eps->converged == EPSConvergedRelative) {
142: blopex->tol.absolute = 0.0;
143: blopex->tol.relative = SlepcDefaultTol(eps->tol);
144: } else { /* EPSConvergedAbsolute */
145: blopex->tol.absolute = SlepcDefaultTol(eps->tol);
146: blopex->tol.relative = 0.0;
147: }
149: SLEPCSetupInterpreter(&blopex->ii);
151: PetscCall(STGetKSP(eps->st,&ksp));
152: PetscCall(PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg));
153: if (!flg) PetscCall(PetscInfo(eps,"Warning: ignoring KSP, should use KSPPREONLY\n"));
155: /* allocate memory */
156: if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
157: PetscCall(PetscObjectTypeCompareAny((PetscObject)eps->V,&flg,BVVECS,BVCONTIGUOUS,""));
158: if (!flg) { /* blopex only works with BVVECS or BVCONTIGUOUS */
159: PetscCall(BVSetType(eps->V,BVCONTIGUOUS));
160: }
161: PetscCall(EPSAllocateSolution(eps,0));
162: if (!blopex->w) PetscCall(BVCreateVec(eps->V,&blopex->w));
164: #if defined(PETSC_USE_COMPLEX)
165: blopex->blap_fn.zpotrf = PETSC_zpotrf_interface;
166: blopex->blap_fn.zhegv = PETSC_zsygv_interface;
167: #else
168: blopex->blap_fn.dpotrf = PETSC_dpotrf_interface;
169: blopex->blap_fn.dsygv = PETSC_dsygv_interface;
170: #endif
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: static PetscErrorCode EPSSolve_BLOPEX(EPS eps)
175: {
176: EPS_BLOPEX *blopex = (EPS_BLOPEX*)eps->data;
177: PetscScalar sigma,*eigr=NULL;
178: PetscReal *errest=NULL;
179: int i,j,info,its,nconv;
180: double *residhist=NULL;
181: mv_MultiVectorPtr eigenvectors,constraints;
182: #if defined(PETSC_USE_COMPLEX)
183: komplex *lambda=NULL,*lambdahist=NULL;
184: #else
185: double *lambda=NULL,*lambdahist=NULL;
186: #endif
188: PetscFunctionBegin;
189: PetscCall(STGetShift(eps->st,&sigma));
190: PetscCall(PetscMalloc1(blopex->bs,&lambda));
191: 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));
193: if (eps->nini<eps->ncv) {
194: /* Complete the initial basis with random vectors */
195: for (i=0;i<eps->nini;i++) { /* in case the initial vectors were also set with VecSetRandom */
196: PetscCall(BVSetRandomColumn(eps->V,eps->nini));
197: }
198: for (i=eps->nini;i<eps->ncv;i++) PetscCall(BVSetRandomColumn(eps->V,i));
199: }
201: while (eps->reason == EPS_CONVERGED_ITERATING) {
203: /* Create multivector of constraints from leading columns of V */
204: PetscCall(PetscObjectComposedDataSetInt((PetscObject)eps->V,slepc_blopex_useconstr,1));
205: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
206: constraints = mv_MultiVectorCreateFromSampleVector(&blopex->ii,eps->nds+eps->nconv,eps->V);
208: /* Create multivector where eigenvectors of this run will be stored */
209: PetscCall(PetscObjectComposedDataSetInt((PetscObject)eps->V,slepc_blopex_useconstr,0));
210: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,eps->nconv+blopex->bs));
211: eigenvectors = mv_MultiVectorCreateFromSampleVector(&blopex->ii,blopex->bs,eps->V);
213: #if defined(PETSC_USE_COMPLEX)
214: info = lobpcg_solve_complex(eigenvectors,blopex,OperatorAMultiVector,
215: eps->isgeneralized?blopex:NULL,eps->isgeneralized?OperatorBMultiVector:NULL,
216: blopex,Precond_FnMultiVector,constraints,
217: blopex->blap_fn,blopex->tol,eps->max_it,0,&its,
218: lambda,lambdahist,blopex->bs,eps->errest+eps->nconv,residhist,blopex->bs);
219: #else
220: info = lobpcg_solve_double(eigenvectors,blopex,OperatorAMultiVector,
221: eps->isgeneralized?blopex:NULL,eps->isgeneralized?OperatorBMultiVector:NULL,
222: blopex,Precond_FnMultiVector,constraints,
223: blopex->blap_fn,blopex->tol,eps->max_it,0,&its,
224: lambda,lambdahist,blopex->bs,eps->errest+eps->nconv,residhist,blopex->bs);
225: #endif
226: PetscCheck(info==0,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"BLOPEX failed with exit code=%d",info);
227: mv_MultiVectorDestroy(constraints);
228: mv_MultiVectorDestroy(eigenvectors);
230: for (j=0;j<blopex->bs;j++) {
231: #if defined(PETSC_USE_COMPLEX)
232: eps->eigr[eps->nconv+j] = PetscCMPLX(lambda[j].real,lambda[j].imag);
233: #else
234: eps->eigr[eps->nconv+j] = lambda[j];
235: #endif
236: }
238: if (eps->numbermonitors>0) {
239: for (i=0;i<its;i++) {
240: nconv = 0;
241: for (j=0;j<blopex->bs;j++) {
242: #if defined(PETSC_USE_COMPLEX)
243: eigr[eps->nconv+j] = PetscCMPLX(lambdahist[j+i*blopex->bs].real,lambdahist[j+i*blopex->bs].imag);
244: #else
245: eigr[eps->nconv+j] = lambdahist[j+i*blopex->bs];
246: #endif
247: errest[eps->nconv+j] = residhist[j+i*blopex->bs];
248: if (residhist[j+i*blopex->bs]<=eps->tol) nconv++;
249: }
250: PetscCall(EPSMonitor(eps,eps->its+i,eps->nconv+nconv,eigr,eps->eigi,errest,eps->nconv+blopex->bs));
251: }
252: }
254: eps->its += its;
255: if (info==-1) {
256: eps->reason = EPS_DIVERGED_ITS;
257: break;
258: } else {
259: for (i=0;i<blopex->bs;i++) {
260: if (sigma != 0.0) eps->eigr[eps->nconv+i] += sigma;
261: }
262: eps->nconv += blopex->bs;
263: if (eps->nconv>=eps->nev) eps->reason = EPS_CONVERGED_TOL;
264: }
265: }
267: PetscCall(PetscFree(lambda));
268: if (eps->numbermonitors>0) PetscCall(PetscFree4(lambdahist,eigr,residhist,errest));
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: static PetscErrorCode EPSBLOPEXSetBlockSize_BLOPEX(EPS eps,PetscInt bs)
273: {
274: EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;
276: PetscFunctionBegin;
277: if (bs==PETSC_DEFAULT || bs==PETSC_DECIDE) {
278: ctx->bs = 0;
279: eps->state = EPS_STATE_INITIAL;
280: } else {
281: PetscCheck(bs>0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block size must be >0");
282: ctx->bs = bs;
283: }
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: /*@
288: EPSBLOPEXSetBlockSize - Sets the block size of the BLOPEX solver.
290: Logically Collective
292: Input Parameters:
293: + eps - the linear eigensolver context
294: - bs - the block size
296: Options Database Key:
297: . -eps_blopex_blocksize bs - sets the block size
299: Level: advanced
301: .seealso: [](ch:eps), `EPSBLOPEX`, `EPSBLOPEXGetBlockSize()`
302: @*/
303: PetscErrorCode EPSBLOPEXSetBlockSize(EPS eps,PetscInt bs)
304: {
305: PetscFunctionBegin;
308: PetscTryMethod(eps,"EPSBLOPEXSetBlockSize_C",(EPS,PetscInt),(eps,bs));
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: static PetscErrorCode EPSBLOPEXGetBlockSize_BLOPEX(EPS eps,PetscInt *bs)
313: {
314: EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;
316: PetscFunctionBegin;
317: *bs = ctx->bs;
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: /*@
322: EPSBLOPEXGetBlockSize - Gets the block size used in the BLOPEX solver.
324: Not Collective
326: Input Parameter:
327: . eps - the linear eigensolver context
329: Output Parameter:
330: . bs - the block size
332: Level: advanced
334: .seealso: [](ch:eps), `EPSBLOPEX`, `EPSBLOPEXSetBlockSize()`
335: @*/
336: PetscErrorCode EPSBLOPEXGetBlockSize(EPS eps,PetscInt *bs)
337: {
338: PetscFunctionBegin;
340: PetscAssertPointer(bs,2);
341: PetscUseMethod(eps,"EPSBLOPEXGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: static PetscErrorCode EPSReset_BLOPEX(EPS eps)
346: {
347: EPS_BLOPEX *blopex = (EPS_BLOPEX*)eps->data;
349: PetscFunctionBegin;
350: PetscCall(VecDestroy(&blopex->w));
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: static PetscErrorCode EPSDestroy_BLOPEX(EPS eps)
355: {
356: PetscFunctionBegin;
357: PetscCall(PetscFree(eps->data));
358: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXSetBlockSize_C",NULL));
359: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXGetBlockSize_C",NULL));
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
363: static PetscErrorCode EPSView_BLOPEX(EPS eps,PetscViewer viewer)
364: {
365: EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;
366: PetscBool isascii;
368: PetscFunctionBegin;
369: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
370: if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer," block size %" PetscInt_FMT "\n",ctx->bs));
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: static PetscErrorCode EPSSetFromOptions_BLOPEX(EPS eps,PetscOptionItems PetscOptionsObject)
375: {
376: PetscBool flg;
377: PetscInt bs;
379: PetscFunctionBegin;
380: PetscOptionsHeadBegin(PetscOptionsObject,"EPS BLOPEX Options");
382: PetscCall(PetscOptionsInt("-eps_blopex_blocksize","Block size","EPSBLOPEXSetBlockSize",20,&bs,&flg));
383: if (flg) PetscCall(EPSBLOPEXSetBlockSize(eps,bs));
385: PetscOptionsHeadEnd();
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: /*MC
390: EPSBLOPEX - EPSBLOPEX = "blopex" - A wrapper to BLOPEX {cite:p}`Kny07`.
392: Note:
393: In principle, the method used here is the same as in `EPSLOBPCG`, but in
394: practice convergence may be different.
396: Level: beginner
398: .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`
399: M*/
400: SLEPC_EXTERN PetscErrorCode EPSCreate_BLOPEX(EPS eps)
401: {
402: EPS_BLOPEX *ctx;
404: PetscFunctionBegin;
405: PetscCall(PetscNew(&ctx));
406: eps->data = (void*)ctx;
408: eps->categ = EPS_CATEGORY_PRECOND;
410: eps->ops->solve = EPSSolve_BLOPEX;
411: eps->ops->setup = EPSSetUp_BLOPEX;
412: eps->ops->setupsort = EPSSetUpSort_Basic;
413: eps->ops->setfromoptions = EPSSetFromOptions_BLOPEX;
414: eps->ops->destroy = EPSDestroy_BLOPEX;
415: eps->ops->reset = EPSReset_BLOPEX;
416: eps->ops->view = EPSView_BLOPEX;
417: eps->ops->backtransform = EPSBackTransform_Default;
418: eps->ops->setdefaultst = EPSSetDefaultST_GMRES;
420: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXSetBlockSize_C",EPSBLOPEXSetBlockSize_BLOPEX));
421: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXGetBlockSize_C",EPSBLOPEXGetBlockSize_BLOPEX));
422: if (slepc_blopex_useconstr < 0) PetscCall(PetscObjectComposedDataRegister(&slepc_blopex_useconstr));
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }