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