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