Actual source code: gd.c

slepc-3.17.2 2022-08-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: PetscErrorCode EPSSetFromOptions_GD(PetscOptionItems *PetscOptionsObject,EPS eps)
 34: {
 35:   PetscBool      flg,flg2,op,orth;
 36:   PetscInt       opi,opi0;

 38:   PetscOptionsHead(PetscOptionsObject,"EPS Generalized Davidson (GD) Options");

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

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

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

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

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

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

 64:   PetscOptionsTail();
 65:   PetscFunctionReturn(0);
 66: }

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

 73:   /* Setup common for all davidson solvers */
 74:   EPSSetUp_XD(eps);

 76:   /* Check some constraints */
 77:   STGetKSP(eps->st,&ksp);
 78:   PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&t);
 80:   PetscFunctionReturn(0);
 81: }

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

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

108: PetscErrorCode EPSDestroy_GD(EPS eps)
109: {
110:   PetscFree(eps->data);
111:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetKrylovStart_C",NULL);
112:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetKrylovStart_C",NULL);
113:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBOrth_C",NULL);
114:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBOrth_C",NULL);
115:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBlockSize_C",NULL);
116:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBlockSize_C",NULL);
117:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetRestart_C",NULL);
118:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetRestart_C",NULL);
119:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetInitialSize_C",NULL);
120:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetInitialSize_C",NULL);
121:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetDoubleExpansion_C",NULL);
122:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetDoubleExpansion_C",NULL);
123:   PetscFunctionReturn(0);
124: }

126: /*@
127:    EPSGDSetKrylovStart - Activates or deactivates starting the searching
128:    subspace with a Krylov basis.

130:    Logically Collective on eps

132:    Input Parameters:
133: +  eps - the eigenproblem solver context
134: -  krylovstart - boolean flag

136:    Options Database Key:
137: .  -eps_gd_krylov_start - Activates starting the searching subspace with a
138:     Krylov basis

140:    Level: advanced

142: .seealso: EPSGDGetKrylovStart()
143: @*/
144: PetscErrorCode EPSGDSetKrylovStart(EPS eps,PetscBool krylovstart)
145: {
148:   PetscTryMethod(eps,"EPSGDSetKrylovStart_C",(EPS,PetscBool),(eps,krylovstart));
149:   PetscFunctionReturn(0);
150: }

152: /*@
153:    EPSGDGetKrylovStart - Returns a flag indicating if the search subspace is started with a
154:    Krylov basis.

156:    Not Collective

158:    Input Parameter:
159: .  eps - the eigenproblem solver context

161:    Output Parameters:
162: .  krylovstart - boolean flag indicating if the search subspace is started
163:    with a Krylov basis

165:    Level: advanced

167: .seealso: EPSGDSetKrylovStart()
168: @*/
169: PetscErrorCode EPSGDGetKrylovStart(EPS eps,PetscBool *krylovstart)
170: {
173:   PetscUseMethod(eps,"EPSGDGetKrylovStart_C",(EPS,PetscBool*),(eps,krylovstart));
174:   PetscFunctionReturn(0);
175: }

177: /*@
178:    EPSGDSetBlockSize - Sets the number of vectors to be added to the searching space
179:    in every iteration.

181:    Logically Collective on eps

183:    Input Parameters:
184: +  eps - the eigenproblem solver context
185: -  blocksize - number of vectors added to the search space in every iteration

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

190:    Level: advanced

192: .seealso: EPSGDSetKrylovStart()
193: @*/
194: PetscErrorCode EPSGDSetBlockSize(EPS eps,PetscInt blocksize)
195: {
198:   PetscTryMethod(eps,"EPSGDSetBlockSize_C",(EPS,PetscInt),(eps,blocksize));
199:   PetscFunctionReturn(0);
200: }

202: /*@
203:    EPSGDGetBlockSize - Returns the number of vectors to be added to the searching space
204:    in every iteration.

206:    Not Collective

208:    Input Parameter:
209: .  eps - the eigenproblem solver context

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

214:    Level: advanced

216: .seealso: EPSGDSetBlockSize()
217: @*/
218: PetscErrorCode EPSGDGetBlockSize(EPS eps,PetscInt *blocksize)
219: {
222:   PetscUseMethod(eps,"EPSGDGetBlockSize_C",(EPS,PetscInt*),(eps,blocksize));
223:   PetscFunctionReturn(0);
224: }

