Actual source code: gd.c

slepc-main 2024-11-09
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  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: "gd"

 13:    Method: Generalized Davidson

 15:    Algorithm:

 17:        Generalized Davidson with various subspace extraction and
 18:        restart techniques.

 20:    References:

 22:        [1] E.R. Davidson, "Super-matrix methods", Comput. Phys. Commun.
 23:            53(2):49-60, 1989.

 25:        [2] E. Romero and J.E. Roman, "A parallel implementation of
 26:            Davidson methods for large-scale eigenvalue problems in
 27:            SLEPc", ACM Trans. Math. Software 40(2), Article 13, 2014.
 28: */

 30: #include <slepc/private/epsimpl.h>
 31: #include <../src/eps/impls/davidson/davidson.h>

 33: static PetscErrorCode EPSSetFromOptions_GD(EPS eps,PetscOptionItems *PetscOptionsObject)
 34: {
 35:   PetscBool      flg,flg2,op,orth;
 36:   PetscInt       opi,opi0;

 38:   PetscFunctionBegin;
 39:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS Generalized Davidson (GD) Options");

 41:     PetscCall(EPSGDGetKrylovStart(eps,&op));
 42:     PetscCall(PetscOptionsBool("-eps_gd_krylov_start","Start the search subspace with a Krylov basis","EPSGDSetKrylovStart",op,&op,&flg));
 43:     if (flg) PetscCall(EPSGDSetKrylovStart(eps,op));

 45:     PetscCall(EPSGDGetBOrth(eps,&orth));
 46:     PetscCall(PetscOptionsBool("-eps_gd_borth","Use B-orthogonalization in the search subspace","EPSGDSetBOrth",op,&op,&flg));
 47:     if (flg) PetscCall(EPSGDSetBOrth(eps,op));

 49:     PetscCall(EPSGDGetBlockSize(eps,&opi));
 50:     PetscCall(PetscOptionsInt("-eps_gd_blocksize","Number of vectors to add to the search subspace","EPSGDSetBlockSize",opi,&opi,&flg));
 51:     if (flg) PetscCall(EPSGDSetBlockSize(eps,opi));

 53:     PetscCall(EPSGDGetRestart(eps,&opi,&opi0));
 54:     PetscCall(PetscOptionsInt("-eps_gd_minv","Size of the search subspace after restarting","EPSGDSetRestart",opi,&opi,&flg));
 55:     PetscCall(PetscOptionsInt("-eps_gd_plusk","Number of eigenvectors saved from the previous iteration when restarting","EPSGDSetRestart",opi0,&opi0,&flg2));
 56:     if (flg || flg2) PetscCall(EPSGDSetRestart(eps,opi,opi0));

 58:     PetscCall(EPSGDGetInitialSize(eps,&opi));
 59:     PetscCall(PetscOptionsInt("-eps_gd_initial_size","Initial size of the search subspace","EPSGDSetInitialSize",opi,&opi,&flg));
 60:     if (flg) PetscCall(EPSGDSetInitialSize(eps,opi));

 62:     PetscCall(PetscOptionsBool("-eps_gd_double_expansion","Use the doble-expansion variant of GD","EPSGDSetDoubleExpansion",PETSC_FALSE,&op,&flg));
 63:     if (flg) PetscCall(EPSGDSetDoubleExpansion(eps,op));

 65:   PetscOptionsHeadEnd();
 66:   PetscFunctionReturn(PETSC_SUCCESS);
 67: }

 69: static PetscErrorCode EPSSetUp_GD(EPS eps)
 70: {
 71:   PetscBool      t;
 72:   KSP            ksp;

 74:   PetscFunctionBegin;
 75:   /* Setup common for all davidson solvers */
 76:   PetscCall(EPSSetUp_XD(eps));

 78:   /* Check some constraints */
 79:   PetscCall(STGetKSP(eps->st,&ksp));
 80:   PetscCall(PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&t));
 81:   PetscCheck(t,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"EPSGD only works with KSPPREONLY");
 82:   PetscFunctionReturn(PETSC_SUCCESS);
 83: }

 85: static PetscErrorCode EPSView_GD(EPS eps,PetscViewer viewer)
 86: {
 87:   PetscBool      isascii,opb;
 88:   PetscInt       opi,opi0;
 89:   PetscBool      borth;
 90:   EPS_DAVIDSON   *data = (EPS_DAVIDSON*)eps->data;

 92:   PetscFunctionBegin;
 93:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
 94:   if (isascii) {
 95:     if (data->doubleexp) PetscCall(PetscViewerASCIIPrintf(viewer,"  using double expansion variant (GD2)\n"));
 96:     PetscCall(EPSXDGetBOrth_XD(eps,&borth));
 97:     if (borth) PetscCall(PetscViewerASCIIPrintf(viewer,"  search subspace is B-orthogonalized\n"));
 98:     else PetscCall(PetscViewerASCIIPrintf(viewer,"  search subspace is orthogonalized\n"));
 99:     PetscCall(EPSXDGetBlockSize_XD(eps,&opi));
100:     PetscCall(PetscViewerASCIIPrintf(viewer,"  block size=%" PetscInt_FMT "\n",opi));
101:     PetscCall(EPSXDGetKrylovStart_XD(eps,&opb));
102:     if (!opb) PetscCall(PetscViewerASCIIPrintf(viewer,"  type of the initial subspace: non-Krylov\n"));
103:     else PetscCall(PetscViewerASCIIPrintf(viewer,"  type of the initial subspace: Krylov\n"));
104:     PetscCall(EPSXDGetRestart_XD(eps,&opi,&opi0));
105:     PetscCall(PetscViewerASCIIPrintf(viewer,"  size of the subspace after restarting: %" PetscInt_FMT "\n",opi));
106:     PetscCall(PetscViewerASCIIPrintf(viewer,"  number of vectors after restarting from the previous iteration: %" PetscInt_FMT "\n",opi0));
107:   }
108:   PetscFunctionReturn(PETSC_SUCCESS);
109: }

