Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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 : */
13 :
14 : #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
15 : #include "blopex.h"
16 : #include <lobpcg.h>
17 : #include <interpreter.h>
18 : #include <multivector.h>
19 : #include <temp_multivector.h>
20 :
21 : PetscInt slepc_blopex_useconstr = -1;
22 :
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;
31 :
32 783 : static void Precond_FnSingleVector(void *data,void *x,void *y)
33 : {
34 783 : EPS_BLOPEX *blopex = (EPS_BLOPEX*)data;
35 783 : MPI_Comm comm = PetscObjectComm((PetscObject)blopex->st);
36 783 : KSP ksp;
37 :
38 783 : PetscFunctionBegin;
39 783 : PetscCallAbort(comm,STGetKSP(blopex->st,&ksp));
40 783 : PetscCallAbort(comm,KSPSolve(ksp,(Vec)x,(Vec)y));
41 783 : PetscFunctionReturnVoid();
42 : }
43 :
44 259 : static void Precond_FnMultiVector(void *data,void *x,void *y)
45 : {
46 259 : EPS_BLOPEX *blopex = (EPS_BLOPEX*)data;
47 :
48 259 : PetscFunctionBegin;
49 259 : blopex->ii.Eval(Precond_FnSingleVector,data,x,y);
50 259 : PetscFunctionReturnVoid();
51 : }
52 :
53 818 : static void OperatorASingleVector(void *data,void *x,void *y)
54 : {
55 818 : EPS_BLOPEX *blopex = (EPS_BLOPEX*)data;
56 818 : MPI_Comm comm = PetscObjectComm((PetscObject)blopex->st);
57 818 : Mat A,B;
58 818 : PetscScalar sigma;
59 818 : PetscInt nmat;
60 :
61 818 : PetscFunctionBegin;
62 818 : PetscCallAbort(comm,STGetNumMatrices(blopex->st,&nmat));
63 818 : PetscCallAbort(comm,STGetMatrix(blopex->st,0,&A));
64 818 : if (nmat>1) PetscCallAbort(comm,STGetMatrix(blopex->st,1,&B));
65 818 : PetscCallAbort(comm,MatMult(A,(Vec)x,(Vec)y));
66 818 : PetscCallAbort(comm,STGetShift(blopex->st,&sigma));
67 818 : if (sigma != 0.0) {
68 181 : if (nmat>1) PetscCallAbort(comm,MatMult(B,(Vec)x,blopex->w));
69 140 : else PetscCallAbort(comm,VecCopy((Vec)x,blopex->w));
70 181 : PetscCallAbort(comm,VecAXPY((Vec)y,-sigma,blopex->w));
71 : }
72 818 : PetscFunctionReturnVoid();
73 : }
74 :
75 268 : static void OperatorAMultiVector(void *data,void *x,void *y)
76 : {
77 268 : EPS_BLOPEX *blopex = (EPS_BLOPEX*)data;
78 :
79 268 : PetscFunctionBegin;
80 268 : blopex->ii.Eval(OperatorASingleVector,data,x,y);
81 268 : PetscFunctionReturnVoid();
82 : }
83 :
84 463 : static void OperatorBSingleVector(void *data,void *x,void *y)
85 : {
86 463 : EPS_BLOPEX *blopex = (EPS_BLOPEX*)data;
87 463 : MPI_Comm comm = PetscObjectComm((PetscObject)blopex->st);
88 463 : Mat B;
89 :
90 463 : PetscFunctionBegin;
91 463 : PetscCallAbort(comm,STGetMatrix(blopex->st,1,&B));
92 463 : PetscCallAbort(comm,MatMult(B,(Vec)x,(Vec)y));
93 463 : PetscFunctionReturnVoid();
94 : }
95 :
96 156 : static void OperatorBMultiVector(void *data,void *x,void *y)
97 : {
98 156 : EPS_BLOPEX *blopex = (EPS_BLOPEX*)data;
99 :
100 156 : PetscFunctionBegin;
101 156 : blopex->ii.Eval(OperatorBSingleVector,data,x,y);
102 156 : PetscFunctionReturnVoid();
103 : }
104 :
105 7 : static PetscErrorCode EPSSetDimensions_BLOPEX(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
106 : {
107 7 : EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;
108 7 : PetscInt k;
109 :
110 7 : PetscFunctionBegin;
111 7 : if (*nev==0) *nev = 1;
112 7 : k = ((*nev-1)/ctx->bs+1)*ctx->bs;
113 7 : if (*ncv!=PETSC_DETERMINE) { /* ncv set */
114 1 : PetscCheck(*ncv>=k,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv is not sufficiently large");
115 6 : } else *ncv = k;
116 7 : if (*mpd==PETSC_DETERMINE) *mpd = *ncv;
117 1 : else PetscCall(PetscInfo(eps,"Warning: given value of mpd ignored\n"));
118 7 : PetscFunctionReturn(PETSC_SUCCESS);
119 : }
120 :
121 7 : static PetscErrorCode EPSSetUp_BLOPEX(EPS eps)
122 : {
123 7 : EPS_BLOPEX *blopex = (EPS_BLOPEX*)eps->data;
124 7 : PetscBool flg;
125 7 : KSP ksp;
126 :
127 7 : PetscFunctionBegin;
128 7 : EPSCheckHermitianDefinite(eps);
129 7 : EPSCheckNotStructured(eps);
130 11 : if (!blopex->bs) blopex->bs = PetscMin(16,eps->nev?eps->nev:1);
131 7 : PetscCall(EPSSetDimensions_BLOPEX(eps,&eps->nev,&eps->ncv,&eps->mpd));
132 7 : if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
133 7 : if (!eps->which) eps->which = EPS_SMALLEST_REAL;
134 7 : PetscCheck(eps->which==EPS_SMALLEST_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only smallest real eigenvalues");
135 7 : EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_STOPPING);
136 7 : EPSCheckIgnored(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_EXTRACTION);
137 :
138 7 : blopex->st = eps->st;
139 :
140 7 : PetscCheck(eps->converged==EPSConvergedRelative || eps->converged==EPSConvergedAbsolute,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Convergence test not supported in this solver");
141 7 : if (eps->converged == EPSConvergedRelative) {
142 4 : blopex->tol.absolute = 0.0;
143 5 : blopex->tol.relative = SlepcDefaultTol(eps->tol);
144 : } else { /* EPSConvergedAbsolute */
145 3 : blopex->tol.absolute = SlepcDefaultTol(eps->tol);
146 3 : blopex->tol.relative = 0.0;
147 : }
148 :
149 7 : SLEPCSetupInterpreter(&blopex->ii);
150 :
151 7 : PetscCall(STGetKSP(eps->st,&ksp));
152 7 : PetscCall(PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg));
153 7 : if (!flg) PetscCall(PetscInfo(eps,"Warning: ignoring KSP, should use KSPPREONLY\n"));
154 :
155 : /* allocate memory */
156 7 : if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
157 7 : PetscCall(PetscObjectTypeCompareAny((PetscObject)eps->V,&flg,BVVECS,BVCONTIGUOUS,""));
158 7 : if (!flg) { /* blopex only works with BVVECS or BVCONTIGUOUS */
159 6 : PetscCall(BVSetType(eps->V,BVCONTIGUOUS));
160 : }
161 7 : PetscCall(EPSAllocateSolution(eps,0));
162 7 : if (!blopex->w) PetscCall(BVCreateVec(eps->V,&blopex->w));
163 :
164 : #if defined(PETSC_USE_COMPLEX)
165 7 : blopex->blap_fn.zpotrf = PETSC_zpotrf_interface;
166 7 : 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 7 : PetscFunctionReturn(PETSC_SUCCESS);
172 : }
173 :
174 7 : static PetscErrorCode EPSSolve_BLOPEX(EPS eps)
175 : {
176 7 : EPS_BLOPEX *blopex = (EPS_BLOPEX*)eps->data;
177 7 : PetscScalar sigma,*eigr=NULL;
178 7 : PetscReal *errest=NULL;
179 7 : int i,j,info,its,nconv;
180 7 : double *residhist=NULL;
181 7 : mv_MultiVectorPtr eigenvectors,constraints;
182 : #if defined(PETSC_USE_COMPLEX)
183 7 : komplex *lambda=NULL,*lambdahist=NULL;
184 : #else
185 : double *lambda=NULL,*lambdahist=NULL;
186 : #endif
187 :
188 7 : PetscFunctionBegin;
189 7 : PetscCall(STGetShift(eps->st,&sigma));
190 7 : PetscCall(PetscMalloc1(blopex->bs,&lambda));
191 7 : 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 :
193 : /* Complete the initial basis with random vectors */
194 8 : for (i=0;i<eps->nini;i++) { /* in case the initial vectors were also set with VecSetRandom */
195 1 : PetscCall(BVSetRandomColumn(eps->V,eps->nini));
196 : }
197 41 : for (i=eps->nini;i<eps->ncv;i++) PetscCall(BVSetRandomColumn(eps->V,i));
198 :
199 16 : while (eps->reason == EPS_CONVERGED_ITERATING) {
200 :
201 : /* Create multivector of constraints from leading columns of V */
202 9 : PetscCall(PetscObjectComposedDataSetInt((PetscObject)eps->V,slepc_blopex_useconstr,1));
203 9 : PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
204 9 : constraints = mv_MultiVectorCreateFromSampleVector(&blopex->ii,eps->nds+eps->nconv,eps->V);
205 :
206 : /* Create multivector where eigenvectors of this run will be stored */
207 9 : PetscCall(PetscObjectComposedDataSetInt((PetscObject)eps->V,slepc_blopex_useconstr,0));
208 9 : PetscCall(BVSetActiveColumns(eps->V,eps->nconv,eps->nconv+blopex->bs));
209 9 : eigenvectors = mv_MultiVectorCreateFromSampleVector(&blopex->ii,blopex->bs,eps->V);
210 :
211 : #if defined(PETSC_USE_COMPLEX)
212 15 : info = lobpcg_solve_complex(eigenvectors,blopex,OperatorAMultiVector,
213 9 : eps->isgeneralized?blopex:NULL,eps->isgeneralized?OperatorBMultiVector:NULL,
214 : blopex,Precond_FnMultiVector,constraints,
215 : blopex->blap_fn,blopex->tol,eps->max_it,0,&its,
216 9 : lambda,lambdahist,blopex->bs,eps->errest+eps->nconv,residhist,blopex->bs);
217 : #else
218 : info = lobpcg_solve_double(eigenvectors,blopex,OperatorAMultiVector,
219 : eps->isgeneralized?blopex:NULL,eps->isgeneralized?OperatorBMultiVector:NULL,
220 : blopex,Precond_FnMultiVector,constraints,
221 : blopex->blap_fn,blopex->tol,eps->max_it,0,&its,
222 : lambda,lambdahist,blopex->bs,eps->errest+eps->nconv,residhist,blopex->bs);
223 : #endif
224 9 : PetscCheck(info==0,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"BLOPEX failed with exit code=%d",info);
225 9 : mv_MultiVectorDestroy(constraints);
226 9 : mv_MultiVectorDestroy(eigenvectors);
227 :
228 53 : for (j=0;j<blopex->bs;j++) {
229 : #if defined(PETSC_USE_COMPLEX)
230 35 : eps->eigr[eps->nconv+j] = PetscCMPLX(lambda[j].real,lambda[j].imag);
231 : #else
232 : eps->eigr[eps->nconv+j] = lambda[j];
233 : #endif
234 : }
235 :
236 9 : if (eps->numbermonitors>0) {
237 0 : for (i=0;i<its;i++) {
238 : nconv = 0;
239 0 : for (j=0;j<blopex->bs;j++) {
240 : #if defined(PETSC_USE_COMPLEX)
241 0 : eigr[eps->nconv+j] = PetscCMPLX(lambdahist[j+i*blopex->bs].real,lambdahist[j+i*blopex->bs].imag);
242 : #else
243 : eigr[eps->nconv+j] = lambdahist[j+i*blopex->bs];
244 : #endif
245 0 : errest[eps->nconv+j] = residhist[j+i*blopex->bs];
246 0 : if (residhist[j+i*blopex->bs]<=eps->tol) nconv++;
247 : }
248 0 : PetscCall(EPSMonitor(eps,eps->its+i,eps->nconv+nconv,eigr,eps->eigi,errest,eps->nconv+blopex->bs));
249 : }
250 : }
251 :
252 9 : eps->its += its;
253 9 : if (info==-1) {
254 0 : eps->reason = EPS_DIVERGED_ITS;
255 0 : break;
256 : } else {
257 44 : for (i=0;i<blopex->bs;i++) {
258 35 : if (sigma != 0.0) eps->eigr[eps->nconv+i] += sigma;
259 : }
260 9 : eps->nconv += blopex->bs;
261 9 : if (eps->nconv>=eps->nev) eps->reason = EPS_CONVERGED_TOL;
262 : }
263 : }
264 :
265 7 : PetscCall(PetscFree(lambda));
266 7 : if (eps->numbermonitors>0) PetscCall(PetscFree4(lambdahist,eigr,residhist,errest));
267 7 : PetscFunctionReturn(PETSC_SUCCESS);
268 : }
269 :
270 2 : static PetscErrorCode EPSBLOPEXSetBlockSize_BLOPEX(EPS eps,PetscInt bs)
271 : {
272 2 : EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;
273 :
274 2 : PetscFunctionBegin;
275 2 : if (bs==PETSC_DEFAULT || bs==PETSC_DECIDE) {
276 0 : ctx->bs = 0;
277 0 : eps->state = EPS_STATE_INITIAL;
278 : } else {
279 2 : PetscCheck(bs>0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block size must be >0");
280 2 : ctx->bs = bs;
281 : }
282 2 : PetscFunctionReturn(PETSC_SUCCESS);
283 : }
284 :
285 : /*@
286 : EPSBLOPEXSetBlockSize - Sets the block size of the BLOPEX solver.
287 :
288 : Logically Collective
289 :
290 : Input Parameters:
291 : + eps - the eigenproblem solver context
292 : - bs - the block size
293 :
294 : Options Database Key:
295 : . -eps_blopex_blocksize - Sets the block size
296 :
297 : Level: advanced
298 :
299 : .seealso: EPSBLOPEXGetBlockSize()
300 : @*/
301 2 : PetscErrorCode EPSBLOPEXSetBlockSize(EPS eps,PetscInt bs)
302 : {
303 2 : PetscFunctionBegin;
304 2 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
305 6 : PetscValidLogicalCollectiveInt(eps,bs,2);
306 2 : PetscTryMethod(eps,"EPSBLOPEXSetBlockSize_C",(EPS,PetscInt),(eps,bs));
307 2 : PetscFunctionReturn(PETSC_SUCCESS);
308 : }
309 :
310 1 : static PetscErrorCode EPSBLOPEXGetBlockSize_BLOPEX(EPS eps,PetscInt *bs)
311 : {
312 1 : EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;
313 :
314 1 : PetscFunctionBegin;
315 1 : *bs = ctx->bs;
316 1 : PetscFunctionReturn(PETSC_SUCCESS);
317 : }
318 :
319 : /*@
320 : EPSBLOPEXGetBlockSize - Gets the block size used in the BLOPEX solver.
321 :
322 : Not Collective
323 :
324 : Input Parameter:
325 : . eps - the eigenproblem solver context
326 :
327 : Output Parameter:
328 : . bs - the block size
329 :
330 : Level: advanced
331 :
332 : .seealso: EPSBLOPEXSetBlockSize()
333 : @*/
334 1 : PetscErrorCode EPSBLOPEXGetBlockSize(EPS eps,PetscInt *bs)
335 : {
336 1 : PetscFunctionBegin;
337 1 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
338 1 : PetscAssertPointer(bs,2);
339 1 : PetscUseMethod(eps,"EPSBLOPEXGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
340 1 : PetscFunctionReturn(PETSC_SUCCESS);
341 : }
342 :
343 6 : static PetscErrorCode EPSReset_BLOPEX(EPS eps)
344 : {
345 6 : EPS_BLOPEX *blopex = (EPS_BLOPEX*)eps->data;
346 :
347 6 : PetscFunctionBegin;
348 6 : PetscCall(VecDestroy(&blopex->w));
349 6 : PetscFunctionReturn(PETSC_SUCCESS);
350 : }
351 :
352 6 : static PetscErrorCode EPSDestroy_BLOPEX(EPS eps)
353 : {
354 6 : PetscFunctionBegin;
355 6 : LOBPCG_DestroyRandomContext();
356 6 : PetscCall(PetscFree(eps->data));
357 6 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXSetBlockSize_C",NULL));
358 6 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXGetBlockSize_C",NULL));
359 6 : PetscFunctionReturn(PETSC_SUCCESS);
360 : }
361 :
362 1 : static PetscErrorCode EPSView_BLOPEX(EPS eps,PetscViewer viewer)
363 : {
364 1 : EPS_BLOPEX *ctx = (EPS_BLOPEX*)eps->data;
365 1 : PetscBool isascii;
366 :
367 1 : PetscFunctionBegin;
368 1 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
369 1 : if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer," block size %" PetscInt_FMT "\n",ctx->bs));
370 1 : PetscFunctionReturn(PETSC_SUCCESS);
371 : }
372 :
373 6 : static PetscErrorCode EPSSetFromOptions_BLOPEX(EPS eps,PetscOptionItems *PetscOptionsObject)
374 : {
375 6 : PetscBool flg;
376 6 : PetscInt bs;
377 :
378 6 : PetscFunctionBegin;
379 6 : PetscOptionsHeadBegin(PetscOptionsObject,"EPS BLOPEX Options");
380 :
381 6 : PetscCall(PetscOptionsInt("-eps_blopex_blocksize","Block size","EPSBLOPEXSetBlockSize",20,&bs,&flg));
382 6 : if (flg) PetscCall(EPSBLOPEXSetBlockSize(eps,bs));
383 :
384 6 : PetscOptionsHeadEnd();
385 :
386 6 : LOBPCG_SetFromOptionsRandomContext();
387 6 : PetscFunctionReturn(PETSC_SUCCESS);
388 : }
389 :
390 6 : SLEPC_EXTERN PetscErrorCode EPSCreate_BLOPEX(EPS eps)
391 : {
392 6 : EPS_BLOPEX *ctx;
393 :
394 6 : PetscFunctionBegin;
395 6 : PetscCall(PetscNew(&ctx));
396 6 : eps->data = (void*)ctx;
397 :
398 6 : eps->categ = EPS_CATEGORY_PRECOND;
399 :
400 6 : eps->ops->solve = EPSSolve_BLOPEX;
401 6 : eps->ops->setup = EPSSetUp_BLOPEX;
402 6 : eps->ops->setupsort = EPSSetUpSort_Basic;
403 6 : eps->ops->setfromoptions = EPSSetFromOptions_BLOPEX;
404 6 : eps->ops->destroy = EPSDestroy_BLOPEX;
405 6 : eps->ops->reset = EPSReset_BLOPEX;
406 6 : eps->ops->view = EPSView_BLOPEX;
407 6 : eps->ops->backtransform = EPSBackTransform_Default;
408 6 : eps->ops->setdefaultst = EPSSetDefaultST_GMRES;
409 :
410 6 : LOBPCG_InitRandomContext(PetscObjectComm((PetscObject)eps),NULL);
411 6 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXSetBlockSize_C",EPSBLOPEXSetBlockSize_BLOPEX));
412 6 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSBLOPEXGetBlockSize_C",EPSBLOPEXGetBlockSize_BLOPEX));
413 6 : if (slepc_blopex_useconstr < 0) PetscCall(PetscObjectComposedDataRegister(&slepc_blopex_useconstr));
414 6 : PetscFunctionReturn(PETSC_SUCCESS);
415 : }
|