226: /*@
227:    EPSGDSetRestart - Sets the number of vectors of the searching space after
228:    restarting and the number of vectors saved from the previous iteration.

230:    Logically Collective on eps

232:    Input Parameters:
233: +  eps - the eigenproblem solver context
234: .  minv - number of vectors of the searching subspace after restarting
235: -  plusk - number of vectors saved from the previous iteration

237:    Options Database Keys:
238: +  -eps_gd_minv - number of vectors of the searching subspace after restarting
239: -  -eps_gd_plusk - number of vectors saved from the previous iteration

241:    Level: advanced

243: .seealso: EPSGDGetRestart()
244: @*/
245: PetscErrorCode EPSGDSetRestart(EPS eps,PetscInt minv,PetscInt plusk)
246: {
250:   PetscTryMethod(eps,"EPSGDSetRestart_C",(EPS,PetscInt,PetscInt),(eps,minv,plusk));
251:   PetscFunctionReturn(0);
252: }

254: /*@
255:    EPSGDGetRestart - Gets the number of vectors of the searching space after
256:    restarting and the number of vectors saved from the previous iteration.

258:    Not Collective

260:    Input Parameter:
261: .  eps - the eigenproblem solver context

263:    Output Parameters:
264: +  minv - number of vectors of the searching subspace after restarting
265: -  plusk - number of vectors saved from the previous iteration

267:    Level: advanced

269: .seealso: EPSGDSetRestart()
270: @*/
271: PetscErrorCode EPSGDGetRestart(EPS eps,PetscInt *minv,PetscInt *plusk)
272: {
274:   PetscUseMethod(eps,"EPSGDGetRestart_C",(EPS,PetscInt*,PetscInt*),(eps,minv,plusk));
275:   PetscFunctionReturn(0);
276: }

278: /*@
279:    EPSGDSetInitialSize - Sets the initial size of the searching space.

281:    Logically Collective on eps

283:    Input Parameters:
284: +  eps - the eigenproblem solver context
285: -  initialsize - number of vectors of the initial searching subspace

287:    Options Database Key:
288: .  -eps_gd_initial_size - number of vectors of the initial searching subspace

290:    Notes:
291:    If EPSGDGetKrylovStart() is PETSC_FALSE and the user provides vectors with
292:    EPSSetInitialSpace(), up to initialsize vectors will be used; and if the
293:    provided vectors are not enough, the solver completes the subspace with
294:    random vectors. In the case of EPSGDGetKrylovStart() being PETSC_TRUE, the solver
295:    gets the first vector provided by the user or, if not available, a random vector,
296:    and expands the Krylov basis up to initialsize vectors.

298:    Level: advanced

300: .seealso: EPSGDGetInitialSize(), EPSGDGetKrylovStart()
301: @*/
302: PetscErrorCode EPSGDSetInitialSize(EPS eps,PetscInt initialsize)
303: {
306:   PetscTryMethod(eps,"EPSGDSetInitialSize_C",(EPS,PetscInt),(eps,initialsize));
307:   PetscFunctionReturn(0);
308: }

310: /*@
311:    EPSGDGetInitialSize - Returns the initial size of the searching space.

313:    Not Collective

315:    Input Parameter:
316: .  eps - the eigenproblem solver context

318:    Output Parameter:
319: .  initialsize - number of vectors of the initial searching subspace

321:    Notes:
322:    If EPSGDGetKrylovStart() is PETSC_FALSE and the user provides vectors with
323:    EPSSetInitialSpace(), up to initialsize vectors will be used; and if the
324:    provided vectors are not enough, the solver completes the subspace with
325:    random vectors. In the case of EPSGDGetKrylovStart() being PETSC_TRUE, the solver
326:    gets the first vector provided by the user or, if not available, a random vector,
327:    and expands the Krylov basis up to initialsize vectors.

329:    Level: advanced

331: .seealso: EPSGDSetInitialSize(), EPSGDGetKrylovStart()
332: @*/
333: PetscErrorCode EPSGDGetInitialSize(EPS eps,PetscInt *initialsize)
334: {
337:   PetscUseMethod(eps,"EPSGDGetInitialSize_C",(EPS,PetscInt*),(eps,initialsize));
338:   PetscFunctionReturn(0);
339: }

341: /*@
342:    EPSGDSetBOrth - Selects the orthogonalization that will be used in the search
343:    subspace in case of generalized Hermitian problems.

345:    Logically Collective on eps

347:    Input Parameters:
348: +  eps   - the eigenproblem solver context
349: -  borth - whether to B-orthogonalize the search subspace

351:    Options Database Key:
352: .  -eps_gd_borth - Set the orthogonalization used in the search subspace

354:    Level: advanced

356: .seealso: EPSGDGetBOrth()
357: @*/
358: PetscErrorCode EPSGDSetBOrth(EPS eps,PetscBool borth)
359: {
362:   PetscTryMethod(eps,"EPSGDSetBOrth_C",(EPS,PetscBool),(eps,borth));
363:   PetscFunctionReturn(0);
364: }