111: static PetscErrorCode EPSDestroy_GD(EPS eps)
112: {
113:   PetscFunctionBegin;
114:   PetscCall(PetscFree(eps->data));
115:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetKrylovStart_C",NULL));
116:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetKrylovStart_C",NULL));
117:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBOrth_C",NULL));
118:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBOrth_C",NULL));
119:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBlockSize_C",NULL));
120:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBlockSize_C",NULL));
121:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetRestart_C",NULL));
122:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetRestart_C",NULL));
123:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetInitialSize_C",NULL));
124:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetInitialSize_C",NULL));
125:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetDoubleExpansion_C",NULL));
126:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetDoubleExpansion_C",NULL));
127:   PetscFunctionReturn(PETSC_SUCCESS);
128: }

130: /*@
131:    EPSGDSetKrylovStart - Activates or deactivates starting the searching
132:    subspace with a Krylov basis.

134:    Logically Collective

136:    Input Parameters:
137: +  eps - the eigenproblem solver context
138: -  krylovstart - boolean flag

140:    Options Database Key:
141: .  -eps_gd_krylov_start - Activates starting the searching subspace with a
142:     Krylov basis

144:    Level: advanced

146: .seealso: EPSGDGetKrylovStart()
147: @*/
148: PetscErrorCode EPSGDSetKrylovStart(EPS eps,PetscBool krylovstart)
149: {
150:   PetscFunctionBegin;
153:   PetscTryMethod(eps,"EPSGDSetKrylovStart_C",(EPS,PetscBool),(eps,krylovstart));
154:   PetscFunctionReturn(PETSC_SUCCESS);
155: }

157: /*@
158:    EPSGDGetKrylovStart - Returns a flag indicating if the search subspace is started with a
159:    Krylov basis.

161:    Not Collective

163:    Input Parameter:
164: .  eps - the eigenproblem solver context

166:    Output Parameters:
167: .  krylovstart - boolean flag indicating if the search subspace is started
168:    with a Krylov basis

170:    Level: advanced

172: .seealso: EPSGDSetKrylovStart()
173: @*/
174: PetscErrorCode EPSGDGetKrylovStart(EPS eps,PetscBool *krylovstart)
175: {
176:   PetscFunctionBegin;
178:   PetscAssertPointer(krylovstart,2);
179:   PetscUseMethod(eps,"EPSGDGetKrylovStart_C",(EPS,PetscBool*),(eps,krylovstart));
180:   PetscFunctionReturn(PETSC_SUCCESS);
181: }

