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