Actual source code: arnoldi.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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: "arnoldi"
13: Method: Explicitly Restarted Arnoldi
15: Algorithm:
17: Arnoldi method with explicit restart and deflation.
19: References:
21: [1] "Arnoldi Methods in SLEPc", SLEPc Technical Report STR-4,
22: available at http://slepc.upv.es.
23: */
25: #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
27: typedef struct {
28: PetscBool delayed;
29: } EPS_ARNOLDI;
31: PetscErrorCode EPSSetUp_Arnoldi(EPS eps)
32: {
36: EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
37: if (eps->ncv>eps->nev+eps->mpd) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must not be larger than nev+mpd");
38: if (!eps->max_it) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
39: if (!eps->which) { EPSSetWhichEigenpairs_Default(eps); }
40: if (eps->which==EPS_ALL || (eps->ishermitian && eps->ispositive && (eps->which==EPS_LARGEST_IMAGINARY || eps->which==EPS_SMALLEST_IMAGINARY))) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Wrong value of eps->which");
42: if (!eps->extraction) {
43: EPSSetExtraction(eps,EPS_RITZ);
44: }
45: if (eps->arbitrary) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Arbitrary selection of eigenpairs not supported in this solver");
47: EPSAllocateSolution(eps,1);
48: EPS_SetInnerProduct(eps);
49: DSSetType(eps->ds,DSNHEP);
50: if (eps->extraction==EPS_REFINED || eps->extraction==EPS_REFINED_HARMONIC) {
51: DSSetRefined(eps->ds,PETSC_TRUE);
52: }
53: DSSetExtraRow(eps->ds,PETSC_TRUE);
54: DSAllocate(eps->ds,eps->ncv+1);
56: if (eps->isgeneralized && eps->ishermitian && !eps->ispositive) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Requested method does not work for indefinite problems");
57: return(0);
58: }
60: PetscErrorCode EPSSolve_Arnoldi(EPS eps)
61: {
62: PetscErrorCode ierr;
63: PetscInt k,nv,ld;
64: Mat U;
65: PetscScalar *H;
66: PetscReal beta,gamma=1.0;
67: PetscBool breakdown,harmonic,refined;
68: BVOrthogRefineType orthog_ref;
69: EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
72: DSGetLeadingDimension(eps->ds,&ld);
73: DSGetRefined(eps->ds,&refined);
74: harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
75: BVGetOrthogonalization(eps->V,NULL,&orthog_ref,NULL,NULL);
77: /* Get the starting Arnoldi vector */
78: EPSGetStartVector(eps,0,NULL);
80: /* Restart loop */
81: while (eps->reason == EPS_CONVERGED_ITERATING) {
82: eps->its++;
84: /* Compute an nv-step Arnoldi factorization */
85: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
86: DSSetDimensions(eps->ds,nv,0,eps->nconv,0);
87: DSGetArray(eps->ds,DS_MAT_A,&H);
88: if (!arnoldi->delayed) {
89: EPSBasicArnoldi(eps,PETSC_FALSE,H,ld,eps->nconv,&nv,&beta,&breakdown);
90: } else if (orthog_ref == BV_ORTHOG_REFINE_NEVER) {
91: EPSDelayedArnoldi1(eps,H,ld,eps->nconv,&nv,&beta,&breakdown);
92: } else {
93: EPSDelayedArnoldi(eps,H,ld,eps->nconv,&nv,&beta,&breakdown);
94: }
95: DSRestoreArray(eps->ds,DS_MAT_A,&H);
96: DSSetState(eps->ds,DS_STATE_INTERMEDIATE);
97: BVSetActiveColumns(eps->V,eps->nconv,nv);
99: /* Compute translation of Krylov decomposition if harmonic extraction used */
100: if (harmonic) {
101: DSTranslateHarmonic(eps->ds,eps->target,beta,PETSC_FALSE,NULL,&gamma);
102: }
104: /* Solve projected problem */
105: DSSolve(eps->ds,eps->eigr,eps->eigi);
106: DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL);
107: DSUpdateExtraRow(eps->ds);
108: DSSynchronize(eps->ds,eps->eigr,eps->eigi);
110: /* Check convergence */
111: EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k);
112: if (refined) {
113: DSGetMat(eps->ds,DS_MAT_X,&U);
114: BVMultInPlace(eps->V,U,eps->nconv,k+1);
115: MatDestroy(&U);
116: BVOrthonormalizeColumn(eps->V,k,PETSC_FALSE,NULL,NULL);
117: } else {
118: DSGetMat(eps->ds,DS_MAT_Q,&U);
119: BVMultInPlace(eps->V,U,eps->nconv,PetscMin(k+1,nv));
120: MatDestroy(&U);
121: }
122: (*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx);
123: if (eps->reason == EPS_CONVERGED_ITERATING && breakdown) {
124: PetscInfo2(eps,"Breakdown in Arnoldi method (it=%D norm=%g)\n",eps->its,(double)beta);
125: EPSGetStartVector(eps,k,&breakdown);
126: if (breakdown) {
127: eps->reason = EPS_DIVERGED_BREAKDOWN;
128: PetscInfo(eps,"Unable to generate more start vectors\n");
129: }
130: }
131: eps->nconv = k;
132: EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,nv);
133: }
135: /* truncate Schur decomposition and change the state to raw so that
136: DSVectors() computes eigenvectors from scratch */
137: DSSetDimensions(eps->ds,eps->nconv,0,0,0);
138: DSSetState(eps->ds,DS_STATE_RAW);
139: return(0);
140: }
142: PetscErrorCode EPSSetFromOptions_Arnoldi(PetscOptionItems *PetscOptionsObject,EPS eps)
143: {
145: PetscBool set,val;
146: EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
149: PetscOptionsHead(PetscOptionsObject,"EPS Arnoldi Options");
151: PetscOptionsBool("-eps_arnoldi_delayed","Use delayed reorthogonalization","EPSArnoldiSetDelayed",arnoldi->delayed,&val,&set);
152: if (set) { EPSArnoldiSetDelayed(eps,val); }
154: PetscOptionsTail();
155: return(0);
156: }
158: static PetscErrorCode EPSArnoldiSetDelayed_Arnoldi(EPS eps,PetscBool delayed)
159: {
160: EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
163: arnoldi->delayed = delayed;
164: return(0);
165: }
167: /*@
168: EPSArnoldiSetDelayed - Activates or deactivates delayed reorthogonalization
169: in the Arnoldi iteration.
171: Logically Collective on EPS
173: Input Parameters:
174: + eps - the eigenproblem solver context
175: - delayed - boolean flag
177: Options Database Key:
178: . -eps_arnoldi_delayed - Activates delayed reorthogonalization in Arnoldi
180: Note:
181: Delayed reorthogonalization is an aggressive optimization for the Arnoldi
182: eigensolver than may provide better scalability, but sometimes makes the
183: solver converge less than the default algorithm.
185: Level: advanced
187: .seealso: EPSArnoldiGetDelayed()
188: @*/
189: PetscErrorCode EPSArnoldiSetDelayed(EPS eps,PetscBool delayed)
190: {
196: PetscTryMethod(eps,"EPSArnoldiSetDelayed_C",(EPS,PetscBool),(eps,delayed));
197: return(0);
198: }
200: static PetscErrorCode EPSArnoldiGetDelayed_Arnoldi(EPS eps,PetscBool *delayed)
201: {
202: EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
205: *delayed = arnoldi->delayed;
206: return(0);
207: }
209: /*@
210: EPSArnoldiGetDelayed - Gets the type of reorthogonalization used during the Arnoldi
211: iteration.
213: Not Collective
215: Input Parameter:
216: . eps - the eigenproblem solver context
218: Input Parameter:
219: . delayed - boolean flag indicating if delayed reorthogonalization has been enabled
221: Level: advanced
223: .seealso: EPSArnoldiSetDelayed()
224: @*/
225: PetscErrorCode EPSArnoldiGetDelayed(EPS eps,PetscBool *delayed)
226: {
232: PetscUseMethod(eps,"EPSArnoldiGetDelayed_C",(EPS,PetscBool*),(eps,delayed));
233: return(0);
234: }
236: PetscErrorCode EPSDestroy_Arnoldi(EPS eps)
237: {
241: PetscFree(eps->data);
242: PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiSetDelayed_C",NULL);
243: PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiGetDelayed_C",NULL);
244: return(0);
245: }
247: PetscErrorCode EPSView_Arnoldi(EPS eps,PetscViewer viewer)
248: {
250: PetscBool isascii;
251: EPS_ARNOLDI *arnoldi = (EPS_ARNOLDI*)eps->data;
254: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
255: if (isascii && arnoldi->delayed) {
256: PetscViewerASCIIPrintf(viewer," using delayed reorthogonalization\n");
257: }
258: return(0);
259: }
261: SLEPC_EXTERN PetscErrorCode EPSCreate_Arnoldi(EPS eps)
262: {
263: EPS_ARNOLDI *ctx;
267: PetscNewLog(eps,&ctx);
268: eps->data = (void*)ctx;
270: eps->useds = PETSC_TRUE;
272: eps->ops->solve = EPSSolve_Arnoldi;
273: eps->ops->setup = EPSSetUp_Arnoldi;
274: eps->ops->setfromoptions = EPSSetFromOptions_Arnoldi;
275: eps->ops->destroy = EPSDestroy_Arnoldi;
276: eps->ops->view = EPSView_Arnoldi;
277: eps->ops->backtransform = EPSBackTransform_Default;
278: eps->ops->computevectors = EPSComputeVectors_Schur;
280: PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiSetDelayed_C",EPSArnoldiSetDelayed_Arnoldi);
281: PetscObjectComposeFunction((PetscObject)eps,"EPSArnoldiGetDelayed_C",EPSArnoldiGetDelayed_Arnoldi);
282: return(0);
283: }