183: /*@
184:    EPSGDSetBlockSize - Sets the number of vectors to be added to the searching space
185:    in every iteration.

187:    Logically Collective

189:    Input Parameters:
190: +  eps - the eigenproblem solver context
191: -  blocksize - number of vectors added to the search space in every iteration

193:    Options Database Key:
194: .  -eps_gd_blocksize - number of vectors added to the search space in every iteration

196:    Level: advanced

198: .seealso: EPSGDSetKrylovStart()
199: @*/
200: PetscErrorCode EPSGDSetBlockSize(EPS eps,PetscInt blocksize)
201: {
202:   PetscFunctionBegin;
205:   PetscTryMethod(eps,"EPSGDSetBlockSize_C",(EPS,PetscInt),(eps,blocksize));
206:   PetscFunctionReturn(PETSC_SUCCESS);
207: }

209: /*@
210:    EPSGDGetBlockSize - Returns the number of vectors to be added to the searching space
211:    in every iteration.

213:    Not Collective

215:    Input Parameter:
216: .  eps - the eigenproblem solver context

218:    Output Parameter:
219: .  blocksize - number of vectors added to the search space in every iteration

221:    Level: advanced

223: .seealso: EPSGDSetBlockSize()
224: @*/
225: PetscErrorCode EPSGDGetBlockSize(EPS eps,PetscInt *blocksize)
226: {
227:   PetscFunctionBegin;
229:   PetscAssertPointer(blocksize,2);
230:   PetscUseMethod(eps,"EPSGDGetBlockSize_C",(EPS,PetscInt*),(eps,blocksize));
231:   PetscFunctionReturn(PETSC_SUCCESS);
232: }

234: /*@
235:    EPSGDSetRestart - Sets the number of vectors of the searching space after
236:    restarting and the number of vectors saved from the previous iteration.

238:    Logically Collective

240:    Input Parameters:
241: +  eps - the eigenproblem solver context
242: .  minv - number of vectors of the searching subspace after restarting
243: -  plusk - number of vectors saved from the previous iteration

245:    Options Database Keys:
246: +  -eps_gd_minv - number of vectors of the searching subspace after restarting
247: -  -eps_gd_plusk - number of vectors saved from the previous iteration

249:    Note:
250:    PETSC_CURRENT can be used to preserve the current value of any of the
251:    arguments, and PETSC_DETERMINE to set them to a default value.

253:    Level: advanced

255: .seealso: EPSGDGetRestart()
256: @*/
257: PetscErrorCode EPSGDSetRestart(EPS eps,PetscInt minv,PetscInt plusk)
258: {
259:   PetscFunctionBegin;
263:   PetscTryMethod(eps,"EPSGDSetRestart_C",(EPS,PetscInt,PetscInt),(eps,minv,plusk));
264:   PetscFunctionReturn(PETSC_SUCCESS);
265: }

267: /*@
268:    EPSGDGetRestart - Gets the number of vectors of the searching space after
269:    restarting and the number of vectors saved from the previous iteration.

271:    Not Collective

273:    Input Parameter:
274: .  eps - the eigenproblem solver context

276:    Output Parameters:
277: +  minv - number of vectors of the searching subspace after restarting
278: -  plusk - number of vectors saved from the previous iteration

280:    Level: advanced

282: .seealso: EPSGDSetRestart()
283: @*/
284: PetscErrorCode EPSGDGetRestart(EPS eps,PetscInt *minv,PetscInt *plusk)
285: {
286:   PetscFunctionBegin;
288:   PetscUseMethod(eps,"EPSGDGetRestart_C",(EPS,PetscInt*,PetscInt*),(eps,minv,plusk));
289:   PetscFunctionReturn(PETSC_SUCCESS);
290: }

