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