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 : NEP routines related to the solution process
12 :
13 : References:
14 :
15 : [1] C. Campos and J.E. Roman, "NEP: a module for the parallel solution
16 : of nonlinear eigenvalue problems in SLEPc", ACM Trans. Math. Soft.
17 : 47(3), 23:1--23:29, 2021.
18 : */
19 :
20 : #include <slepc/private/nepimpl.h> /*I "slepcnep.h" I*/
21 : #include <slepc/private/bvimpl.h>
22 : #include <petscdraw.h>
23 :
24 : static PetscBool cited = PETSC_FALSE;
25 : static const char citation[] =
26 : "@Article{slepc-nep,\n"
27 : " author = \"C. Campos and J. E. Roman\",\n"
28 : " title = \"{NEP}: a module for the parallel solution of nonlinear eigenvalue problems in {SLEPc}\",\n"
29 : " journal = \"{ACM} Trans. Math. Software\",\n"
30 : " volume = \"47\",\n"
31 : " number = \"3\",\n"
32 : " pages = \"23:1--23:29\",\n"
33 : " year = \"2021\",\n"
34 : " doi = \"10.1145/3447544\"\n"
35 : "}\n";
36 :
37 768 : PetscErrorCode NEPComputeVectors(NEP nep)
38 : {
39 768 : PetscFunctionBegin;
40 768 : NEPCheckSolved(nep,1);
41 768 : if (nep->state==NEP_STATE_SOLVED) PetscTryTypeMethod(nep,computevectors);
42 768 : nep->state = NEP_STATE_EIGENVECTORS;
43 768 : PetscFunctionReturn(PETSC_SUCCESS);
44 : }
45 :
46 : /*@
47 : NEPSolve - Solves the nonlinear eigensystem.
48 :
49 : Collective
50 :
51 : Input Parameter:
52 : . nep - eigensolver context obtained from NEPCreate()
53 :
54 : Options Database Keys:
55 : + -nep_view - print information about the solver used
56 : . -nep_view_vectors - view the computed eigenvectors
57 : . -nep_view_values - view the computed eigenvalues
58 : . -nep_converged_reason - print reason for convergence, and number of iterations
59 : . -nep_error_absolute - print absolute errors of each eigenpair
60 : . -nep_error_relative - print relative errors of each eigenpair
61 : - -nep_error_backward - print backward errors of each eigenpair
62 :
63 : Notes:
64 : All the command-line options listed above admit an optional argument specifying
65 : the viewer type and options. For instance, use '-nep_view_vectors binary:myvecs.bin'
66 : to save the eigenvectors to a binary file, '-nep_view_values draw' to draw the computed
67 : eigenvalues graphically, or '-nep_error_relative :myerr.m:ascii_matlab' to save
68 : the errors in a file that can be executed in Matlab.
69 :
70 : Level: beginner
71 :
72 : .seealso: NEPCreate(), NEPSetUp(), NEPDestroy(), NEPSetTolerances()
73 : @*/
74 163 : PetscErrorCode NEPSolve(NEP nep)
75 : {
76 163 : PetscInt i;
77 :
78 163 : PetscFunctionBegin;
79 163 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
80 163 : if (nep->state>=NEP_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
81 163 : PetscCall(PetscCitationsRegister(citation,&cited));
82 163 : PetscCall(PetscLogEventBegin(NEP_Solve,nep,0,0,0));
83 :
84 : /* call setup */
85 163 : PetscCall(NEPSetUp(nep));
86 163 : nep->nconv = 0;
87 163 : nep->its = 0;
88 11757 : for (i=0;i<nep->ncv;i++) {
89 11594 : nep->eigr[i] = 0.0;
90 11594 : nep->eigi[i] = 0.0;
91 11594 : nep->errest[i] = 0.0;
92 11594 : nep->perm[i] = i;
93 : }
94 163 : PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view_pre"));
95 163 : PetscCall(RGViewFromOptions(nep->rg,NULL,"-rg_view"));
96 :
97 : /* call solver */
98 163 : PetscUseTypeMethod(nep,solve);
99 163 : PetscCheck(nep->reason,PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
100 163 : nep->state = NEP_STATE_SOLVED;
101 :
102 : /* Only the first nconv columns contain useful information */
103 163 : PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
104 163 : if (nep->twosided) PetscCall(BVSetActiveColumns(nep->W,0,nep->nconv));
105 :
106 163 : if (nep->refine==NEP_REFINE_SIMPLE && nep->rits>0 && nep->nconv>0) {
107 13 : PetscCall(NEPComputeVectors(nep));
108 13 : PetscCall(NEPNewtonRefinementSimple(nep,&nep->rits,nep->rtol,nep->nconv));
109 13 : nep->state = NEP_STATE_EIGENVECTORS;
110 : }
111 :
112 : /* sort eigenvalues according to nep->which parameter */
113 163 : PetscCall(SlepcSortEigenvalues(nep->sc,nep->nconv,nep->eigr,nep->eigi,nep->perm));
114 163 : PetscCall(PetscLogEventEnd(NEP_Solve,nep,0,0,0));
115 :
116 : /* various viewers */
117 163 : PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view"));
118 163 : PetscCall(NEPConvergedReasonViewFromOptions(nep));
119 163 : PetscCall(NEPErrorViewFromOptions(nep));
120 163 : PetscCall(NEPValuesViewFromOptions(nep));
121 163 : PetscCall(NEPVectorsViewFromOptions(nep));
122 :
123 : /* Remove the initial subspace */
124 163 : nep->nini = 0;
125 :
126 : /* Reset resolvent information */
127 163 : PetscCall(MatDestroy(&nep->resolvent));
128 163 : PetscFunctionReturn(PETSC_SUCCESS);
129 : }
130 :
131 : /*@
132 : NEPProjectOperator - Computes the projection of the nonlinear operator.
133 :
134 : Collective
135 :
136 : Input Parameters:
137 : + nep - the nonlinear eigensolver context
138 : . j0 - initial index
139 : - j1 - final index
140 :
141 : Notes:
142 : This is available for split operator only.
143 :
144 : The nonlinear operator T(lambda) is projected onto span(V), where V is
145 : an orthonormal basis built internally by the solver. The projected
146 : operator is equal to sum_i V'*A_i*V*f_i(lambda), so this function
147 : computes all matrices Ei = V'*A_i*V, and stores them in the extra
148 : matrices inside DS. Only rows/columns in the range [j0,j1-1] are computed,
149 : the previous ones are assumed to be available already.
150 :
151 : Level: developer
152 :
153 : .seealso: NEPSetSplitOperator()
154 : @*/
155 1 : PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)
156 : {
157 1 : PetscInt k;
158 1 : Mat G;
159 :
160 1 : PetscFunctionBegin;
161 1 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
162 4 : PetscValidLogicalCollectiveInt(nep,j0,2);
163 4 : PetscValidLogicalCollectiveInt(nep,j1,3);
164 1 : NEPCheckProblem(nep,1);
165 1 : NEPCheckSplit(nep,1);
166 1 : PetscCall(BVSetActiveColumns(nep->V,j0,j1));
167 4 : for (k=0;k<nep->nt;k++) {
168 3 : PetscCall(DSGetMat(nep->ds,DSMatExtra[k],&G));
169 3 : PetscCall(BVMatProject(nep->V,nep->A[k],nep->V,G));
170 3 : PetscCall(DSRestoreMat(nep->ds,DSMatExtra[k],&G));
171 : }
172 1 : PetscFunctionReturn(PETSC_SUCCESS);
173 : }
174 :
175 : /*@
176 : NEPApplyFunction - Applies the nonlinear function T(lambda) to a given vector.
177 :
178 : Collective
179 :
180 : Input Parameters:
181 : + nep - the nonlinear eigensolver context
182 : . lambda - scalar argument
183 : . x - vector to be multiplied against
184 : - v - workspace vector (used only in the case of split form)
185 :
186 : Output Parameters:
187 : + y - result vector
188 : . A - (optional) Function matrix, for callback interface only
189 : - B - (unused) preconditioning matrix
190 :
191 : Note:
192 : If the nonlinear operator is represented in split form, the result
193 : y = T(lambda)*x is computed without building T(lambda) explicitly. In
194 : that case, parameters A and B are not used. Otherwise, the matrix
195 : T(lambda) is built and the effect is the same as a call to
196 : NEPComputeFunction() followed by a MatMult().
197 :
198 : Level: developer
199 :
200 : .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyAdjoint()
201 : @*/
202 862 : PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
203 : {
204 862 : PetscInt i;
205 862 : PetscScalar alpha;
206 :
207 862 : PetscFunctionBegin;
208 862 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
209 3448 : PetscValidLogicalCollectiveScalar(nep,lambda,2);
210 862 : PetscValidHeaderSpecific(x,VEC_CLASSID,3);
211 862 : if (v) PetscValidHeaderSpecific(v,VEC_CLASSID,4);
212 862 : PetscValidHeaderSpecific(y,VEC_CLASSID,5);
213 862 : if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,6);
214 862 : if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,7);
215 :
216 862 : if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
217 706 : PetscCall(VecSet(y,0.0));
218 2776 : for (i=0;i<nep->nt;i++) {
219 2070 : PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
220 2070 : PetscCall(MatMult(nep->A[i],x,v));
221 2070 : PetscCall(VecAXPY(y,alpha,v));
222 : }
223 : } else {
224 156 : if (!A) A = nep->function;
225 156 : PetscCall(NEPComputeFunction(nep,lambda,A,A));
226 156 : PetscCall(MatMult(A,x,y));
227 : }
228 862 : PetscFunctionReturn(PETSC_SUCCESS);
229 : }
230 :
231 : /*@
232 : NEPApplyAdjoint - Applies the adjoint nonlinear function T(lambda)^* to a given vector.
233 :
234 : Collective
235 :
236 : Input Parameters:
237 : + nep - the nonlinear eigensolver context
238 : . lambda - scalar argument
239 : . x - vector to be multiplied against
240 : - v - workspace vector (used only in the case of split form)
241 :
242 : Output Parameters:
243 : + y - result vector
244 : . A - (optional) Function matrix, for callback interface only
245 : - B - (unused) preconditioning matrix
246 :
247 : Level: developer
248 :
249 : .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyFunction()
250 : @*/
251 21 : PetscErrorCode NEPApplyAdjoint(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
252 : {
253 21 : PetscInt i;
254 21 : PetscScalar alpha;
255 21 : Vec w;
256 :
257 21 : PetscFunctionBegin;
258 21 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
259 84 : PetscValidLogicalCollectiveScalar(nep,lambda,2);
260 21 : PetscValidHeaderSpecific(x,VEC_CLASSID,3);
261 21 : if (v) PetscValidHeaderSpecific(v,VEC_CLASSID,4);
262 21 : PetscValidHeaderSpecific(y,VEC_CLASSID,5);
263 21 : if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,6);
264 21 : if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,7);
265 :
266 21 : PetscCall(VecDuplicate(x,&w));
267 21 : PetscCall(VecCopy(x,w));
268 21 : PetscCall(VecConjugate(w));
269 21 : if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
270 13 : PetscCall(VecSet(y,0.0));
271 52 : for (i=0;i<nep->nt;i++) {
272 39 : PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
273 39 : PetscCall(MatMultTranspose(nep->A[i],w,v));
274 39 : PetscCall(VecAXPY(y,alpha,v));
275 : }
276 : } else {
277 8 : if (!A) A = nep->function;
278 8 : PetscCall(NEPComputeFunction(nep,lambda,A,A));
279 8 : PetscCall(MatMultTranspose(A,w,y));
280 : }
281 21 : PetscCall(VecDestroy(&w));
282 21 : PetscCall(VecConjugate(y));
283 21 : PetscFunctionReturn(PETSC_SUCCESS);
284 : }
285 :
286 : /*@
287 : NEPApplyJacobian - Applies the nonlinear Jacobian T'(lambda) to a given vector.
288 :
289 : Collective
290 :
291 : Input Parameters:
292 : + nep - the nonlinear eigensolver context
293 : . lambda - scalar argument
294 : . x - vector to be multiplied against
295 : - v - workspace vector (used only in the case of split form)
296 :
297 : Output Parameters:
298 : + y - result vector
299 : - A - (optional) Jacobian matrix, for callback interface only
300 :
301 : Note:
302 : If the nonlinear operator is represented in split form, the result
303 : y = T'(lambda)*x is computed without building T'(lambda) explicitly. In
304 : that case, parameter A is not used. Otherwise, the matrix
305 : T'(lambda) is built and the effect is the same as a call to
306 : NEPComputeJacobian() followed by a MatMult().
307 :
308 : Level: developer
309 :
310 : .seealso: NEPSetSplitOperator(), NEPComputeJacobian()
311 : @*/
312 32 : PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)
313 : {
314 32 : PetscInt i;
315 32 : PetscScalar alpha;
316 :
317 32 : PetscFunctionBegin;
318 32 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
319 128 : PetscValidLogicalCollectiveScalar(nep,lambda,2);
320 32 : PetscValidHeaderSpecific(x,VEC_CLASSID,3);
321 32 : if (v) PetscValidHeaderSpecific(v,VEC_CLASSID,4);
322 32 : PetscValidHeaderSpecific(y,VEC_CLASSID,5);
323 32 : if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,6);
324 :
325 32 : if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
326 32 : PetscCall(VecSet(y,0.0));
327 128 : for (i=0;i<nep->nt;i++) {
328 96 : PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
329 96 : PetscCall(MatMult(nep->A[i],x,v));
330 96 : PetscCall(VecAXPY(y,alpha,v));
331 : }
332 : } else {
333 0 : if (!A) A = nep->jacobian;
334 0 : PetscCall(NEPComputeJacobian(nep,lambda,A));
335 0 : PetscCall(MatMult(A,x,y));
336 : }
337 32 : PetscFunctionReturn(PETSC_SUCCESS);
338 : }
339 :
340 : /*@
341 : NEPGetIterationNumber - Gets the current iteration number. If the
342 : call to NEPSolve() is complete, then it returns the number of iterations
343 : carried out by the solution method.
344 :
345 : Not Collective
346 :
347 : Input Parameter:
348 : . nep - the nonlinear eigensolver context
349 :
350 : Output Parameter:
351 : . its - number of iterations
352 :
353 : Note:
354 : During the i-th iteration this call returns i-1. If NEPSolve() is
355 : complete, then parameter "its" contains either the iteration number at
356 : which convergence was successfully reached, or failure was detected.
357 : Call NEPGetConvergedReason() to determine if the solver converged or
358 : failed and why.
359 :
360 : Level: intermediate
361 :
362 : .seealso: NEPGetConvergedReason(), NEPSetTolerances()
363 : @*/
364 9 : PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)
365 : {
366 9 : PetscFunctionBegin;
367 9 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
368 9 : PetscAssertPointer(its,2);
369 9 : *its = nep->its;
370 9 : PetscFunctionReturn(PETSC_SUCCESS);
371 : }
372 :
373 : /*@
374 : NEPGetConverged - Gets the number of converged eigenpairs.
375 :
376 : Not Collective
377 :
378 : Input Parameter:
379 : . nep - the nonlinear eigensolver context
380 :
381 : Output Parameter:
382 : . nconv - number of converged eigenpairs
383 :
384 : Note:
385 : This function should be called after NEPSolve() has finished.
386 :
387 : Level: beginner
388 :
389 : .seealso: NEPSetDimensions(), NEPSolve(), NEPGetEigenpair()
390 : @*/
391 10 : PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)
392 : {
393 10 : PetscFunctionBegin;
394 10 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
395 10 : PetscAssertPointer(nconv,2);
396 10 : NEPCheckSolved(nep,1);
397 10 : *nconv = nep->nconv;
398 10 : PetscFunctionReturn(PETSC_SUCCESS);
399 : }
400 :
401 : /*@
402 : NEPGetConvergedReason - Gets the reason why the NEPSolve() iteration was
403 : stopped.
404 :
405 : Not Collective
406 :
407 : Input Parameter:
408 : . nep - the nonlinear eigensolver context
409 :
410 : Output Parameter:
411 : . reason - negative value indicates diverged, positive value converged
412 :
413 : Options Database Key:
414 : . -nep_converged_reason - print the reason to a viewer
415 :
416 : Notes:
417 : Possible values for reason are
418 : + NEP_CONVERGED_TOL - converged up to tolerance
419 : . NEP_CONVERGED_USER - converged due to a user-defined condition
420 : . NEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
421 : . NEP_DIVERGED_BREAKDOWN - generic breakdown in method
422 : . NEP_DIVERGED_LINEAR_SOLVE - inner linear solve failed
423 : - NEP_DIVERGED_SUBSPACE_EXHAUSTED - run out of space for the basis in an
424 : unrestarted solver
425 :
426 : Can only be called after the call to NEPSolve() is complete.
427 :
428 : Level: intermediate
429 :
430 : .seealso: NEPSetTolerances(), NEPSolve(), NEPConvergedReason
431 : @*/
432 1 : PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)
433 : {
434 1 : PetscFunctionBegin;
435 1 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
436 1 : PetscAssertPointer(reason,2);
437 1 : NEPCheckSolved(nep,1);
438 1 : *reason = nep->reason;
439 1 : PetscFunctionReturn(PETSC_SUCCESS);
440 : }
441 :
442 : /*@C
443 : NEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
444 : NEPSolve(). The solution consists in both the eigenvalue and the eigenvector.
445 :
446 : Collective
447 :
448 : Input Parameters:
449 : + nep - nonlinear eigensolver context
450 : - i - index of the solution
451 :
452 : Output Parameters:
453 : + eigr - real part of eigenvalue
454 : . eigi - imaginary part of eigenvalue
455 : . Vr - real part of eigenvector
456 : - Vi - imaginary part of eigenvector
457 :
458 : Notes:
459 : It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
460 : required. Otherwise, the caller must provide valid Vec objects, i.e.,
461 : they must be created by the calling program with e.g. MatCreateVecs().
462 :
463 : If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
464 : configured with complex scalars the eigenvalue is stored
465 : directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
466 : set to zero). In any case, the user can pass NULL in Vr or Vi if one of
467 : them is not required.
468 :
469 : The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
470 : Eigenpairs are indexed according to the ordering criterion established
471 : with NEPSetWhichEigenpairs().
472 :
473 : Level: beginner
474 :
475 : .seealso: NEPSolve(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPGetLeftEigenvector()
476 : @*/
477 727 : PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
478 : {
479 727 : PetscInt k;
480 :
481 727 : PetscFunctionBegin;
482 727 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
483 2908 : PetscValidLogicalCollectiveInt(nep,i,2);
484 727 : if (Vr) { PetscValidHeaderSpecific(Vr,VEC_CLASSID,5); PetscCheckSameComm(nep,1,Vr,5); }
485 727 : if (Vi) { PetscValidHeaderSpecific(Vi,VEC_CLASSID,6); PetscCheckSameComm(nep,1,Vi,6); }
486 727 : NEPCheckSolved(nep,1);
487 727 : PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
488 727 : PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
489 :
490 727 : PetscCall(NEPComputeVectors(nep));
491 727 : k = nep->perm[i];
492 :
493 : /* eigenvalue */
494 : #if defined(PETSC_USE_COMPLEX)
495 727 : if (eigr) *eigr = nep->eigr[k];
496 727 : if (eigi) *eigi = 0;
497 : #else
498 : if (eigr) *eigr = nep->eigr[k];
499 : if (eigi) *eigi = nep->eigi[k];
500 : #endif
501 :
502 : /* eigenvector */
503 727 : PetscCall(BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi));
504 727 : PetscFunctionReturn(PETSC_SUCCESS);
505 : }
506 :
507 : /*@
508 : NEPGetLeftEigenvector - Gets the i-th left eigenvector as computed by NEPSolve().
509 :
510 : Collective
511 :
512 : Input Parameters:
513 : + nep - eigensolver context
514 : - i - index of the solution
515 :
516 : Output Parameters:
517 : + Wr - real part of left eigenvector
518 : - Wi - imaginary part of left eigenvector
519 :
520 : Notes:
521 : The caller must provide valid Vec objects, i.e., they must be created
522 : by the calling program with e.g. MatCreateVecs().
523 :
524 : If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
525 : configured with complex scalars the eigenvector is stored directly in Wr
526 : (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if one of
527 : them is not required.
528 :
529 : The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
530 : Eigensolutions are indexed according to the ordering criterion established
531 : with NEPSetWhichEigenpairs().
532 :
533 : Left eigenvectors are available only if the twosided flag was set, see
534 : NEPSetTwoSided().
535 :
536 : Level: intermediate
537 :
538 : .seealso: NEPGetEigenpair(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPSetTwoSided()
539 : @*/
540 21 : PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)
541 : {
542 21 : PetscInt k;
543 :
544 21 : PetscFunctionBegin;
545 21 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
546 84 : PetscValidLogicalCollectiveInt(nep,i,2);
547 21 : if (Wr) { PetscValidHeaderSpecific(Wr,VEC_CLASSID,3); PetscCheckSameComm(nep,1,Wr,3); }
548 21 : if (Wi) { PetscValidHeaderSpecific(Wi,VEC_CLASSID,4); PetscCheckSameComm(nep,1,Wi,4); }
549 21 : NEPCheckSolved(nep,1);
550 21 : PetscCheck(nep->twosided,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with NEPSetTwoSided");
551 21 : PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
552 21 : PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
553 21 : PetscCall(NEPComputeVectors(nep));
554 21 : k = nep->perm[i];
555 21 : PetscCall(BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi));
556 21 : PetscFunctionReturn(PETSC_SUCCESS);
557 : }
558 :
559 : /*@
560 : NEPGetErrorEstimate - Returns the error estimate associated to the i-th
561 : computed eigenpair.
562 :
563 : Not Collective
564 :
565 : Input Parameters:
566 : + nep - nonlinear eigensolver context
567 : - i - index of eigenpair
568 :
569 : Output Parameter:
570 : . errest - the error estimate
571 :
572 : Notes:
573 : This is the error estimate used internally by the eigensolver. The actual
574 : error bound can be computed with NEPComputeError().
575 :
576 : Level: advanced
577 :
578 : .seealso: NEPComputeError()
579 : @*/
580 4 : PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)
581 : {
582 4 : PetscFunctionBegin;
583 4 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
584 4 : PetscAssertPointer(errest,3);
585 4 : NEPCheckSolved(nep,1);
586 4 : PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
587 4 : PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
588 4 : *errest = nep->errest[nep->perm[i]];
589 4 : PetscFunctionReturn(PETSC_SUCCESS);
590 : }
591 :
592 : /*
593 : NEPComputeResidualNorm_Private - Computes the norm of the residual vector
594 : associated with an eigenpair.
595 :
596 : Input Parameters:
597 : adj - whether the adjoint T^* must be used instead of T
598 : lambda - eigenvalue
599 : x - eigenvector
600 : w - array of work vectors (two vectors in split form, one vector otherwise)
601 : */
602 875 : PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)
603 : {
604 875 : Vec y,z=NULL;
605 :
606 875 : PetscFunctionBegin;
607 875 : y = w[0];
608 875 : if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
609 875 : if (adj) PetscCall(NEPApplyAdjoint(nep,lambda,x,z,y,NULL,NULL));
610 862 : else PetscCall(NEPApplyFunction(nep,lambda,x,z,y,NULL,NULL));
611 875 : PetscCall(VecNorm(y,NORM_2,norm));
612 875 : PetscFunctionReturn(PETSC_SUCCESS);
613 : }
614 :
615 : /*@
616 : NEPComputeError - Computes the error (based on the residual norm) associated
617 : with the i-th computed eigenpair.
618 :
619 : Collective
620 :
621 : Input Parameters:
622 : + nep - the nonlinear eigensolver context
623 : . i - the solution index
624 : - type - the type of error to compute
625 :
626 : Output Parameter:
627 : . error - the error
628 :
629 : Notes:
630 : The error can be computed in various ways, all of them based on the residual
631 : norm computed as ||T(lambda)x||_2 where lambda is the eigenvalue and x is the
632 : eigenvector.
633 :
634 : Level: beginner
635 :
636 : .seealso: NEPErrorType, NEPSolve(), NEPGetErrorEstimate()
637 : @*/
638 695 : PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)
639 : {
640 695 : Vec xr,xi=NULL;
641 695 : PetscInt j,nwork,issplit=0;
642 695 : PetscScalar kr,ki,s;
643 695 : PetscReal er,z=0.0,errorl,nrm;
644 695 : PetscBool flg;
645 :
646 695 : PetscFunctionBegin;
647 695 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
648 2780 : PetscValidLogicalCollectiveInt(nep,i,2);
649 2780 : PetscValidLogicalCollectiveEnum(nep,type,3);
650 695 : PetscAssertPointer(error,4);
651 695 : NEPCheckSolved(nep,1);
652 :
653 : /* allocate work vectors */
654 : #if defined(PETSC_USE_COMPLEX)
655 695 : nwork = 2;
656 : #else
657 : nwork = 3;
658 : #endif
659 695 : if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
660 587 : issplit = 1;
661 587 : nwork++; /* need an extra work vector for NEPComputeResidualNorm_Private */
662 : }
663 695 : PetscCall(NEPSetWorkVecs(nep,nwork));
664 695 : xr = nep->work[issplit+1];
665 : #if !defined(PETSC_USE_COMPLEX)
666 : xi = nep->work[issplit+2];
667 : #endif
668 :
669 : /* compute residual norms */
670 695 : PetscCall(NEPGetEigenpair(nep,i,&kr,&ki,xr,xi));
671 : #if !defined(PETSC_USE_COMPLEX)
672 : PetscCheck(ki==0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Not implemented for complex eigenvalues with real scalars");
673 : #endif
674 695 : PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error));
675 695 : PetscCall(VecNorm(xr,NORM_2,&er));
676 :
677 : /* if two-sided, compute left residual norm and take the maximum */
678 695 : if (nep->twosided) {
679 13 : PetscCall(NEPGetLeftEigenvector(nep,i,xr,xi));
680 13 : PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_TRUE,kr,xr,nep->work,&errorl));
681 24 : *error = PetscMax(*error,errorl);
682 : }
683 :
684 : /* compute error */
685 695 : switch (type) {
686 : case NEP_ERROR_ABSOLUTE:
687 : break;
688 602 : case NEP_ERROR_RELATIVE:
689 602 : *error /= PetscAbsScalar(kr)*er;
690 602 : break;
691 92 : case NEP_ERROR_BACKWARD:
692 92 : if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
693 42 : PetscCall(NEPComputeFunction(nep,kr,nep->function,nep->function));
694 42 : PetscCall(MatHasOperation(nep->function,MATOP_NORM,&flg));
695 42 : PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
696 42 : PetscCall(MatNorm(nep->function,NORM_INFINITY,&nrm));
697 42 : *error /= nrm*er;
698 42 : break;
699 : }
700 : /* initialization of matrix norms */
701 50 : if (!nep->nrma[0]) {
702 46 : for (j=0;j<nep->nt;j++) {
703 32 : PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&flg));
704 32 : PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
705 32 : PetscCall(MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]));
706 : }
707 : }
708 161 : for (j=0;j<nep->nt;j++) {
709 111 : PetscCall(FNEvaluateFunction(nep->f[j],kr,&s));
710 111 : z = z + nep->nrma[j]*PetscAbsScalar(s);
711 : }
712 50 : *error /= z*er;
713 50 : break;
714 0 : default:
715 0 : SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
716 : }
717 695 : PetscFunctionReturn(PETSC_SUCCESS);
718 : }
719 :
720 : /*@
721 : NEPComputeFunction - Computes the function matrix T(lambda) that has been
722 : set with NEPSetFunction().
723 :
724 : Collective
725 :
726 : Input Parameters:
727 : + nep - the NEP context
728 : - lambda - the scalar argument
729 :
730 : Output Parameters:
731 : + A - Function matrix
732 : - B - optional preconditioning matrix
733 :
734 : Notes:
735 : NEPComputeFunction() is typically used within nonlinear eigensolvers
736 : implementations, so most users would not generally call this routine
737 : themselves.
738 :
739 : Level: developer
740 :
741 : .seealso: NEPSetFunction(), NEPGetFunction()
742 : @*/
743 3418 : PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)
744 : {
745 3418 : PetscInt i;
746 3418 : PetscScalar alpha;
747 :
748 3418 : PetscFunctionBegin;
749 3418 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
750 3418 : NEPCheckProblem(nep,1);
751 3418 : switch (nep->fui) {
752 2040 : case NEP_USER_INTERFACE_CALLBACK:
753 2040 : PetscCheck(nep->computefunction,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
754 2040 : PetscCall(PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0));
755 2040 : PetscCallBack("NEP user Function function",(*nep->computefunction)(nep,lambda,A,B,nep->functionctx));
756 2040 : PetscCall(PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0));
757 : break;
758 1378 : case NEP_USER_INTERFACE_SPLIT:
759 1378 : PetscCall(MatZeroEntries(A));
760 1378 : if (A != B) PetscCall(MatZeroEntries(B));
761 5442 : for (i=0;i<nep->nt;i++) {
762 4064 : PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
763 4064 : PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
764 4064 : if (A != B) PetscCall(MatAXPY(B,alpha,nep->P[i],nep->mstrp));
765 : }
766 : break;
767 : }
768 3418 : PetscFunctionReturn(PETSC_SUCCESS);
769 : }
770 :
771 : /*@
772 : NEPComputeJacobian - Computes the Jacobian matrix T'(lambda) that has been
773 : set with NEPSetJacobian().
774 :
775 : Collective
776 :
777 : Input Parameters:
778 : + nep - the NEP context
779 : - lambda - the scalar argument
780 :
781 : Output Parameters:
782 : . A - Jacobian matrix
783 :
784 : Notes:
785 : Most users should not need to explicitly call this routine, as it
786 : is used internally within the nonlinear eigensolvers.
787 :
788 : Level: developer
789 :
790 : .seealso: NEPSetJacobian(), NEPGetJacobian()
791 : @*/
792 1743 : PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)
793 : {
794 1743 : PetscInt i;
795 1743 : PetscScalar alpha;
796 :
797 1743 : PetscFunctionBegin;
798 1743 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
799 1743 : NEPCheckProblem(nep,1);
800 1743 : switch (nep->fui) {
801 641 : case NEP_USER_INTERFACE_CALLBACK:
802 641 : PetscCheck(nep->computejacobian,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
803 641 : PetscCall(PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0));
804 641 : PetscCallBack("NEP user Jacobian function",(*nep->computejacobian)(nep,lambda,A,nep->jacobianctx));
805 641 : PetscCall(PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0));
806 : break;
807 1102 : case NEP_USER_INTERFACE_SPLIT:
808 1102 : PetscCall(MatZeroEntries(A));
809 4340 : for (i=0;i<nep->nt;i++) {
810 3238 : PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
811 3238 : PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
812 : }
813 : break;
814 : }
815 1743 : PetscFunctionReturn(PETSC_SUCCESS);
816 : }
|
