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 : SLEPc eigensolver: "lobpcg"
12 :
13 : Method: Locally Optimal Block Preconditioned Conjugate Gradient
14 :
15 : Algorithm:
16 :
17 : LOBPCG with soft and hard locking. Follows the implementation
18 : in BLOPEX [2].
19 :
20 : References:
21 :
22 : [1] A. V. Knyazev, "Toward the optimal preconditioned eigensolver:
23 : locally optimal block preconditioned conjugate gradient method",
24 : SIAM J. Sci. Comput. 23(2):517-541, 2001.
25 :
26 : [2] A. V. Knyazev et al., "Block Locally Optimal Preconditioned
27 : Eigenvalue Xolvers (BLOPEX) in Hypre and PETSc", SIAM J. Sci.
28 : Comput. 29(5):2224-2239, 2007.
29 : */
30 :
31 : #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
32 :
33 : typedef struct {
34 : PetscInt bs; /* block size */
35 : PetscBool lock; /* soft locking active/inactive */
36 : PetscReal restart; /* restart parameter */
37 : PetscInt guard; /* number of guard vectors */
38 : } EPS_LOBPCG;
39 :
40 29 : static PetscErrorCode EPSSetDimensions_LOBPCG(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
41 : {
42 29 : EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
43 29 : PetscInt k;
44 :
45 29 : PetscFunctionBegin;
46 29 : k = PetscMax(3*ctx->bs,((eps->nev-1)/ctx->bs+3)*ctx->bs);
47 29 : if (*ncv!=PETSC_DETERMINE) { /* ncv set */
48 9 : PetscCheck(*ncv>=k,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv is not sufficiently large");
49 20 : } else *ncv = k;
50 29 : if (*mpd==PETSC_DETERMINE) *mpd = 3*ctx->bs;
51 6 : else PetscCheck(*mpd==3*ctx->bs,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"This solver does not allow a value of mpd different from 3*blocksize");
52 29 : PetscFunctionReturn(PETSC_SUCCESS);
53 : }
54 :
55 29 : static PetscErrorCode EPSSetUp_LOBPCG(EPS eps)
56 : {
57 29 : EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
58 :
59 29 : PetscFunctionBegin;
60 29 : EPSCheckHermitianDefinite(eps);
61 29 : EPSCheckNotStructured(eps);
62 29 : if (!ctx->bs) ctx->bs = PetscMin(16,eps->nev);
63 29 : PetscCheck(eps->n-eps->nds>=5*ctx->bs,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The problem size is too small relative to the block size");
64 29 : PetscCall(EPSSetDimensions_LOBPCG(eps,eps->nev,&eps->ncv,&eps->mpd));
65 29 : if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
66 29 : if (!eps->which) eps->which = EPS_SMALLEST_REAL;
67 29 : PetscCheck(eps->which==EPS_SMALLEST_REAL || eps->which==EPS_LARGEST_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only smallest real or largest real eigenvalues");
68 29 : EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION);
69 29 : EPSCheckIgnored(eps,EPS_FEATURE_BALANCE);
70 :
71 29 : if (!ctx->restart) ctx->restart = 0.9;
72 :
73 : /* number of guard vectors */
74 29 : if (ctx->bs==1) ctx->guard = 0;
75 25 : else ctx->guard = PetscMin((PetscInt)((1.0-ctx->restart)*ctx->bs+0.45),ctx->bs-1);
76 :
77 29 : PetscCall(EPSAllocateSolution(eps,0));
78 29 : PetscCall(EPS_SetInnerProduct(eps));
79 29 : PetscCall(DSSetType(eps->ds,DSGHEP));
80 29 : PetscCall(DSAllocate(eps->ds,eps->mpd));
81 29 : PetscCall(EPSSetWorkVecs(eps,1));
82 29 : PetscFunctionReturn(PETSC_SUCCESS);
83 : }
84 :
85 29 : static PetscErrorCode EPSSolve_LOBPCG(EPS eps)
86 : {
87 29 : EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
88 29 : PetscInt i,j,k,nv,ini,nmat,nc,nconv,locked,its,prev=0;
89 29 : PetscReal norm;
90 29 : PetscScalar *eigr,dot;
91 29 : PetscBool breakdown,countc,flip=PETSC_FALSE,checkprecond=PETSC_FALSE;
92 29 : Mat A,B,M,V=NULL,W=NULL;
93 29 : Vec v,z,w=eps->work[0];
94 29 : BV X,Y=NULL,Z,R,P,AX,BX;
95 29 : SlepcSC sc;
96 :
97 29 : PetscFunctionBegin;
98 29 : PetscCall(STGetNumMatrices(eps->st,&nmat));
99 29 : PetscCall(STGetMatrix(eps->st,0,&A));
100 29 : if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
101 21 : else B = NULL;
102 :
103 29 : if (eps->which==EPS_LARGEST_REAL) { /* flip spectrum */
104 2 : flip = PETSC_TRUE;
105 2 : PetscCall(DSGetSlepcSC(eps->ds,&sc));
106 2 : sc->comparison = SlepcCompareSmallestReal;
107 : }
108 :
109 : /* undocumented option to check for a positive-definite preconditioner (turn-off by default) */
110 29 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-eps_lobpcg_checkprecond",&checkprecond,NULL));
111 :
112 : /* 1. Allocate memory */
113 29 : PetscCall(PetscCalloc1(3*ctx->bs,&eigr));
114 29 : PetscCall(BVDuplicateResize(eps->V,3*ctx->bs,&Z));
115 29 : PetscCall(BVDuplicateResize(eps->V,ctx->bs,&X));
116 29 : PetscCall(BVDuplicateResize(eps->V,ctx->bs,&R));
117 29 : PetscCall(BVDuplicateResize(eps->V,ctx->bs,&P));
118 29 : PetscCall(BVDuplicateResize(eps->V,ctx->bs,&AX));
119 29 : if (B) PetscCall(BVDuplicateResize(eps->V,ctx->bs,&BX));
120 29 : nc = eps->nds;
121 29 : if (nc>0 || eps->nev>ctx->bs-ctx->guard) PetscCall(BVDuplicateResize(eps->V,nc+eps->nev,&Y));
122 29 : if (nc>0) {
123 5 : for (j=0;j<nc;j++) {
124 3 : PetscCall(BVGetColumn(eps->V,-nc+j,&v));
125 3 : PetscCall(BVInsertVec(Y,j,v));
126 3 : PetscCall(BVRestoreColumn(eps->V,-nc+j,&v));
127 : }
128 2 : PetscCall(BVSetActiveColumns(Y,0,nc));
129 : }
130 :
131 : /* 2. Apply the constraints to the initial vectors */
132 : /* 3. B-orthogonalize initial vectors */
133 289 : for (k=eps->nini;k<eps->ncv-ctx->bs;k++) { /* Generate more initial vectors if necessary */
134 260 : PetscCall(BVSetRandomColumn(eps->V,k));
135 260 : PetscCall(BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL));
136 : }
137 29 : nv = ctx->bs;
138 29 : PetscCall(BVSetActiveColumns(eps->V,0,nv));
139 29 : PetscCall(BVSetActiveColumns(Z,0,nv));
140 29 : PetscCall(BVCopy(eps->V,Z));
141 29 : PetscCall(BVCopy(Z,X));
142 :
143 : /* 4. Compute initial Ritz vectors */
144 29 : PetscCall(BVMatMult(X,A,AX));
145 29 : PetscCall(DSSetDimensions(eps->ds,nv,0,0));
146 29 : PetscCall(DSGetMat(eps->ds,DS_MAT_A,&M));
147 29 : PetscCall(BVMatProject(AX,NULL,X,M));
148 29 : if (flip) PetscCall(MatScale(M,-1.0));
149 29 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&M));
150 29 : PetscCall(DSSetIdentity(eps->ds,DS_MAT_B));
151 29 : PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
152 29 : PetscCall(DSSolve(eps->ds,eigr,NULL));
153 29 : PetscCall(DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL));
154 29 : PetscCall(DSSynchronize(eps->ds,eigr,NULL));
155 138 : for (j=0;j<nv;j++) eps->eigr[j] = flip? -eigr[j]: eigr[j];
156 29 : PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
157 29 : PetscCall(DSGetMat(eps->ds,DS_MAT_X,&M));
158 29 : PetscCall(BVMultInPlace(X,M,0,nv));
159 29 : PetscCall(BVMultInPlace(AX,M,0,nv));
160 29 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&M));
161 :
162 : /* 5. Initialize range of active iterates */
163 : locked = 0; /* hard-locked vectors, the leading locked columns of V are eigenvectors */
164 : nconv = 0; /* number of converged eigenvalues in the current block */
165 : its = 0; /* iterations for the current block */
166 :
167 : /* 6. Main loop */
168 1354 : while (eps->reason == EPS_CONVERGED_ITERATING) {
169 :
170 1354 : if (ctx->lock) {
171 1003 : PetscCall(BVSetActiveColumns(R,nconv,ctx->bs));
172 1003 : PetscCall(BVSetActiveColumns(AX,nconv,ctx->bs));
173 1003 : if (B) PetscCall(BVSetActiveColumns(BX,nconv,ctx->bs));
174 : }
175 :
176 : /* 7. Compute residuals */
177 1354 : ini = (ctx->lock)? nconv: 0;
178 1354 : PetscCall(BVCopy(AX,R));
179 1354 : if (B) PetscCall(BVMatMult(X,B,BX));
180 5337 : for (j=ini;j<ctx->bs;j++) {
181 3983 : PetscCall(BVGetColumn(R,j,&v));
182 3983 : PetscCall(BVGetColumn(B?BX:X,j,&z));
183 3983 : PetscCall(VecAXPY(v,-eps->eigr[locked+j],z));
184 3983 : PetscCall(BVRestoreColumn(R,j,&v));
185 3983 : PetscCall(BVRestoreColumn(B?BX:X,j,&z));
186 : }
187 :
188 : /* 8. Compute residual norms and update index set of active iterates */
189 : k = ini;
190 : countc = PETSC_TRUE;
191 1580 : for (j=ini;j<ctx->bs;j++) {
192 1554 : i = locked+j;
193 1554 : PetscCall(BVGetColumn(R,j,&v));
194 1554 : PetscCall(VecNorm(v,NORM_2,&norm));
195 1554 : PetscCall(BVRestoreColumn(R,j,&v));
196 1554 : PetscCall((*eps->converged)(eps,eps->eigr[i],eps->eigi[i],norm,&eps->errest[i],eps->convergedctx));
197 1554 : if (countc) {
198 1554 : if (eps->errest[i] < eps->tol) k++;
199 : else countc = PETSC_FALSE;
200 : }
201 1554 : if (!countc && !eps->trackall) break;
202 : }
203 1354 : nconv = k;
204 1354 : eps->nconv = locked + nconv;
205 1354 : if (its) PetscCall(EPSMonitor(eps,eps->its+its,eps->nconv,eps->eigr,eps->eigi,eps->errest,locked+ctx->bs));
206 1354 : PetscCall((*eps->stopping)(eps,eps->its+its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
207 1354 : if (eps->reason != EPS_CONVERGED_ITERATING || nconv >= ctx->bs-ctx->guard) {
208 35 : PetscCall(BVSetActiveColumns(eps->V,locked,eps->nconv));
209 35 : PetscCall(BVSetActiveColumns(X,0,nconv));
210 35 : PetscCall(BVCopy(X,eps->V));
211 : }
212 1354 : if (eps->reason != EPS_CONVERGED_ITERATING) {
213 : break;
214 1325 : } else if (nconv >= ctx->bs-ctx->guard) {
215 6 : eps->its += its-1;
216 6 : its = 0;
217 1319 : } else its++;
218 :
219 1325 : if (nconv >= ctx->bs-ctx->guard) { /* force hard locking of vectors and compute new R */
220 :
221 : /* extend constraints */
222 6 : PetscCall(BVSetActiveColumns(Y,nc+locked,nc+locked+nconv));
223 6 : PetscCall(BVCopy(X,Y));
224 6 : PetscCall(BVSetActiveColumns(Y,0,nc+locked+nconv));
225 :
226 : /* shift work BV's */
227 9 : for (j=nconv;j<ctx->bs;j++) {
228 3 : PetscCall(BVCopyColumn(X,j,j-nconv));
229 3 : PetscCall(BVCopyColumn(R,j,j-nconv));
230 3 : PetscCall(BVCopyColumn(P,j,j-nconv));
231 3 : PetscCall(BVCopyColumn(AX,j,j-nconv));
232 3 : if (B) PetscCall(BVCopyColumn(BX,j,j-nconv));
233 : }
234 :
235 : /* set new initial vectors */
236 6 : PetscCall(BVSetActiveColumns(eps->V,locked+ctx->bs,locked+ctx->bs+nconv));
237 6 : PetscCall(BVSetActiveColumns(X,ctx->bs-nconv,ctx->bs));
238 6 : PetscCall(BVCopy(eps->V,X));
239 28 : for (j=ctx->bs-nconv;j<ctx->bs;j++) {
240 22 : PetscCall(BVGetColumn(X,j,&v));
241 22 : PetscCall(BVOrthogonalizeVec(Y,v,NULL,&norm,&breakdown));
242 22 : if (norm>0.0 && !breakdown) PetscCall(VecScale(v,1.0/norm));
243 : else {
244 0 : PetscCall(PetscInfo(eps,"Orthogonalization of initial vector failed\n"));
245 0 : eps->reason = EPS_DIVERGED_BREAKDOWN;
246 0 : goto diverged;
247 : }
248 22 : PetscCall(BVRestoreColumn(X,j,&v));
249 : }
250 6 : locked += nconv;
251 6 : nconv = 0;
252 6 : PetscCall(BVSetActiveColumns(X,nconv,ctx->bs));
253 :
254 : /* B-orthogonalize initial vectors */
255 6 : PetscCall(BVOrthogonalize(X,NULL));
256 6 : PetscCall(BVSetActiveColumns(Z,nconv,ctx->bs));
257 6 : PetscCall(BVSetActiveColumns(AX,nconv,ctx->bs));
258 6 : PetscCall(BVCopy(X,Z));
259 :
260 : /* compute initial Ritz vectors */
261 6 : nv = ctx->bs;
262 6 : PetscCall(BVMatMult(X,A,AX));
263 6 : PetscCall(DSSetDimensions(eps->ds,nv,0,0));
264 6 : PetscCall(DSGetMat(eps->ds,DS_MAT_A,&M));
265 6 : PetscCall(BVMatProject(AX,NULL,X,M));
266 6 : if (flip) PetscCall(MatScale(M,-1.0));
267 6 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&M));
268 6 : PetscCall(DSSetIdentity(eps->ds,DS_MAT_B));
269 6 : PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
270 6 : PetscCall(DSSolve(eps->ds,eigr,NULL));
271 6 : PetscCall(DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL));
272 6 : PetscCall(DSSynchronize(eps->ds,eigr,NULL));
273 31 : for (j=0;j<nv;j++) if (locked+j<eps->ncv) eps->eigr[locked+j] = flip? -eigr[j]: eigr[j];
274 6 : PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
275 6 : PetscCall(DSGetMat(eps->ds,DS_MAT_X,&M));
276 6 : PetscCall(BVMultInPlace(X,M,0,nv));
277 6 : PetscCall(BVMultInPlace(AX,M,0,nv));
278 6 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&M));
279 :
280 6 : continue; /* skip the rest of the iteration */
281 : }
282 :
283 1319 : ini = (ctx->lock)? nconv: 0;
284 1319 : if (ctx->lock) {
285 972 : PetscCall(BVSetActiveColumns(R,nconv,ctx->bs));
286 972 : PetscCall(BVSetActiveColumns(P,nconv,ctx->bs));
287 972 : PetscCall(BVSetActiveColumns(AX,nconv,ctx->bs));
288 972 : if (B) PetscCall(BVSetActiveColumns(BX,nconv,ctx->bs));
289 : }
290 :
291 : /* 9. Apply preconditioner to the residuals */
292 1319 : PetscCall(BVGetMat(R,&V));
293 1319 : if (prev != ctx->bs-ini) {
294 90 : prev = ctx->bs-ini;
295 90 : PetscCall(MatDestroy(&W));
296 90 : PetscCall(MatDuplicate(V,MAT_SHARE_NONZERO_PATTERN,&W));
297 : }
298 1319 : PetscCall(STApplyMat(eps->st,V,W));
299 1319 : if (checkprecond) {
300 0 : for (j=ini;j<ctx->bs;j++) {
301 0 : PetscCall(MatDenseGetColumnVecRead(V,j-ini,&v));
302 0 : PetscCall(MatDenseGetColumnVecRead(W,j-ini,&w));
303 0 : PetscCall(VecDot(v,w,&dot));
304 0 : PetscCall(MatDenseRestoreColumnVecRead(W,j-ini,&w));
305 0 : PetscCall(MatDenseRestoreColumnVecRead(V,j-ini,&v));
306 0 : if (PetscRealPart(dot)<0.0) {
307 0 : PetscCall(PetscInfo(eps,"The preconditioner is not positive-definite\n"));
308 0 : eps->reason = EPS_DIVERGED_BREAKDOWN;
309 0 : goto diverged;
310 : }
311 : }
312 : }
313 1319 : if (nc+locked>0) {
314 1381 : for (j=ini;j<ctx->bs;j++) {
315 1060 : PetscCall(MatDenseGetColumnVecWrite(W,j-ini,&w));
316 1060 : PetscCall(BVOrthogonalizeVec(Y,w,NULL,&norm,&breakdown));
317 1060 : if (norm>0.0 && !breakdown) PetscCall(VecScale(w,1.0/norm));
318 1060 : PetscCall(MatDenseRestoreColumnVecWrite(W,j-ini,&w));
319 1060 : if (norm<=0.0 || breakdown) {
320 0 : PetscCall(PetscInfo(eps,"Orthogonalization of preconditioned residual failed\n"));
321 0 : eps->reason = EPS_DIVERGED_BREAKDOWN;
322 0 : goto diverged;
323 : }
324 : }
325 : }
326 1319 : PetscCall(MatCopy(W,V,SAME_NONZERO_PATTERN));
327 1319 : PetscCall(BVRestoreMat(R,&V));
328 :
329 : /* 11. B-orthonormalize preconditioned residuals */
330 1319 : PetscCall(BVOrthogonalize(R,NULL));
331 :
332 : /* 13-16. B-orthonormalize conjugate directions */
333 1319 : if (its>1) PetscCall(BVOrthogonalize(P,NULL));
334 :
335 : /* 17-23. Compute symmetric Gram matrices */
336 1319 : PetscCall(BVSetActiveColumns(Z,0,ctx->bs));
337 1319 : PetscCall(BVSetActiveColumns(X,0,ctx->bs));
338 1319 : PetscCall(BVCopy(X,Z));
339 1319 : PetscCall(BVSetActiveColumns(Z,ctx->bs,2*ctx->bs-ini));
340 1319 : PetscCall(BVCopy(R,Z));
341 1319 : if (its>1) {
342 1285 : PetscCall(BVSetActiveColumns(Z,2*ctx->bs-ini,3*ctx->bs-2*ini));
343 1285 : PetscCall(BVCopy(P,Z));
344 : }
345 :
346 1319 : if (its>1) nv = 3*ctx->bs-2*ini;
347 34 : else nv = 2*ctx->bs-ini;
348 :
349 1319 : PetscCall(BVSetActiveColumns(Z,0,nv));
350 1319 : PetscCall(DSSetDimensions(eps->ds,nv,0,0));
351 1319 : PetscCall(DSGetMat(eps->ds,DS_MAT_A,&M));
352 1319 : PetscCall(BVMatProject(Z,A,Z,M));
353 1319 : if (flip) PetscCall(MatScale(M,-1.0));
354 1319 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&M));
355 1319 : PetscCall(DSGetMat(eps->ds,DS_MAT_B,&M));
356 1319 : PetscCall(BVMatProject(Z,B,Z,M)); /* covers also the case B=NULL */
357 1319 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_B,&M));
358 :
359 : /* 24. Solve the generalized eigenvalue problem */
360 1319 : PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
361 1319 : PetscCall(DSSolve(eps->ds,eigr,NULL));
362 1319 : PetscCall(DSSort(eps->ds,eigr,NULL,NULL,NULL,NULL));
363 1319 : PetscCall(DSSynchronize(eps->ds,eigr,NULL));
364 13472 : for (j=0;j<nv;j++) if (locked+j<eps->ncv) eps->eigr[locked+j] = flip? -eigr[j]: eigr[j];
365 1319 : PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
366 :
367 : /* 25-33. Compute Ritz vectors */
368 1319 : PetscCall(DSGetMat(eps->ds,DS_MAT_X,&M));
369 1319 : PetscCall(BVSetActiveColumns(Z,ctx->bs,nv));
370 1319 : if (ctx->lock) PetscCall(BVSetActiveColumns(P,0,ctx->bs));
371 1319 : PetscCall(BVMult(P,1.0,0.0,Z,M));
372 1319 : PetscCall(BVCopy(P,X));
373 1319 : if (ctx->lock) PetscCall(BVSetActiveColumns(P,nconv,ctx->bs));
374 1319 : PetscCall(BVSetActiveColumns(Z,0,ctx->bs));
375 1319 : PetscCall(BVMult(X,1.0,1.0,Z,M));
376 1319 : if (ctx->lock) PetscCall(BVSetActiveColumns(X,nconv,ctx->bs));
377 1319 : PetscCall(BVMatMult(X,A,AX));
378 2673 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&M));
379 : }
380 :
381 29 : diverged:
382 29 : eps->its += its;
383 :
384 29 : if (flip) sc->comparison = SlepcCompareLargestReal;
385 29 : PetscCall(PetscFree(eigr));
386 29 : PetscCall(MatDestroy(&W));
387 29 : if (V) PetscCall(BVRestoreMat(R,&V)); /* only needed when goto diverged is reached */
388 29 : PetscCall(BVDestroy(&Z));
389 29 : PetscCall(BVDestroy(&X));
390 29 : PetscCall(BVDestroy(&R));
391 29 : PetscCall(BVDestroy(&P));
392 29 : PetscCall(BVDestroy(&AX));
393 29 : if (B) PetscCall(BVDestroy(&BX));
394 29 : if (nc>0 || eps->nev>ctx->bs-ctx->guard) PetscCall(BVDestroy(&Y));
395 29 : PetscFunctionReturn(PETSC_SUCCESS);
396 : }
397 :
398 8 : static PetscErrorCode EPSLOBPCGSetBlockSize_LOBPCG(EPS eps,PetscInt bs)
399 : {
400 8 : EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
401 :
402 8 : PetscFunctionBegin;
403 8 : if (bs == PETSC_DEFAULT || bs == PETSC_DECIDE) bs = 0;
404 8 : else PetscCheck(bs>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size %" PetscInt_FMT,bs);
405 8 : if (ctx->bs != bs) {
406 8 : ctx->bs = bs;
407 8 : eps->state = EPS_STATE_INITIAL;
408 : }
409 8 : PetscFunctionReturn(PETSC_SUCCESS);
410 : }
411 :
412 : /*@
413 : EPSLOBPCGSetBlockSize - Sets the block size of the LOBPCG method.
414 :
415 : Logically Collective
416 :
417 : Input Parameters:
418 : + eps - the eigenproblem solver context
419 : - bs - the block size
420 :
421 : Options Database Key:
422 : . -eps_lobpcg_blocksize - Sets the block size
423 :
424 : Level: advanced
425 :
426 : .seealso: EPSLOBPCGGetBlockSize()
427 : @*/
428 8 : PetscErrorCode EPSLOBPCGSetBlockSize(EPS eps,PetscInt bs)
429 : {
430 8 : PetscFunctionBegin;
431 8 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
432 24 : PetscValidLogicalCollectiveInt(eps,bs,2);
433 8 : PetscTryMethod(eps,"EPSLOBPCGSetBlockSize_C",(EPS,PetscInt),(eps,bs));
434 8 : PetscFunctionReturn(PETSC_SUCCESS);
435 : }
436 :
437 1 : static PetscErrorCode EPSLOBPCGGetBlockSize_LOBPCG(EPS eps,PetscInt *bs)
438 : {
439 1 : EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
440 :
441 1 : PetscFunctionBegin;
442 1 : *bs = ctx->bs;
443 1 : PetscFunctionReturn(PETSC_SUCCESS);
444 : }
445 :
446 : /*@
447 : EPSLOBPCGGetBlockSize - Gets the block size used in the LOBPCG method.
448 :
449 : Not Collective
450 :
451 : Input Parameter:
452 : . eps - the eigenproblem solver context
453 :
454 : Output Parameter:
455 : . bs - the block size
456 :
457 : Level: advanced
458 :
459 : .seealso: EPSLOBPCGSetBlockSize()
460 : @*/
461 1 : PetscErrorCode EPSLOBPCGGetBlockSize(EPS eps,PetscInt *bs)
462 : {
463 1 : PetscFunctionBegin;
464 1 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
465 1 : PetscAssertPointer(bs,2);
466 1 : PetscUseMethod(eps,"EPSLOBPCGGetBlockSize_C",(EPS,PetscInt*),(eps,bs));
467 1 : PetscFunctionReturn(PETSC_SUCCESS);
468 : }
469 :
470 1 : static PetscErrorCode EPSLOBPCGSetRestart_LOBPCG(EPS eps,PetscReal restart)
471 : {
472 1 : EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
473 :
474 1 : PetscFunctionBegin;
475 1 : if (restart==(PetscReal)PETSC_DEFAULT || restart==(PetscReal)PETSC_DECIDE) restart = 0.9;
476 1 : PetscCheck(restart>=0.1 && restart<=1.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The restart argument %g must be in the range [0.1,1.0]",(double)restart);
477 1 : if (restart != ctx->restart) {
478 1 : ctx->restart = restart;
479 1 : eps->state = EPS_STATE_INITIAL;
480 : }
481 1 : PetscFunctionReturn(PETSC_SUCCESS);
482 : }
483 :
484 : /*@
485 : EPSLOBPCGSetRestart - Sets the restart parameter for the LOBPCG method.
486 :
487 : Logically Collective
488 :
489 : Input Parameters:
490 : + eps - the eigenproblem solver context
491 : - restart - the percentage of the block of vectors to force a restart
492 :
493 : Options Database Key:
494 : . -eps_lobpcg_restart - Sets the restart parameter
495 :
496 : Notes:
497 : The meaning of this parameter is the proportion of vectors within the
498 : current block iterate that must have converged in order to force a
499 : restart with hard locking.
500 : Allowed values are in the range [0.1,1.0]. The default is 0.9.
501 :
502 : Level: advanced
503 :
504 : .seealso: EPSLOBPCGGetRestart()
505 : @*/
506 1 : PetscErrorCode EPSLOBPCGSetRestart(EPS eps,PetscReal restart)
507 : {
508 1 : PetscFunctionBegin;
509 1 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
510 3 : PetscValidLogicalCollectiveReal(eps,restart,2);
511 1 : PetscTryMethod(eps,"EPSLOBPCGSetRestart_C",(EPS,PetscReal),(eps,restart));
512 1 : PetscFunctionReturn(PETSC_SUCCESS);
513 : }
514 :
515 1 : static PetscErrorCode EPSLOBPCGGetRestart_LOBPCG(EPS eps,PetscReal *restart)
516 : {
517 1 : EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
518 :
519 1 : PetscFunctionBegin;
520 1 : *restart = ctx->restart;
521 1 : PetscFunctionReturn(PETSC_SUCCESS);
522 : }
523 :
524 : /*@
525 : EPSLOBPCGGetRestart - Gets the restart parameter used in the LOBPCG method.
526 :
527 : Not Collective
528 :
529 : Input Parameter:
530 : . eps - the eigenproblem solver context
531 :
532 : Output Parameter:
533 : . restart - the restart parameter
534 :
535 : Level: advanced
536 :
537 : .seealso: EPSLOBPCGSetRestart()
538 : @*/
539 1 : PetscErrorCode EPSLOBPCGGetRestart(EPS eps,PetscReal *restart)
540 : {
541 1 : PetscFunctionBegin;
542 1 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
543 1 : PetscAssertPointer(restart,2);
544 1 : PetscUseMethod(eps,"EPSLOBPCGGetRestart_C",(EPS,PetscReal*),(eps,restart));
545 1 : PetscFunctionReturn(PETSC_SUCCESS);
546 : }
547 :
548 1 : static PetscErrorCode EPSLOBPCGSetLocking_LOBPCG(EPS eps,PetscBool lock)
549 : {
550 1 : EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
551 :
552 1 : PetscFunctionBegin;
553 1 : ctx->lock = lock;
554 1 : PetscFunctionReturn(PETSC_SUCCESS);
555 : }
556 :
557 : /*@
558 : EPSLOBPCGSetLocking - Choose between locking and non-locking variants of
559 : the LOBPCG method.
560 :
561 : Logically Collective
562 :
563 : Input Parameters:
564 : + eps - the eigenproblem solver context
565 : - lock - true if the locking variant must be selected
566 :
567 : Options Database Key:
568 : . -eps_lobpcg_locking - Sets the locking flag
569 :
570 : Notes:
571 : This flag refers to soft locking (converged vectors within the current
572 : block iterate), since hard locking is always used (when nev is larger
573 : than the block size).
574 :
575 : Level: advanced
576 :
577 : .seealso: EPSLOBPCGGetLocking()
578 : @*/
579 1 : PetscErrorCode EPSLOBPCGSetLocking(EPS eps,PetscBool lock)
580 : {
581 1 : PetscFunctionBegin;
582 1 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
583 3 : PetscValidLogicalCollectiveBool(eps,lock,2);
584 1 : PetscTryMethod(eps,"EPSLOBPCGSetLocking_C",(EPS,PetscBool),(eps,lock));
585 1 : PetscFunctionReturn(PETSC_SUCCESS);
586 : }
587 :
588 1 : static PetscErrorCode EPSLOBPCGGetLocking_LOBPCG(EPS eps,PetscBool *lock)
589 : {
590 1 : EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
591 :
592 1 : PetscFunctionBegin;
593 1 : *lock = ctx->lock;
594 1 : PetscFunctionReturn(PETSC_SUCCESS);
595 : }
596 :
597 : /*@
598 : EPSLOBPCGGetLocking - Gets the locking flag used in the LOBPCG method.
599 :
600 : Not Collective
601 :
602 : Input Parameter:
603 : . eps - the eigenproblem solver context
604 :
605 : Output Parameter:
606 : . lock - the locking flag
607 :
608 : Level: advanced
609 :
610 : .seealso: EPSLOBPCGSetLocking()
611 : @*/
612 1 : PetscErrorCode EPSLOBPCGGetLocking(EPS eps,PetscBool *lock)
613 : {
614 1 : PetscFunctionBegin;
615 1 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
616 1 : PetscAssertPointer(lock,2);
617 1 : PetscUseMethod(eps,"EPSLOBPCGGetLocking_C",(EPS,PetscBool*),(eps,lock));
618 1 : PetscFunctionReturn(PETSC_SUCCESS);
619 : }
620 :
621 0 : static PetscErrorCode EPSView_LOBPCG(EPS eps,PetscViewer viewer)
622 : {
623 0 : EPS_LOBPCG *ctx = (EPS_LOBPCG*)eps->data;
624 0 : PetscBool isascii;
625 :
626 0 : PetscFunctionBegin;
627 0 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
628 0 : if (isascii) {
629 0 : PetscCall(PetscViewerASCIIPrintf(viewer," block size %" PetscInt_FMT "\n",ctx->bs));
630 0 : PetscCall(PetscViewerASCIIPrintf(viewer," restart parameter=%g (using %" PetscInt_FMT " guard vectors)\n",(double)ctx->restart,ctx->guard));
631 0 : PetscCall(PetscViewerASCIIPrintf(viewer," soft locking %sactivated\n",ctx->lock?"":"de"));
632 : }
633 0 : PetscFunctionReturn(PETSC_SUCCESS);
634 : }
635 :
636 23 : static PetscErrorCode EPSSetFromOptions_LOBPCG(EPS eps,PetscOptionItems *PetscOptionsObject)
637 : {
638 23 : PetscBool lock,flg;
639 23 : PetscInt bs;
640 23 : PetscReal restart;
641 :
642 23 : PetscFunctionBegin;
643 23 : PetscOptionsHeadBegin(PetscOptionsObject,"EPS LOBPCG Options");
644 :
645 23 : PetscCall(PetscOptionsInt("-eps_lobpcg_blocksize","Block size","EPSLOBPCGSetBlockSize",20,&bs,&flg));
646 23 : if (flg) PetscCall(EPSLOBPCGSetBlockSize(eps,bs));
647 :
648 23 : PetscCall(PetscOptionsReal("-eps_lobpcg_restart","Percentage of the block of vectors to force a restart","EPSLOBPCGSetRestart",0.5,&restart,&flg));
649 23 : if (flg) PetscCall(EPSLOBPCGSetRestart(eps,restart));
650 :
651 23 : PetscCall(PetscOptionsBool("-eps_lobpcg_locking","Choose between locking and non-locking variants","EPSLOBPCGSetLocking",PETSC_TRUE,&lock,&flg));
652 23 : if (flg) PetscCall(EPSLOBPCGSetLocking(eps,lock));
653 :
654 23 : PetscOptionsHeadEnd();
655 23 : PetscFunctionReturn(PETSC_SUCCESS);
656 : }
657 :
658 23 : static PetscErrorCode EPSDestroy_LOBPCG(EPS eps)
659 : {
660 23 : PetscFunctionBegin;
661 23 : PetscCall(PetscFree(eps->data));
662 23 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",NULL));
663 23 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",NULL));
664 23 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",NULL));
665 23 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",NULL));
666 23 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",NULL));
667 23 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",NULL));
668 23 : PetscFunctionReturn(PETSC_SUCCESS);
669 : }
670 :
671 23 : SLEPC_EXTERN PetscErrorCode EPSCreate_LOBPCG(EPS eps)
672 : {
673 23 : EPS_LOBPCG *lobpcg;
674 :
675 23 : PetscFunctionBegin;
676 23 : PetscCall(PetscNew(&lobpcg));
677 23 : eps->data = (void*)lobpcg;
678 23 : lobpcg->lock = PETSC_TRUE;
679 :
680 23 : eps->useds = PETSC_TRUE;
681 23 : eps->categ = EPS_CATEGORY_PRECOND;
682 :
683 23 : eps->ops->solve = EPSSolve_LOBPCG;
684 23 : eps->ops->setup = EPSSetUp_LOBPCG;
685 23 : eps->ops->setupsort = EPSSetUpSort_Default;
686 23 : eps->ops->setfromoptions = EPSSetFromOptions_LOBPCG;
687 23 : eps->ops->destroy = EPSDestroy_LOBPCG;
688 23 : eps->ops->view = EPSView_LOBPCG;
689 23 : eps->ops->backtransform = EPSBackTransform_Default;
690 23 : eps->ops->setdefaultst = EPSSetDefaultST_GMRES;
691 :
692 23 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetBlockSize_C",EPSLOBPCGSetBlockSize_LOBPCG));
693 23 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetBlockSize_C",EPSLOBPCGGetBlockSize_LOBPCG));
694 23 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetRestart_C",EPSLOBPCGSetRestart_LOBPCG));
695 23 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetRestart_C",EPSLOBPCGGetRestart_LOBPCG));
696 23 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGSetLocking_C",EPSLOBPCGSetLocking_LOBPCG));
697 23 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLOBPCGGetLocking_C",EPSLOBPCGGetLocking_LOBPCG));
698 23 : PetscFunctionReturn(PETSC_SUCCESS);
699 : }
|