366: /*@
367:    EPSGDGetBOrth - Returns the orthogonalization used in the search
368:    subspace in case of generalized Hermitian problems.

370:    Not Collective

372:    Input Parameter:
373: .  eps - the eigenproblem solver context

375:    Output Parameters:
376: .  borth - whether to B-orthogonalize the search subspace

378:    Level: advanced

380: .seealso: EPSGDSetBOrth()
381: @*/
382: PetscErrorCode EPSGDGetBOrth(EPS eps,PetscBool *borth)
383: {
386:   PetscUseMethod(eps,"EPSGDGetBOrth_C",(EPS,PetscBool*),(eps,borth));
387:   PetscFunctionReturn(0);
388: }

390: static PetscErrorCode EPSGDSetDoubleExpansion_GD(EPS eps,PetscBool doubleexp)
391: {
392:   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;

394:   data->doubleexp = doubleexp;
395:   PetscFunctionReturn(0);
396: }

398: /*@
399:    EPSGDSetDoubleExpansion - Activate the double expansion variant of GD.

401:    Logically Collective on eps

403:    Input Parameters:
404: +  eps - the eigenproblem solver context
405: -  doubleexp - the boolean flag

407:    Options Database Keys:
408: .  -eps_gd_double_expansion - activate the double-expansion variant of GD

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

415:    Level: advanced

417: .seealso: EPSGDGetDoubleExpansion()
418: @*/
419: PetscErrorCode EPSGDSetDoubleExpansion(EPS eps,PetscBool doubleexp)
420: {
423:   PetscTryMethod(eps,"EPSGDSetDoubleExpansion_C",(EPS,PetscBool),(eps,doubleexp));
424:   PetscFunctionReturn(0);
425: }

427: static PetscErrorCode EPSGDGetDoubleExpansion_GD(EPS eps,PetscBool *doubleexp)
428: {
429:   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;

431:   *doubleexp = data->doubleexp;
432:   PetscFunctionReturn(0);
433: }

435: /*@
436:    EPSGDGetDoubleExpansion - Gets a flag indicating whether the double
437:    expansion variant has been activated or not.

439:    Not Collective

441:    Input Parameter:
442: .  eps - the eigenproblem solver context

444:    Output Parameter:
445: .  doubleexp - the flag

447:    Level: advanced

449: .seealso: EPSGDSetDoubleExpansion()
450: @*/
451: PetscErrorCode EPSGDGetDoubleExpansion(EPS eps,PetscBool *doubleexp)
452: {
455:   PetscUseMethod(eps,"EPSGDGetDoubleExpansion_C",(EPS,PetscBool*),(eps,doubleexp));
456:   PetscFunctionReturn(0);
457: }

459: SLEPC_EXTERN PetscErrorCode EPSCreate_GD(EPS eps)
460: {
461:   EPS_DAVIDSON    *data;

463:   PetscNewLog(eps,&data);
464:   eps->data = (void*)data;

466:   data->blocksize   = 1;
467:   data->initialsize = 0;
468:   data->minv        = 0;
469:   data->plusk       = PETSC_DEFAULT;
470:   data->ipB         = PETSC_TRUE;
471:   data->fix         = 0.0;
472:   data->krylovstart = PETSC_FALSE;
473:   data->dynamic     = PETSC_FALSE;

475:   eps->useds = PETSC_TRUE;
476:   eps->categ = EPS_CATEGORY_PRECOND;

478:   eps->ops->solve          = EPSSolve_XD;
479:   eps->ops->setup          = EPSSetUp_GD;
480:   eps->ops->setupsort      = EPSSetUpSort_Default;
481:   eps->ops->setfromoptions = EPSSetFromOptions_GD;
482:   eps->ops->destroy        = EPSDestroy_GD;
483:   eps->ops->reset          = EPSReset_XD;
484:   eps->ops->view           = EPSView_GD;
485:   eps->ops->backtransform  = EPSBackTransform_Default;
486:   eps->ops->computevectors = EPSComputeVectors_XD;
487:   eps->ops->setdefaultst   = EPSSetDefaultST_Precond;

489:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetKrylovStart_C",EPSXDSetKrylovStart_XD);
490:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetKrylovStart_C",EPSXDGetKrylovStart_XD);
491:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBOrth_C",EPSXDSetBOrth_XD);
492:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBOrth_C",EPSXDGetBOrth_XD);
493:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBlockSize_C",EPSXDSetBlockSize_XD);
494:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBlockSize_C",EPSXDGetBlockSize_XD);
495:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetRestart_C",EPSXDSetRestart_XD);
496:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetRestart_C",EPSXDGetRestart_XD);
497:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetInitialSize_C",EPSXDSetInitialSize_XD);
498:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetInitialSize_C",EPSXDGetInitialSize_XD);
499:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetDoubleExpansion_C",EPSGDSetDoubleExpansion_GD);
500:   PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetDoubleExpansion_C",EPSGDGetDoubleExpansion_GD);
501:   PetscFunctionReturn(0);
502: }