292: /*@
293:    EPSGDSetInitialSize - Sets the initial size of the searching space.

295:    Logically Collective

297:    Input Parameters:
298: +  eps - the eigenproblem solver context
299: -  initialsize - number of vectors of the initial searching subspace

301:    Options Database Key:
302: .  -eps_gd_initial_size - number of vectors of the initial searching subspace

304:    Notes:
305:    If EPSGDGetKrylovStart() is PETSC_FALSE and the user provides vectors with
306:    EPSSetInitialSpace(), up to initialsize vectors will be used; and if the
307:    provided vectors are not enough, the solver completes the subspace with
308:    random vectors. In the case of EPSGDGetKrylovStart() being PETSC_TRUE, the solver
309:    gets the first vector provided by the user or, if not available, a random vector,
310:    and expands the Krylov basis up to initialsize vectors.

312:    Level: advanced

314: .seealso: EPSGDGetInitialSize(), EPSGDGetKrylovStart()
315: @*/
316: PetscErrorCode EPSGDSetInitialSize(EPS eps,PetscInt initialsize)
317: {
318:   PetscFunctionBegin;
321:   PetscTryMethod(eps,"EPSGDSetInitialSize_C",(EPS,PetscInt),(eps,initialsize));
322:   PetscFunctionReturn(PETSC_SUCCESS);
323: }

325: /*@
326:    EPSGDGetInitialSize - Returns the initial size of the searching space.

328:    Not Collective

330:    Input Parameter:
331: .  eps - the eigenproblem solver context

333:    Output Parameter:
334: .  initialsize - number of vectors of the initial searching subspace

336:    Notes:
337:    If EPSGDGetKrylovStart() is PETSC_FALSE and the user provides vectors with
338:    EPSSetInitialSpace(), up to initialsize vectors will be used; and if the
339:    provided vectors are not enough, the solver completes the subspace with
340:    random vectors. In the case of EPSGDGetKrylovStart() being PETSC_TRUE, the solver
341:    gets the first vector provided by the user or, if not available, a random vector,
342:    and expands the Krylov basis up to initialsize vectors.

344:    Level: advanced

346: .seealso: EPSGDSetInitialSize(), EPSGDGetKrylovStart()
347: @*/
348: PetscErrorCode EPSGDGetInitialSize(EPS eps,PetscInt *initialsize)
349: {
350:   PetscFunctionBegin;
352:   PetscAssertPointer(initialsize,2);
353:   PetscUseMethod(eps,"EPSGDGetInitialSize_C",(EPS,PetscInt*),(eps,initialsize));
354:   PetscFunctionReturn(PETSC_SUCCESS);
355: }

357: /*@
358:    EPSGDSetBOrth - Selects the orthogonalization that will be used in the search
359:    subspace in case of generalized Hermitian problems.

361:    Logically Collective

363:    Input Parameters:
364: +  eps   - the eigenproblem solver context
365: -  borth - whether to B-orthogonalize the search subspace

367:    Options Database Key:
368: .  -eps_gd_borth - Set the orthogonalization used in the search subspace

370:    Level: advanced

372: .seealso: EPSGDGetBOrth()
373: @*/
374: PetscErrorCode EPSGDSetBOrth(EPS eps,PetscBool borth)
375: {
376:   PetscFunctionBegin;
379:   PetscTryMethod(eps,"EPSGDSetBOrth_C",(EPS,PetscBool),(eps,borth));
380:   PetscFunctionReturn(PETSC_SUCCESS);
381: }

383: /*@
384:    EPSGDGetBOrth - Returns the orthogonalization used in the search
385:    subspace in case of generalized Hermitian problems.

387:    Not Collective

389:    Input Parameter:
390: .  eps - the eigenproblem solver context

392:    Output Parameters:
393: .  borth - whether to B-orthogonalize the search subspace

395:    Level: advanced

397: .seealso: EPSGDSetBOrth()
398: @*/
399: PetscErrorCode EPSGDGetBOrth(EPS eps,PetscBool *borth)
400: {
401:   PetscFunctionBegin;
403:   PetscAssertPointer(borth,2);
404:   PetscUseMethod(eps,"EPSGDGetBOrth_C",(EPS,PetscBool*),(eps,borth));
405:   PetscFunctionReturn(PETSC_SUCCESS);
406: }

408: static PetscErrorCode EPSGDSetDoubleExpansion_GD(EPS eps,PetscBool doubleexp)
409: {
410:   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;

412:   PetscFunctionBegin;
413:   data->doubleexp = doubleexp;
414:   PetscFunctionReturn(PETSC_SUCCESS);
415: }

