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 571 : PetscErrorCode NEPComputeVectors(NEP nep)
38 : {
39 571 : PetscFunctionBegin;
40 571 : NEPCheckSolved(nep,1);
41 571 : if (nep->state==NEP_STATE_SOLVED) PetscTryTypeMethod(nep,computevectors);
42 571 : nep->state = NEP_STATE_EIGENVECTORS;
43 571 : 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 120 : PetscErrorCode NEPSolve(NEP nep)
75 : {
76 120 : PetscInt i;
77 :
78 120 : PetscFunctionBegin;
79 120 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
80 120 : if (nep->state>=NEP_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
81 120 : PetscCall(PetscCitationsRegister(citation,&cited));
82 120 : PetscCall(PetscLogEventBegin(NEP_Solve,nep,0,0,0));
83 :
84 : /* call setup */
85 120 : PetscCall(NEPSetUp(nep));
86 120 : nep->nconv = 0;
87 120 : nep->its = 0;
88 1522 : for (i=0;i<nep->ncv;i++) {
89 1402 : nep->eigr[i] = 0.0;
90 1402 : nep->eigi[i] = 0.0;
91 1402 : nep->errest[i] = 0.0;
92 1402 : nep->perm[i] = i;
93 : }
94 120 : PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view_pre"));
95 120 : PetscCall(RGViewFromOptions(nep->rg,NULL,"-rg_view"));
96 :
97 : /* call solver */
98 120 : PetscUseTypeMethod(nep,solve);
99 120 : PetscCheck(nep->reason,PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
100 120 : nep->state = NEP_STATE_SOLVED;
101 :
102 : /* Only the first nconv columns contain useful information */
103 120 : PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
104 120 : if (nep->twosided) PetscCall(BVSetActiveColumns(nep->W,0,nep->nconv));
105 :
106 120 : 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 120 : PetscCall(SlepcSortEigenvalues(nep->sc,nep->nconv,nep->eigr,nep->eigi,nep->perm));
114 120 : PetscCall(PetscLogEventEnd(NEP_Solve,nep,0,0,0));
115 :
116 : /* various viewers */
117 120 : PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view"));
118 120 : PetscCall(NEPConvergedReasonViewFromOptions(nep));
119 120 : PetscCall(NEPErrorViewFromOptions(nep));
120 120 : PetscCall(NEPValuesViewFromOptions(nep));
121 120 : PetscCall(NEPVectorsViewFromOptions(nep));
122 :
123 : /* Remove the initial subspace */
124 120 : nep->nini = 0;
125 :
126 : /* Reset resolvent information */
127 120 : PetscCall(MatDestroy(&nep->resolvent));
128 120 : 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 499 : PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
203 : {
204 499 : PetscInt i;
205 499 : PetscScalar alpha;
206 :
207 499 : PetscFunctionBegin;
208 499 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
209 1996 : PetscValidLogicalCollectiveScalar(nep,lambda,2);
210 499 : PetscValidHeaderSpecific(x,VEC_CLASSID,3);
211 499 : if (v) PetscValidHeaderSpecific(v,VEC_CLASSID,4);
212 499 : PetscValidHeaderSpecific(y,VEC_CLASSID,5);
213 499 : if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,6);
214 499 : if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,7);
215 :
216 499 : if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
217 435 : PetscCall(VecSet(y,0.0));
218 1716 : for (i=0;i<nep->nt;i++) {
219 1281 : PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
220 1281 : PetscCall(MatMult(nep->A[i],x,v));
221 1281 : PetscCall(VecAXPY(y,alpha,v));
222 : }
223 : } else {
224 64 : if (!A) A = nep->function;
225 64 : PetscCall(NEPComputeFunction(nep,lambda,A,A));
226 64 : PetscCall(MatMult(A,x,y));
227 : }
228 499 : 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 8 : PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)
365 : {
366 8 : PetscFunctionBegin;
367 8 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
368 8 : PetscAssertPointer(its,2);
369 8 : *its = nep->its;
370 8 : 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 9 : PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)
392 : {
393 9 : PetscFunctionBegin;
394 9 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
395 9 : PetscAssertPointer(nconv,2);
396 9 : NEPCheckSolved(nep,1);
397 9 : *nconv = nep->nconv;
398 9 : 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 530 : PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
478 : {
479 530 : PetscInt k;
480 :
481 530 : PetscFunctionBegin;
482 530 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
483 2120 : PetscValidLogicalCollectiveInt(nep,i,2);
484 530 : if (Vr) { PetscValidHeaderSpecific(Vr,VEC_CLASSID,5); PetscCheckSameComm(nep,1,Vr,5); }
485 530 : if (Vi) { PetscValidHeaderSpecific(Vi,VEC_CLASSID,6); PetscCheckSameComm(nep,1,Vi,6); }
486 530 : NEPCheckSolved(nep,1);
487 530 : PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
488 530 : PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
489 :
490 530 : PetscCall(NEPComputeVectors(nep));
491 530 : k = nep->perm[i];
492 :
493 : /* eigenvalue */
494 : #if defined(PETSC_USE_COMPLEX)
495 : if (eigr) *eigr = nep->eigr[k];
496 : if (eigi) *eigi = 0;
497 : #else
498 530 : if (eigr) *eigr = nep->eigr[k];
499 530 : if (eigi) *eigi = nep->eigi[k];
500 : #endif
501 :
502 : /* eigenvector */
503 530 : PetscCall(BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi));
504 530 : 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 3 : PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)
581 : {
582 3 : PetscFunctionBegin;
583 3 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
584 3 : PetscAssertPointer(errest,3);
585 3 : NEPCheckSolved(nep,1);
586 3 : PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
587 3 : PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
588 3 : *errest = nep->errest[nep->perm[i]];
589 3 : 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 512 : PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)
603 : {
604 512 : Vec y,z=NULL;
605 :
606 512 : PetscFunctionBegin;
607 512 : y = w[0];
608 512 : if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
609 512 : if (adj) PetscCall(NEPApplyAdjoint(nep,lambda,x,z,y,NULL,NULL));
610 499 : else PetscCall(NEPApplyFunction(nep,lambda,x,z,y,NULL,NULL));
611 512 : PetscCall(VecNorm(y,NORM_2,norm));
612 512 : 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 499 : PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)
639 : {
640 499 : Vec xr,xi=NULL;
641 499 : PetscInt j,nwork,issplit=0;
642 499 : PetscScalar kr,ki,s;
643 499 : PetscReal er,z=0.0,errorl,nrm;
644 499 : PetscBool flg;
645 :
646 499 : PetscFunctionBegin;
647 499 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
648 1996 : PetscValidLogicalCollectiveInt(nep,i,2);
649 1996 : PetscValidLogicalCollectiveEnum(nep,type,3);
650 499 : PetscAssertPointer(error,4);
651 499 : NEPCheckSolved(nep,1);
652 :
653 : /* allocate work vectors */
654 : #if defined(PETSC_USE_COMPLEX)
655 : nwork = 2;
656 : #else
657 499 : nwork = 3;
658 : #endif
659 499 : if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
660 435 : issplit = 1;
661 435 : nwork++; /* need an extra work vector for NEPComputeResidualNorm_Private */
662 : }
663 499 : PetscCall(NEPSetWorkVecs(nep,nwork));
664 499 : xr = nep->work[issplit+1];
665 : #if !defined(PETSC_USE_COMPLEX)
666 499 : xi = nep->work[issplit+2];
667 : #endif
668 :
669 : /* compute residual norms */
670 499 : PetscCall(NEPGetEigenpair(nep,i,&kr,&ki,xr,xi));
671 : #if !defined(PETSC_USE_COMPLEX)
672 499 : PetscCheck(ki==0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Not implemented for complex eigenvalues with real scalars");
673 : #endif
674 499 : PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error));
675 499 : PetscCall(VecNorm(xr,NORM_2,&er));
676 :
677 : /* if two-sided, compute left residual norm and take the maximum */
678 499 : 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 18 : *error = PetscMax(*error,errorl);
682 : }
683 :
684 : /* compute error */
685 499 : switch (type) {
686 : case NEP_ERROR_ABSOLUTE:
687 : break;
688 470 : case NEP_ERROR_RELATIVE:
689 470 : *error /= PetscAbsScalar(kr)*er;
690 470 : break;
691 28 : case NEP_ERROR_BACKWARD:
692 28 : if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
693 6 : PetscCall(NEPComputeFunction(nep,kr,nep->function,nep->function));
694 6 : PetscCall(MatHasOperation(nep->function,MATOP_NORM,&flg));
695 6 : PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
696 6 : PetscCall(MatNorm(nep->function,NORM_INFINITY,&nrm));
697 6 : *error /= nrm*er;
698 6 : break;
699 : }
700 : /* initialization of matrix norms */
701 22 : if (!nep->nrma[0]) {
702 24 : for (j=0;j<nep->nt;j++) {
703 17 : PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&flg));
704 17 : PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
705 17 : PetscCall(MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]));
706 : }
707 : }
708 73 : for (j=0;j<nep->nt;j++) {
709 51 : PetscCall(FNEvaluateFunction(nep->f[j],kr,&s));
710 51 : z = z + nep->nrma[j]*PetscAbsScalar(s);
711 : }
712 22 : *error /= z*er;
713 22 : break;
714 0 : default:
715 0 : SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
716 : }
717 499 : 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 2245 : PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)
744 : {
745 2245 : PetscInt i;
746 2245 : PetscScalar alpha;
747 :
748 2245 : PetscFunctionBegin;
749 2245 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
750 2245 : NEPCheckProblem(nep,1);
751 2245 : switch (nep->fui) {
752 1299 : case NEP_USER_INTERFACE_CALLBACK:
753 1299 : PetscCheck(nep->computefunction,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
754 1299 : PetscCall(PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0));
755 1299 : PetscCallBack("NEP user Function function",(*nep->computefunction)(nep,lambda,A,B,nep->functionctx));
756 1299 : PetscCall(PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0));
757 : break;
758 946 : case NEP_USER_INTERFACE_SPLIT:
759 946 : PetscCall(MatZeroEntries(A));
760 946 : if (A != B) PetscCall(MatZeroEntries(B));
761 3778 : for (i=0;i<nep->nt;i++) {
762 2832 : PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
763 2832 : PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
764 2832 : if (A != B) PetscCall(MatAXPY(B,alpha,nep->P[i],nep->mstrp));
765 : }
766 : break;
767 : }
768 2245 : 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 923 : PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)
793 : {
794 923 : PetscInt i;
795 923 : PetscScalar alpha;
796 :
797 923 : PetscFunctionBegin;
798 923 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
799 923 : NEPCheckProblem(nep,1);
800 923 : switch (nep->fui) {
801 322 : case NEP_USER_INTERFACE_CALLBACK:
802 322 : PetscCheck(nep->computejacobian,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
803 322 : PetscCall(PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0));
804 322 : PetscCallBack("NEP user Jacobian function",(*nep->computejacobian)(nep,lambda,A,nep->jacobianctx));
805 322 : PetscCall(PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0));
806 : break;
807 601 : case NEP_USER_INTERFACE_SPLIT:
808 601 : PetscCall(MatZeroEntries(A));
809 2400 : for (i=0;i<nep->nt;i++) {
810 1799 : PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
811 1799 : PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
812 : }
813 : break;
814 : }
815 923 : PetscFunctionReturn(PETSC_SUCCESS);
816 : }
|