Actual source code: nepdefault.c
slepc-3.21.0 2024-03-30
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: Simple default routines for common NEP operations
12: */
14: #include <slepc/private/nepimpl.h>
16: /*@
17: NEPSetWorkVecs - Sets a number of work vectors into a NEP object
19: Collective
21: Input Parameters:
22: + nep - nonlinear eigensolver context
23: - nw - number of work vectors to allocate
25: Developer Notes:
26: This is SLEPC_EXTERN because it may be required by user plugin NEP
27: implementations.
29: Level: developer
31: .seealso: NEPSetUp()
32: @*/
33: PetscErrorCode NEPSetWorkVecs(NEP nep,PetscInt nw)
34: {
35: Vec t;
37: PetscFunctionBegin;
40: PetscCheck(nw>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %" PetscInt_FMT,nw);
41: if (nep->nwork < nw) {
42: PetscCall(VecDestroyVecs(nep->nwork,&nep->work));
43: nep->nwork = nw;
44: PetscCall(BVGetColumn(nep->V,0,&t));
45: PetscCall(VecDuplicateVecs(t,nw,&nep->work));
46: PetscCall(BVRestoreColumn(nep->V,0,&t));
47: }
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: /*
52: NEPGetDefaultShift - Return the value of sigma to start the nonlinear iteration.
53: */
54: PetscErrorCode NEPGetDefaultShift(NEP nep,PetscScalar *sigma)
55: {
56: PetscFunctionBegin;
57: PetscAssertPointer(sigma,2);
58: switch (nep->which) {
59: case NEP_LARGEST_MAGNITUDE:
60: case NEP_LARGEST_IMAGINARY:
61: case NEP_ALL:
62: case NEP_WHICH_USER:
63: *sigma = 1.0; /* arbitrary value */
64: break;
65: case NEP_SMALLEST_MAGNITUDE:
66: case NEP_SMALLEST_IMAGINARY:
67: *sigma = 0.0;
68: break;
69: case NEP_LARGEST_REAL:
70: *sigma = PETSC_MAX_REAL;
71: break;
72: case NEP_SMALLEST_REAL:
73: *sigma = PETSC_MIN_REAL;
74: break;
75: case NEP_TARGET_MAGNITUDE:
76: case NEP_TARGET_REAL:
77: case NEP_TARGET_IMAGINARY:
78: *sigma = nep->target;
79: break;
80: }
81: PetscFunctionReturn(PETSC_SUCCESS);
82: }
84: /*
85: NEPConvergedRelative - Checks convergence relative to the eigenvalue.
86: */
87: PetscErrorCode NEPConvergedRelative(NEP nep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
88: {
89: PetscReal w;
91: PetscFunctionBegin;
92: w = SlepcAbsEigenvalue(eigr,eigi);
93: *errest = res/w;
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: /*
98: NEPConvergedAbsolute - Checks convergence absolutely.
99: */
100: PetscErrorCode NEPConvergedAbsolute(NEP nep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
101: {
102: PetscFunctionBegin;
103: *errest = res;
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: /*
108: NEPConvergedNorm - Checks convergence relative to the matrix norms.
109: */
110: PetscErrorCode NEPConvergedNorm(NEP nep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
111: {
112: PetscScalar s;
113: PetscReal w=0.0;
114: PetscInt j;
115: PetscBool flg;
117: PetscFunctionBegin;
118: if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
119: PetscCall(NEPComputeFunction(nep,eigr,nep->function,nep->function));
120: PetscCall(MatHasOperation(nep->function,MATOP_NORM,&flg));
121: PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
122: PetscCall(MatNorm(nep->function,NORM_INFINITY,&w));
123: } else {
124: /* initialization of matrix norms */
125: if (!nep->nrma[0]) {
126: for (j=0;j<nep->nt;j++) {
127: PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&flg));
128: PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The convergence test related to the matrix norms requires a matrix norm operation");
129: PetscCall(MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]));
130: }
131: }
132: for (j=0;j<nep->nt;j++) {
133: PetscCall(FNEvaluateFunction(nep->f[j],eigr,&s));
134: w = w + nep->nrma[j]*PetscAbsScalar(s);
135: }
136: }
137: *errest = res/w;
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: /*@C
142: NEPStoppingBasic - Default routine to determine whether the outer eigensolver
143: iteration must be stopped.
145: Collective
147: Input Parameters:
148: + nep - nonlinear eigensolver context obtained from NEPCreate()
149: . its - current number of iterations
150: . max_it - maximum number of iterations
151: . nconv - number of currently converged eigenpairs
152: . nev - number of requested eigenpairs
153: - ctx - context (not used here)
155: Output Parameter:
156: . reason - result of the stopping test
158: Notes:
159: A positive value of reason indicates that the iteration has finished successfully
160: (converged), and a negative value indicates an error condition (diverged). If
161: the iteration needs to be continued, reason must be set to NEP_CONVERGED_ITERATING
162: (zero).
164: NEPStoppingBasic() will stop if all requested eigenvalues are converged, or if
165: the maximum number of iterations has been reached.
167: Use NEPSetStoppingTest() to provide your own test instead of using this one.
169: Level: advanced
171: .seealso: NEPSetStoppingTest(), NEPConvergedReason, NEPGetConvergedReason()
172: @*/
173: PetscErrorCode NEPStoppingBasic(NEP nep,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,NEPConvergedReason *reason,void *ctx)
174: {
175: PetscFunctionBegin;
176: *reason = NEP_CONVERGED_ITERATING;
177: if (nconv >= nev) {
178: PetscCall(PetscInfo(nep,"Nonlinear eigensolver finished successfully: %" PetscInt_FMT " eigenpairs converged at iteration %" PetscInt_FMT "\n",nconv,its));
179: *reason = NEP_CONVERGED_TOL;
180: } else if (its >= max_it) {
181: *reason = NEP_DIVERGED_ITS;
182: PetscCall(PetscInfo(nep,"Nonlinear eigensolver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
183: }
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: PetscErrorCode NEPComputeVectors_Schur(NEP nep)
188: {
189: Mat Z;
191: PetscFunctionBegin;
192: PetscCall(DSVectors(nep->ds,DS_MAT_X,NULL,NULL));
193: PetscCall(DSGetMat(nep->ds,DS_MAT_X,&Z));
194: PetscCall(BVMultInPlace(nep->V,Z,0,nep->nconv));
195: PetscCall(DSRestoreMat(nep->ds,DS_MAT_X,&Z));
196: PetscCall(BVNormalize(nep->V,nep->eigi));
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }