1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2018, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: NEP routines related to the solution process
12: */
14: #include <slepc/private/nepimpl.h> /*I "slepcnep.h" I*/
15: #include <petscdraw.h>
17: PetscErrorCode NEPComputeVectors(NEP nep) 18: {
22: NEPCheckSolved(nep,1);
23: if (nep->state==NEP_STATE_SOLVED && nep->ops->computevectors) {
24: (*nep->ops->computevectors)(nep);
25: }
26: nep->state = NEP_STATE_EIGENVECTORS;
27: return(0);
28: }
30: /*@
31: NEPSolve - Solves the nonlinear eigensystem.
33: Collective on NEP 35: Input Parameter:
36: . nep - eigensolver context obtained from NEPCreate()
38: Options Database Keys:
39: + -nep_view - print information about the solver used
40: . -nep_view_vectors binary - save the computed eigenvectors to the default binary viewer
41: . -nep_view_values - print computed eigenvalues
42: . -nep_converged_reason - print reason for convergence, and number of iterations
43: . -nep_error_absolute - print absolute errors of each eigenpair
44: - -nep_error_relative - print relative errors of each eigenpair
46: Level: beginner
48: .seealso: NEPCreate(), NEPSetUp(), NEPDestroy(), NEPSetTolerances()
49: @*/
50: PetscErrorCode NEPSolve(NEP nep) 51: {
53: PetscInt i;
57: if (nep->state>=NEP_STATE_SOLVED) return(0);
58: PetscLogEventBegin(NEP_Solve,nep,0,0,0);
60: /* call setup */
61: NEPSetUp(nep);
62: nep->nconv = 0;
63: nep->its = 0;
64: for (i=0;i<nep->ncv;i++) {
65: nep->eigr[i] = 0.0;
66: nep->eigi[i] = 0.0;
67: nep->errest[i] = 0.0;
68: nep->perm[i] = i;
69: }
70: NEPViewFromOptions(nep,NULL,"-nep_view_pre");
71: RGViewFromOptions(nep->rg,NULL,"-rg_view");
73: (*nep->ops->solve)(nep);
74: nep->state = NEP_STATE_SOLVED;
76: if (!nep->reason) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
78: if (nep->refine==NEP_REFINE_SIMPLE && nep->rits>0 && nep->nconv>0) {
79: NEPComputeVectors(nep);
80: NEPNewtonRefinementSimple(nep,&nep->rits,nep->rtol,nep->nconv);
81: nep->state = NEP_STATE_EIGENVECTORS;
82: }
84: /* sort eigenvalues according to nep->which parameter */
85: SlepcSortEigenvalues(nep->sc,nep->nconv,nep->eigr,nep->eigi,nep->perm);
86: PetscLogEventEnd(NEP_Solve,nep,0,0,0);
88: /* various viewers */
89: NEPViewFromOptions(nep,NULL,"-nep_view");
90: NEPReasonViewFromOptions(nep);
91: NEPErrorViewFromOptions(nep);
92: NEPValuesViewFromOptions(nep);
93: NEPVectorsViewFromOptions(nep);
95: /* Remove the initial subspace */
96: nep->nini = 0;
97: return(0);
98: }
100: /*@
101: NEPProjectOperator - Computes the projection of the nonlinear operator.
103: Collective on NEP105: Input Parameters:
106: + nep - the nonlinear eigensolver context
107: . j0 - initial index
108: - j1 - final index
110: Notes:
111: This is available for split operator only.
113: The nonlinear operator T(lambda) is projected onto span(V), where V is
114: an orthonormal basis built internally by the solver. The projected
115: operator is equal to sum_i V'*A_i*V*f_i(lambda), so this function
116: computes all matrices Ei = V'*A_i*V, and stores them in the extra
117: matrices inside DS. Only rows/columns in the range [j0,j1-1] are computed,
118: the previous ones are assumed to be available already.
120: Level: developer
122: .seealso: NEPSetSplitOperator()
123: @*/
124: PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)125: {
127: PetscInt k;
128: Mat G;
134: NEPCheckProblem(nep,1);
135: NEPCheckSplit(nep,1);
136: BVSetActiveColumns(nep->V,j0,j1);
137: for (k=0;k<nep->nt;k++) {
138: DSGetMat(nep->ds,DSMatExtra[k],&G);
139: BVMatProject(nep->V,nep->A[k],nep->V,G);
140: DSRestoreMat(nep->ds,DSMatExtra[k],&G);
141: }
142: return(0);
143: }
145: /*@
146: NEPApplyFunction - Applies the nonlinear function T(lambda) to a given vector.
148: Collective on NEP150: Input Parameters:
151: + nep - the nonlinear eigensolver context
152: . lambda - scalar argument
153: . x - vector to be multiplied against
154: - v - workspace vector (used only in the case of split form)
156: Output Parameters:
157: + y - result vector
158: . A - Function matrix
159: - B - optional preconditioning matrix
161: Note:
162: If the nonlinear operator is represented in split form, the result
163: y = T(lambda)*x is computed without building T(lambda) explicitly. In
164: that case, parameters A and B are not used. Otherwise, the matrix
165: T(lambda) is built and the effect is the same as a call to
166: NEPComputeFunction() followed by a MatMult().
168: Level: developer
170: .seealso: NEPSetSplitOperator(), NEPComputeFunction()
171: @*/
172: PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)173: {
175: PetscInt i;
176: PetscScalar alpha;
187: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
188: VecSet(y,0.0);
189: for (i=0;i<nep->nt;i++) {
190: FNEvaluateFunction(nep->f[i],lambda,&alpha);
191: MatMult(nep->A[i],x,v);
192: VecAXPY(y,alpha,v);
193: }
194: } else {
195: NEPComputeFunction(nep,lambda,A,B);
196: MatMult(A,x,y);
197: }
198: return(0);
199: }
201: /*@
202: NEPApplyJacobian - Applies the nonlinear Jacobian T'(lambda) to a given vector.
204: Collective on NEP206: Input Parameters:
207: + nep - the nonlinear eigensolver context
208: . lambda - scalar argument
209: . x - vector to be multiplied against
210: - v - workspace vector (used only in the case of split form)
212: Output Parameters:
213: + y - result vector
214: - A - Jacobian matrix
216: Note:
217: If the nonlinear operator is represented in split form, the result
218: y = T'(lambda)*x is computed without building T'(lambda) explicitly. In
219: that case, parameter A is not used. Otherwise, the matrix
220: T'(lambda) is built and the effect is the same as a call to
221: NEPComputeJacobian() followed by a MatMult().
223: Level: developer
225: .seealso: NEPSetSplitOperator(), NEPComputeJacobian()
226: @*/
227: PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)228: {
230: PetscInt i;
231: PetscScalar alpha;
241: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
242: VecSet(y,0.0);
243: for (i=0;i<nep->nt;i++) {
244: FNEvaluateDerivative(nep->f[i],lambda,&alpha);
245: MatMult(nep->A[i],x,v);
246: VecAXPY(y,alpha,v);
247: }
248: } else {
249: NEPComputeJacobian(nep,lambda,A);
250: MatMult(A,x,y);
251: }
252: return(0);
253: }
255: /*@
256: NEPGetIterationNumber - Gets the current iteration number. If the
257: call to NEPSolve() is complete, then it returns the number of iterations
258: carried out by the solution method.
260: Not Collective
262: Input Parameter:
263: . nep - the nonlinear eigensolver context
265: Output Parameter:
266: . its - number of iterations
268: Note:
269: During the i-th iteration this call returns i-1. If NEPSolve() is
270: complete, then parameter "its" contains either the iteration number at
271: which convergence was successfully reached, or failure was detected.
272: Call NEPGetConvergedReason() to determine if the solver converged or
273: failed and why.
275: Level: intermediate
277: .seealso: NEPGetConvergedReason(), NEPSetTolerances()
278: @*/
279: PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)280: {
284: *its = nep->its;
285: return(0);
286: }
288: /*@
289: NEPGetConverged - Gets the number of converged eigenpairs.
291: Not Collective
293: Input Parameter:
294: . nep - the nonlinear eigensolver context
296: Output Parameter:
297: . nconv - number of converged eigenpairs
299: Note:
300: This function should be called after NEPSolve() has finished.
302: Level: beginner
304: .seealso: NEPSetDimensions(), NEPSolve()
305: @*/
306: PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)307: {
311: NEPCheckSolved(nep,1);
312: *nconv = nep->nconv;
313: return(0);
314: }
316: /*@
317: NEPGetConvergedReason - Gets the reason why the NEPSolve() iteration was
318: stopped.
320: Not Collective
322: Input Parameter:
323: . nep - the nonlinear eigensolver context
325: Output Parameter:
326: . reason - negative value indicates diverged, positive value converged
328: Notes:
330: Possible values for reason are
331: + NEP_CONVERGED_TOL - converged up to tolerance
332: . NEP_CONVERGED_USER - converged due to a user-defined condition
333: . NEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
334: . NEP_DIVERGED_BREAKDOWN - generic breakdown in method
335: . NEP_DIVERGED_LINEAR_SOLVE - inner linear solve failed
336: - NEP_DIVERGED_SUBSPACE_EXHAUSTED - run out of space for the basis in an
337: unrestarted solver
339: Can only be called after the call to NEPSolve() is complete.
341: Level: intermediate
343: .seealso: NEPSetTolerances(), NEPSolve(), NEPConvergedReason344: @*/
345: PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)346: {
350: NEPCheckSolved(nep,1);
351: *reason = nep->reason;
352: return(0);
353: }
355: /*@C
356: NEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
357: NEPSolve(). The solution consists in both the eigenvalue and the eigenvector.
359: Logically Collective on NEP361: Input Parameters:
362: + nep - nonlinear eigensolver context
363: - i - index of the solution
365: Output Parameters:
366: + eigr - real part of eigenvalue
367: . eigi - imaginary part of eigenvalue
368: . Vr - real part of eigenvector
369: - Vi - imaginary part of eigenvector
371: Notes:
372: It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
373: required. Otherwise, the caller must provide valid Vec objects, i.e.,
374: they must be created by the calling program with e.g. MatCreateVecs().
376: If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
377: configured with complex scalars the eigenvalue is stored
378: directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
379: set to zero). In both cases, the user can pass NULL in eigi and Vi.
381: The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
382: Eigenpairs are indexed according to the ordering criterion established
383: with NEPSetWhichEigenpairs().
385: Level: beginner
387: .seealso: NEPSolve(), NEPGetConverged(), NEPSetWhichEigenpairs()
388: @*/
389: PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)390: {
391: PetscInt k;
399: NEPCheckSolved(nep,1);
400: if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
402: NEPComputeVectors(nep);
403: k = nep->perm[i];
405: /* eigenvalue */
406: #if defined(PETSC_USE_COMPLEX)
407: if (eigr) *eigr = nep->eigr[k];
408: if (eigi) *eigi = 0;
409: #else
410: if (eigr) *eigr = nep->eigr[k];
411: if (eigi) *eigi = nep->eigi[k];
412: #endif
414: /* eigenvector */
415: #if defined(PETSC_USE_COMPLEX)
416: if (Vr) { BVCopyVec(nep->V,k,Vr); }
417: if (Vi) { VecSet(Vi,0.0); }
418: #else
419: if (nep->eigi[k]>0) { /* first value of conjugate pair */
420: if (Vr) { BVCopyVec(nep->V,k,Vr); }
421: if (Vi) { BVCopyVec(nep->V,k+1,Vi); }
422: } else if (nep->eigi[k]<0) { /* second value of conjugate pair */
423: if (Vr) { BVCopyVec(nep->V,k-1,Vr); }
424: if (Vi) {
425: BVCopyVec(nep->V,k,Vi);
426: VecScale(Vi,-1.0);
427: }
428: } else { /* real eigenvalue */
429: if (Vr) { BVCopyVec(nep->V,k,Vr); }
430: if (Vi) { VecSet(Vi,0.0); }
431: }
432: #endif
433: return(0);
434: }
436: /*@
437: NEPGetErrorEstimate - Returns the error estimate associated to the i-th
438: computed eigenpair.
440: Not Collective
442: Input Parameter:
443: + nep - nonlinear eigensolver context
444: - i - index of eigenpair
446: Output Parameter:
447: . errest - the error estimate
449: Notes:
450: This is the error estimate used internally by the eigensolver. The actual
451: error bound can be computed with NEPComputeRelativeError().
453: Level: advanced
455: .seealso: NEPComputeRelativeError()
456: @*/
457: PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)458: {
462: NEPCheckSolved(nep,1);
463: if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
464: *errest = nep->errest[nep->perm[i]];
465: return(0);
466: }
468: /*
469: NEPComputeResidualNorm_Private - Computes the norm of the residual vector
470: associated with an eigenpair.
472: Input Parameters:
473: lambda - eigenvalue
474: x - eigenvector
475: w - array of work vectors (two vectors in split form, one vector otherwise)
476: */
477: PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)478: {
480: Vec y,z=NULL;
483: y = w[0];
484: if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
485: NEPApplyFunction(nep,lambda,x,z,y,nep->function,nep->function_pre);
486: VecNorm(y,NORM_2,norm);
487: return(0);
488: }
490: /*@
491: NEPComputeError - Computes the error (based on the residual norm) associated
492: with the i-th computed eigenpair.
494: Collective on NEP496: Input Parameter:
497: + nep - the nonlinear eigensolver context
498: . i - the solution index
499: - type - the type of error to compute
501: Output Parameter:
502: . error - the error
504: Notes:
505: The error can be computed in various ways, all of them based on the residual
506: norm computed as ||T(lambda)x||_2 where lambda is the eigenvalue and x is the
507: eigenvector.
509: Level: beginner
511: .seealso: NEPErrorType, NEPSolve(), NEPGetErrorEstimate()
512: @*/
513: PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)514: {
516: Vec xr,xi=NULL;
517: PetscInt j,nwork,issplit=0;
518: PetscScalar kr,ki,s;
519: PetscReal er,z=0.0;
520: PetscBool flg;
527: NEPCheckSolved(nep,1);
529: /* allocate work vectors */
530: #if defined(PETSC_USE_COMPLEX)
531: nwork = 2;
532: #else
533: nwork = 3;
534: #endif
535: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
536: issplit = 1;
537: nwork++; /* need an extra work vector for NEPComputeResidualNorm_Private */
538: }
539: NEPSetWorkVecs(nep,nwork);
540: xr = nep->work[issplit+1];
541: #if !defined(PETSC_USE_COMPLEX)
542: xi = nep->work[issplit+2];
543: #endif
545: /* compute residual norms */
546: NEPGetEigenpair(nep,i,&kr,&ki,xr,xi);
547: #if !defined(PETSC_USE_COMPLEX)
548: if (ki) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Not implemented for complex eigenvalues with real scalars");
549: #endif
550: NEPComputeResidualNorm_Private(nep,kr,xr,nep->work,error);
551: VecNorm(xr,NORM_2,&er);
553: /* compute error */
554: switch (type) {
555: case NEP_ERROR_ABSOLUTE:
556: break;
557: case NEP_ERROR_RELATIVE:
558: *error /= PetscAbsScalar(kr)*er;
559: break;
560: case NEP_ERROR_BACKWARD:
561: if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
562: *error = 0.0;
563: PetscInfo(nep,"Backward error only available in split form\n");
564: break;
565: }
566: /* initialization of matrix norms */
567: if (!nep->nrma[0]) {
568: for (j=0;j<nep->nt;j++) {
569: MatHasOperation(nep->A[j],MATOP_NORM,&flg);
570: if (!flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
571: MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]);
572: }
573: }
574: for (j=0;j<nep->nt;j++) {
575: FNEvaluateFunction(nep->f[j],kr,&s);
576: z = z + nep->nrma[j]*PetscAbsScalar(s);
577: }
578: *error /= z;
579: break;
580: default:581: SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
582: }
583: return(0);
584: }
586: /*@
587: NEPComputeFunction - Computes the function matrix T(lambda) that has been
588: set with NEPSetFunction().
590: Collective on NEP and Mat
592: Input Parameters:
593: + nep - the NEP context
594: - lambda - the scalar argument
596: Output Parameters:
597: + A - Function matrix
598: - B - optional preconditioning matrix
600: Notes:
601: NEPComputeFunction() is typically used within nonlinear eigensolvers
602: implementations, so most users would not generally call this routine
603: themselves.
605: Level: developer
607: .seealso: NEPSetFunction(), NEPGetFunction()
608: @*/
609: PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)610: {
612: PetscInt i;
613: PetscScalar alpha;
617: NEPCheckProblem(nep,1);
618: switch (nep->fui) {
619: case NEP_USER_INTERFACE_CALLBACK:
620: if (!nep->computefunction) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
621: PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0);
622: PetscStackPush("NEP user Function function");
623: (*nep->computefunction)(nep,lambda,A,B,nep->functionctx);
624: PetscStackPop;
625: PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0);
626: break;
627: case NEP_USER_INTERFACE_SPLIT:
628: MatZeroEntries(A);
629: for (i=0;i<nep->nt;i++) {
630: FNEvaluateFunction(nep->f[i],lambda,&alpha);
631: MatAXPY(A,alpha,nep->A[i],nep->mstr);
632: }
633: if (A != B) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Not implemented");
634: break;
635: case NEP_USER_INTERFACE_DERIVATIVES:
636: PetscLogEventBegin(NEP_DerivativesEval,nep,A,B,0);
637: PetscStackPush("NEP user Derivatives function");
638: (*nep->computederivatives)(nep,lambda,0,A,nep->derivativesctx);
639: PetscStackPop;
640: PetscLogEventEnd(NEP_DerivativesEval,nep,A,B,0);
641: break;
642: }
643: return(0);
644: }
646: /*@
647: NEPComputeJacobian - Computes the Jacobian matrix T'(lambda) that has been
648: set with NEPSetJacobian().
650: Collective on NEP and Mat
652: Input Parameters:
653: + nep - the NEP context
654: - lambda - the scalar argument
656: Output Parameters:
657: . A - Jacobian matrix
659: Notes:
660: Most users should not need to explicitly call this routine, as it
661: is used internally within the nonlinear eigensolvers.
663: Level: developer
665: .seealso: NEPSetJacobian(), NEPGetJacobian()
666: @*/
667: PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)668: {
670: PetscInt i;
671: PetscScalar alpha;
675: NEPCheckProblem(nep,1);
676: switch (nep->fui) {
677: case NEP_USER_INTERFACE_CALLBACK:
678: if (!nep->computejacobian) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
679: PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0);
680: PetscStackPush("NEP user Jacobian function");
681: (*nep->computejacobian)(nep,lambda,A,nep->jacobianctx);
682: PetscStackPop;
683: PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0);
684: break;
685: case NEP_USER_INTERFACE_SPLIT:
686: MatZeroEntries(A);
687: for (i=0;i<nep->nt;i++) {
688: FNEvaluateDerivative(nep->f[i],lambda,&alpha);
689: MatAXPY(A,alpha,nep->A[i],nep->mstr);
690: }
691: break;
692: case NEP_USER_INTERFACE_DERIVATIVES:
693: PetscLogEventBegin(NEP_DerivativesEval,nep,A,0,0);
694: PetscStackPush("NEP user Derivatives function");
695: (*nep->computederivatives)(nep,lambda,1,A,nep->derivativesctx);
696: PetscStackPop;
697: PetscLogEventEnd(NEP_DerivativesEval,nep,A,0,0);
698: break;
699: }
700: return(0);
701: }