417: /*@
418:    EPSGDSetDoubleExpansion - Activate the double expansion variant of GD.

420:    Logically Collective

422:    Input Parameters:
423: +  eps - the eigenproblem solver context
424: -  doubleexp - the boolean flag

426:    Options Database Keys:
427: .  -eps_gd_double_expansion - activate the double-expansion variant of GD

429:    Notes:
430:    In the double expansion variant the search subspace is expanded with K*[A*x B*x]
431:    instead of the classic K*r, where K is the preconditioner, x the selected
432:    approximate eigenvector and r its associated residual vector.

434:    Level: advanced

436: .seealso: EPSGDGetDoubleExpansion()
437: @*/
438: PetscErrorCode EPSGDSetDoubleExpansion(EPS eps,PetscBool doubleexp)
439: {
440:   PetscFunctionBegin;
443:   PetscTryMethod(eps,"EPSGDSetDoubleExpansion_C",(EPS,PetscBool),(eps,doubleexp));
444:   PetscFunctionReturn(PETSC_SUCCESS);
445: }

447: static PetscErrorCode EPSGDGetDoubleExpansion_GD(EPS eps,PetscBool *doubleexp)
448: {
449:   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;

451:   PetscFunctionBegin;
452:   *doubleexp = data->doubleexp;
453:   PetscFunctionReturn(PETSC_SUCCESS);
454: }

456: /*@
457:    EPSGDGetDoubleExpansion - Gets a flag indicating whether the double
458:    expansion variant has been activated or not.

460:    Not Collective

462:    Input Parameter:
463: .  eps - the eigenproblem solver context

465:    Output Parameter:
466: .  doubleexp - the flag

468:    Level: advanced

470: .seealso: EPSGDSetDoubleExpansion()
471: @*/
472: PetscErrorCode EPSGDGetDoubleExpansion(EPS eps,PetscBool *doubleexp)
473: {
474:   PetscFunctionBegin;
476:   PetscAssertPointer(doubleexp,2);
477:   PetscUseMethod(eps,"EPSGDGetDoubleExpansion_C",(EPS,PetscBool*),(eps,doubleexp));
478:   PetscFunctionReturn(PETSC_SUCCESS);
479: }

481: SLEPC_EXTERN PetscErrorCode EPSCreate_GD(EPS eps)
482: {
483:   EPS_DAVIDSON    *data;

485:   PetscFunctionBegin;
486:   PetscCall(PetscNew(&data));
487:   eps->data = (void*)data;

489:   data->blocksize   = 1;
490:   data->initialsize = 0;
491:   data->minv        = 0;
492:   data->plusk       = PETSC_DETERMINE;
493:   data->ipB         = PETSC_TRUE;
494:   data->fix         = 0.0;
495:   data->krylovstart = PETSC_FALSE;
496:   data->dynamic     = PETSC_FALSE;

498:   eps->useds = PETSC_TRUE;
499:   eps->categ = EPS_CATEGORY_PRECOND;

501:   eps->ops->solve          = EPSSolve_XD;
502:   eps->ops->setup          = EPSSetUp_GD;
503:   eps->ops->setupsort      = EPSSetUpSort_Default;
504:   eps->ops->setfromoptions = EPSSetFromOptions_GD;
505:   eps->ops->destroy        = EPSDestroy_GD;
506:   eps->ops->reset          = EPSReset_XD;
507:   eps->ops->view           = EPSView_GD;
508:   eps->ops->backtransform  = EPSBackTransform_Default;
509:   eps->ops->computevectors = EPSComputeVectors_XD;
510:   eps->ops->setdefaultst   = EPSSetDefaultST_Precond;

512:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetKrylovStart_C",EPSXDSetKrylovStart_XD));
513:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetKrylovStart_C",EPSXDGetKrylovStart_XD));
514:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBOrth_C",EPSXDSetBOrth_XD));
515:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBOrth_C",EPSXDGetBOrth_XD));
516:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBlockSize_C",EPSXDSetBlockSize_XD));
517:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBlockSize_C",EPSXDGetBlockSize_XD));
518:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetRestart_C",EPSXDSetRestart_XD));
519:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetRestart_C",EPSXDGetRestart_XD));
520:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetInitialSize_C",EPSXDSetInitialSize_XD));
521:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetInitialSize_C",EPSXDGetInitialSize_XD));
522:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetDoubleExpansion_C",EPSGDSetDoubleExpansion_GD));
523:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetDoubleExpansion_C",EPSGDGetDoubleExpansion_GD));
524:   PetscFunctionReturn(PETSC_SUCCESS);
525: }