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: "rqcg"
12 :
13 : Method: Rayleigh Quotient Conjugate Gradient
14 :
15 : Algorithm:
16 :
17 : Conjugate Gradient minimization of the Rayleigh quotient with
18 : periodic Rayleigh-Ritz acceleration.
19 :
20 : References:
21 :
22 : [1] L. Bergamaschi et al., "Parallel preconditioned conjugate gradient
23 : optimization of the Rayleigh quotient for the solution of sparse
24 : eigenproblems", Appl. Math. Comput. 175(2):1694-1715, 2006.
25 : */
26 :
27 : #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
28 :
29 : static PetscErrorCode EPSSolve_RQCG(EPS);
30 :
31 : typedef struct {
32 : PetscInt nrest; /* user-provided reset parameter */
33 : PetscInt allocsize; /* number of columns of work BV's allocated at setup */
34 : BV AV,W,P,G;
35 : } EPS_RQCG;
36 :
37 14 : static PetscErrorCode EPSSetUp_RQCG(EPS eps)
38 : {
39 14 : PetscInt nmat;
40 14 : EPS_RQCG *ctx = (EPS_RQCG*)eps->data;
41 :
42 14 : PetscFunctionBegin;
43 14 : EPSCheckHermitianDefinite(eps);
44 14 : PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
45 14 : if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
46 14 : if (!eps->which) eps->which = EPS_SMALLEST_REAL;
47 14 : PetscCheck(eps->which==EPS_SMALLEST_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only smallest real eigenvalues");
48 14 : EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION);
49 14 : EPSCheckIgnored(eps,EPS_FEATURE_BALANCE);
50 :
51 14 : if (!ctx->nrest) ctx->nrest = 20;
52 :
53 14 : PetscCall(EPSAllocateSolution(eps,0));
54 14 : PetscCall(EPS_SetInnerProduct(eps));
55 :
56 14 : PetscCall(STGetNumMatrices(eps->st,&nmat));
57 14 : if (!ctx->allocsize) {
58 11 : ctx->allocsize = eps->mpd;
59 11 : PetscCall(BVDuplicateResize(eps->V,eps->mpd,&ctx->AV));
60 11 : if (nmat>1) PetscCall(BVDuplicate(ctx->AV,&ctx->W));
61 11 : PetscCall(BVDuplicate(ctx->AV,&ctx->P));
62 11 : PetscCall(BVDuplicate(ctx->AV,&ctx->G));
63 3 : } else if (ctx->allocsize!=eps->mpd) {
64 1 : ctx->allocsize = eps->mpd;
65 1 : PetscCall(BVResize(ctx->AV,eps->mpd,PETSC_FALSE));
66 1 : if (nmat>1) PetscCall(BVResize(ctx->W,eps->mpd,PETSC_FALSE));
67 1 : PetscCall(BVResize(ctx->P,eps->mpd,PETSC_FALSE));
68 1 : PetscCall(BVResize(ctx->G,eps->mpd,PETSC_FALSE));
69 : }
70 14 : PetscCall(DSSetType(eps->ds,DSHEP));
71 14 : PetscCall(DSAllocate(eps->ds,eps->ncv));
72 14 : PetscCall(EPSSetWorkVecs(eps,1));
73 14 : PetscFunctionReturn(PETSC_SUCCESS);
74 : }
75 :
76 14 : static PetscErrorCode EPSSolve_RQCG(EPS eps)
77 : {
78 14 : EPS_RQCG *ctx = (EPS_RQCG*)eps->data;
79 14 : PetscInt i,j,k,ld,nv,ncv = eps->ncv,kini,nmat;
80 14 : PetscScalar *C,*gamma,g,pap,pbp,pbx,pax,nu,mu,alpha,beta;
81 14 : PetscReal resnorm,a,b,c,d,disc,t;
82 14 : PetscBool reset;
83 14 : Mat A,B,Q,Q1;
84 14 : Vec v,av,bv,p,w=eps->work[0];
85 :
86 14 : PetscFunctionBegin;
87 14 : PetscCall(DSGetLeadingDimension(eps->ds,&ld));
88 14 : PetscCall(STGetNumMatrices(eps->st,&nmat));
89 14 : PetscCall(STGetMatrix(eps->st,0,&A));
90 14 : if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
91 12 : else B = NULL;
92 14 : PetscCall(PetscMalloc1(eps->mpd,&gamma));
93 :
94 14 : kini = eps->nini;
95 593 : while (eps->reason == EPS_CONVERGED_ITERATING) {
96 579 : eps->its++;
97 579 : nv = PetscMin(eps->nconv+eps->mpd,ncv);
98 579 : PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,0));
99 857 : for (;kini<nv;kini++) { /* Generate more initial vectors if necessary */
100 278 : PetscCall(BVSetRandomColumn(eps->V,kini));
101 278 : PetscCall(BVOrthonormalizeColumn(eps->V,kini,PETSC_TRUE,NULL,NULL));
102 : }
103 579 : reset = (eps->its>1 && (eps->its-1)%ctx->nrest==0)? PETSC_TRUE: PETSC_FALSE;
104 :
105 579 : if (reset) {
106 : /* Prevent BVDotVec below to use B-product, restored at the end */
107 36 : PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
108 :
109 : /* Compute Rayleigh quotient */
110 36 : PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
111 36 : PetscCall(BVSetActiveColumns(ctx->AV,0,nv-eps->nconv));
112 36 : PetscCall(BVMatMult(eps->V,A,ctx->AV));
113 36 : PetscCall(DSGetArray(eps->ds,DS_MAT_A,&C));
114 804 : for (i=eps->nconv;i<nv;i++) {
115 768 : PetscCall(BVSetActiveColumns(eps->V,eps->nconv,i+1));
116 768 : PetscCall(BVGetColumn(ctx->AV,i-eps->nconv,&av));
117 768 : PetscCall(BVDotVec(eps->V,av,C+eps->nconv+i*ld));
118 768 : PetscCall(BVRestoreColumn(ctx->AV,i-eps->nconv,&av));
119 8962 : for (j=eps->nconv;j<i-1;j++) C[i+j*ld] = PetscConj(C[j+i*ld]);
120 : }
121 36 : PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&C));
122 36 : PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
123 :
124 : /* Solve projected problem */
125 36 : PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
126 36 : PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
127 36 : PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
128 :
129 : /* Update vectors V(:,idx) = V * Y(:,idx) */
130 36 : PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
131 36 : PetscCall(BVMultInPlace(eps->V,Q,eps->nconv,nv));
132 36 : PetscCall(MatDenseGetSubMatrix(Q,eps->nconv,PETSC_DECIDE,eps->nconv,PETSC_DECIDE,&Q1));
133 36 : PetscCall(BVMultInPlace(ctx->AV,Q1,0,nv-eps->nconv));
134 36 : PetscCall(MatDenseRestoreSubMatrix(Q,&Q1));
135 36 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
136 36 : if (B) PetscCall(BVSetMatrix(eps->V,B,PETSC_FALSE));
137 : } else {
138 : /* No need to do Rayleigh-Ritz, just take diag(V'*A*V) */
139 9938 : for (i=eps->nconv;i<nv;i++) {
140 9395 : PetscCall(BVGetColumn(eps->V,i,&v));
141 9395 : PetscCall(BVGetColumn(ctx->AV,i-eps->nconv,&av));
142 9395 : PetscCall(MatMult(A,v,av));
143 9395 : PetscCall(VecDot(av,v,eps->eigr+i));
144 9395 : PetscCall(BVRestoreColumn(eps->V,i,&v));
145 9395 : PetscCall(BVRestoreColumn(ctx->AV,i-eps->nconv,&av));
146 : }
147 : }
148 :
149 : /* Compute gradient and check convergence */
150 579 : k = -1;
151 10742 : for (i=eps->nconv;i<nv;i++) {
152 10163 : PetscCall(BVGetColumn(eps->V,i,&v));
153 10163 : PetscCall(BVGetColumn(ctx->AV,i-eps->nconv,&av));
154 10163 : PetscCall(BVGetColumn(ctx->G,i-eps->nconv,&p));
155 10163 : if (B) {
156 1415 : PetscCall(BVGetColumn(ctx->W,i-eps->nconv,&bv));
157 1415 : PetscCall(MatMult(B,v,bv));
158 1415 : PetscCall(VecWAXPY(p,-eps->eigr[i],bv,av));
159 1415 : PetscCall(BVRestoreColumn(ctx->W,i-eps->nconv,&bv));
160 8748 : } else PetscCall(VecWAXPY(p,-eps->eigr[i],v,av));
161 10163 : PetscCall(BVRestoreColumn(eps->V,i,&v));
162 10163 : PetscCall(BVRestoreColumn(ctx->AV,i-eps->nconv,&av));
163 10163 : PetscCall(VecNorm(p,NORM_2,&resnorm));
164 10163 : PetscCall(BVRestoreColumn(ctx->G,i-eps->nconv,&p));
165 10163 : PetscCall((*eps->converged)(eps,eps->eigr[i],0.0,resnorm,&eps->errest[i],eps->convergedctx));
166 10163 : if (k==-1 && eps->errest[i] >= eps->tol) k = i;
167 : }
168 579 : if (k==-1) k = nv;
169 579 : PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
170 :
171 : /* The next lines are necessary to avoid DS zeroing eigr */
172 579 : PetscCall(DSGetArray(eps->ds,DS_MAT_A,&C));
173 719 : for (i=eps->nconv;i<k;i++) C[i+i*ld] = eps->eigr[i];
174 579 : PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&C));
175 :
176 579 : if (eps->reason == EPS_CONVERGED_ITERATING) {
177 :
178 : /* Search direction */
179 10471 : for (i=0;i<nv-eps->nconv;i++) {
180 9906 : PetscCall(BVGetColumn(ctx->G,i,&v));
181 9906 : PetscCall(STApply(eps->st,v,w));
182 9906 : PetscCall(VecDot(w,v,&g));
183 9906 : PetscCall(BVRestoreColumn(ctx->G,i,&v));
184 9906 : beta = (!reset && eps->its>1)? g/gamma[i]: 0.0;
185 9906 : gamma[i] = g;
186 9906 : PetscCall(BVGetColumn(ctx->P,i,&v));
187 9906 : PetscCall(VecAXPBY(v,1.0,beta,w));
188 9906 : if (i+eps->nconv>0) {
189 9690 : PetscCall(BVSetActiveColumns(eps->V,0,i+eps->nconv));
190 9690 : PetscCall(BVOrthogonalizeVec(eps->V,v,NULL,NULL,NULL));
191 : }
192 9906 : PetscCall(BVRestoreColumn(ctx->P,i,&v));
193 : }
194 :
195 : /* Minimization problem */
196 10471 : for (i=eps->nconv;i<nv;i++) {
197 9906 : PetscCall(BVGetColumn(eps->V,i,&v));
198 9906 : PetscCall(BVGetColumn(ctx->AV,i-eps->nconv,&av));
199 9906 : PetscCall(BVGetColumn(ctx->P,i-eps->nconv,&p));
200 9906 : PetscCall(VecDot(av,v,&nu));
201 9906 : PetscCall(VecDot(av,p,&pax));
202 9906 : PetscCall(MatMult(A,p,w));
203 9906 : PetscCall(VecDot(w,p,&pap));
204 9906 : if (B) {
205 1381 : PetscCall(BVGetColumn(ctx->W,i-eps->nconv,&bv));
206 1381 : PetscCall(VecDot(bv,v,&mu));
207 1381 : PetscCall(VecDot(bv,p,&pbx));
208 1381 : PetscCall(BVRestoreColumn(ctx->W,i-eps->nconv,&bv));
209 1381 : PetscCall(MatMult(B,p,w));
210 1381 : PetscCall(VecDot(w,p,&pbp));
211 : } else {
212 8525 : PetscCall(VecDot(v,v,&mu));
213 8525 : PetscCall(VecDot(v,p,&pbx));
214 8525 : PetscCall(VecDot(p,p,&pbp));
215 : }
216 9906 : PetscCall(BVRestoreColumn(ctx->AV,i-eps->nconv,&av));
217 9906 : a = PetscRealPart(pap*pbx-pax*pbp);
218 9906 : b = PetscRealPart(nu*pbp-mu*pap);
219 9906 : c = PetscRealPart(mu*pax-nu*pbx);
220 10493 : t = PetscMax(PetscMax(PetscAbsReal(a),PetscAbsReal(b)),PetscAbsReal(c));
221 9906 : if (t!=0.0) { a /= t; b /= t; c /= t; }
222 9906 : disc = b*b-4.0*a*c;
223 9906 : d = PetscSqrtReal(PetscAbsReal(disc));
224 9906 : if (b>=0.0 && a!=0.0) alpha = (b+d)/(2.0*a);
225 9066 : else if (b!=d) alpha = 2.0*c/(b-d);
226 : else alpha = 0;
227 : /* Next iterate */
228 9906 : if (alpha!=0.0) PetscCall(VecAXPY(v,alpha,p));
229 9906 : PetscCall(BVRestoreColumn(eps->V,i,&v));
230 9906 : PetscCall(BVRestoreColumn(ctx->P,i-eps->nconv,&p));
231 9906 : PetscCall(BVOrthonormalizeColumn(eps->V,i,PETSC_TRUE,NULL,NULL));
232 : }
233 : }
234 :
235 579 : PetscCall(EPSMonitor(eps,eps->its,k,eps->eigr,eps->eigi,eps->errest,nv));
236 579 : eps->nconv = k;
237 : }
238 :
239 14 : PetscCall(PetscFree(gamma));
240 14 : PetscFunctionReturn(PETSC_SUCCESS);
241 : }
242 :
243 2 : static PetscErrorCode EPSRQCGSetReset_RQCG(EPS eps,PetscInt nrest)
244 : {
245 2 : EPS_RQCG *ctx = (EPS_RQCG*)eps->data;
246 :
247 2 : PetscFunctionBegin;
248 2 : if (nrest==PETSC_DEFAULT) {
249 0 : ctx->nrest = 0;
250 0 : eps->state = EPS_STATE_INITIAL;
251 : } else {
252 2 : PetscCheck(nrest>0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reset parameter must be >0");
253 2 : ctx->nrest = nrest;
254 : }
255 2 : PetscFunctionReturn(PETSC_SUCCESS);
256 : }
257 :
258 : /*@
259 : EPSRQCGSetReset - Sets the reset parameter of the RQCG iteration. Every
260 : nrest iterations, the solver performs a Rayleigh-Ritz projection step.
261 :
262 : Logically Collective
263 :
264 : Input Parameters:
265 : + eps - the eigenproblem solver context
266 : - nrest - the number of iterations between resets
267 :
268 : Options Database Key:
269 : . -eps_rqcg_reset - Sets the reset parameter
270 :
271 : Level: advanced
272 :
273 : .seealso: EPSRQCGGetReset()
274 : @*/
275 2 : PetscErrorCode EPSRQCGSetReset(EPS eps,PetscInt nrest)
276 : {
277 2 : PetscFunctionBegin;
278 2 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
279 8 : PetscValidLogicalCollectiveInt(eps,nrest,2);
280 2 : PetscTryMethod(eps,"EPSRQCGSetReset_C",(EPS,PetscInt),(eps,nrest));
281 2 : PetscFunctionReturn(PETSC_SUCCESS);
282 : }
283 :
284 1 : static PetscErrorCode EPSRQCGGetReset_RQCG(EPS eps,PetscInt *nrest)
285 : {
286 1 : EPS_RQCG *ctx = (EPS_RQCG*)eps->data;
287 :
288 1 : PetscFunctionBegin;
289 1 : *nrest = ctx->nrest;
290 1 : PetscFunctionReturn(PETSC_SUCCESS);
291 : }
292 :
293 : /*@
294 : EPSRQCGGetReset - Gets the reset parameter used in the RQCG method.
295 :
296 : Not Collective
297 :
298 : Input Parameter:
299 : . eps - the eigenproblem solver context
300 :
301 : Output Parameter:
302 : . nrest - the reset parameter
303 :
304 : Level: advanced
305 :
306 : .seealso: EPSRQCGSetReset()
307 : @*/
308 1 : PetscErrorCode EPSRQCGGetReset(EPS eps,PetscInt *nrest)
309 : {
310 1 : PetscFunctionBegin;
311 1 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
312 1 : PetscAssertPointer(nrest,2);
313 1 : PetscUseMethod(eps,"EPSRQCGGetReset_C",(EPS,PetscInt*),(eps,nrest));
314 1 : PetscFunctionReturn(PETSC_SUCCESS);
315 : }
316 :
317 11 : static PetscErrorCode EPSReset_RQCG(EPS eps)
318 : {
319 11 : EPS_RQCG *ctx = (EPS_RQCG*)eps->data;
320 :
321 11 : PetscFunctionBegin;
322 11 : PetscCall(BVDestroy(&ctx->AV));
323 11 : PetscCall(BVDestroy(&ctx->W));
324 11 : PetscCall(BVDestroy(&ctx->P));
325 11 : PetscCall(BVDestroy(&ctx->G));
326 11 : ctx->allocsize = 0;
327 11 : PetscFunctionReturn(PETSC_SUCCESS);
328 : }
329 :
330 10 : static PetscErrorCode EPSSetFromOptions_RQCG(EPS eps,PetscOptionItems *PetscOptionsObject)
331 : {
332 10 : PetscBool flg;
333 10 : PetscInt nrest;
334 :
335 10 : PetscFunctionBegin;
336 10 : PetscOptionsHeadBegin(PetscOptionsObject,"EPS RQCG Options");
337 :
338 10 : PetscCall(PetscOptionsInt("-eps_rqcg_reset","Reset parameter","EPSRQCGSetReset",20,&nrest,&flg));
339 10 : if (flg) PetscCall(EPSRQCGSetReset(eps,nrest));
340 :
341 10 : PetscOptionsHeadEnd();
342 10 : PetscFunctionReturn(PETSC_SUCCESS);
343 : }
344 :
345 11 : static PetscErrorCode EPSDestroy_RQCG(EPS eps)
346 : {
347 11 : PetscFunctionBegin;
348 11 : PetscCall(PetscFree(eps->data));
349 11 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSRQCGSetReset_C",NULL));
350 11 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSRQCGGetReset_C",NULL));
351 11 : PetscFunctionReturn(PETSC_SUCCESS);
352 : }
353 :
354 1 : static PetscErrorCode EPSView_RQCG(EPS eps,PetscViewer viewer)
355 : {
356 1 : EPS_RQCG *ctx = (EPS_RQCG*)eps->data;
357 1 : PetscBool isascii;
358 :
359 1 : PetscFunctionBegin;
360 1 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
361 1 : if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer," reset every %" PetscInt_FMT " iterations\n",ctx->nrest));
362 1 : PetscFunctionReturn(PETSC_SUCCESS);
363 : }
364 :
365 11 : SLEPC_EXTERN PetscErrorCode EPSCreate_RQCG(EPS eps)
366 : {
367 11 : EPS_RQCG *rqcg;
368 :
369 11 : PetscFunctionBegin;
370 11 : PetscCall(PetscNew(&rqcg));
371 11 : eps->data = (void*)rqcg;
372 :
373 11 : eps->useds = PETSC_TRUE;
374 11 : eps->categ = EPS_CATEGORY_PRECOND;
375 :
376 11 : eps->ops->solve = EPSSolve_RQCG;
377 11 : eps->ops->setup = EPSSetUp_RQCG;
378 11 : eps->ops->setupsort = EPSSetUpSort_Default;
379 11 : eps->ops->setfromoptions = EPSSetFromOptions_RQCG;
380 11 : eps->ops->destroy = EPSDestroy_RQCG;
381 11 : eps->ops->reset = EPSReset_RQCG;
382 11 : eps->ops->view = EPSView_RQCG;
383 11 : eps->ops->backtransform = EPSBackTransform_Default;
384 11 : eps->ops->setdefaultst = EPSSetDefaultST_GMRES;
385 :
386 11 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSRQCGSetReset_C",EPSRQCGSetReset_RQCG));
387 11 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSRQCGGetReset_C",EPSRQCGGetReset_RQCG));
388 11 : PetscFunctionReturn(PETSC_SUCCESS);
389 : }
|