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