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