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