1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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:    Simple default routines for common NEP operations
 12: */

 14: #include <slepc/private/nepimpl.h>

 16: /*@
 17:    NEPSetWorkVecs - Sets a number of work vectors into a NEP object

 19:    Collective

 21:    Input Parameters:
 22: +  nep - nonlinear eigensolver context
 23: -  nw  - number of work vectors to allocate

 25:    Developer Notes:
 26:    This is SLEPC_EXTERN because it may be required by user plugin NEP
 27:    implementations.

 29:    Level: developer

 31: .seealso: NEPSetUp()
 32: @*/
 33: PetscErrorCode NEPSetWorkVecs(NEP nep,PetscInt nw)
 34: {
 35:   Vec            t;

 37:   PetscFunctionBegin;
 40:   PetscCheck(nw>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %" PetscInt_FMT,nw);
 41:   if (nep->nwork < nw) {
 42:     PetscCall(VecDestroyVecs(nep->nwork,&nep->work));
 43:     nep->nwork = nw;
 44:     PetscCall(BVGetColumn(nep->V,0,&t));
 45:     PetscCall(VecDuplicateVecs(t,nw,&nep->work));
 46:     PetscCall(BVRestoreColumn(nep->V,0,&t));
 47:   }
 48:   PetscFunctionReturn(PETSC_SUCCESS);
 49: }

 51: /*
 52:   NEPGetDefaultShift - Return the value of sigma to start the nonlinear iteration.
 53:  */
 54: PetscErrorCode NEPGetDefaultShift(NEP nep,PetscScalar *sigma)
 55: {
 56:   PetscFunctionBegin;
 57:   PetscAssertPointer(sigma,2);
 58:   switch (nep->which) {
 59:     case NEP_LARGEST_MAGNITUDE:
 60:     case NEP_LARGEST_IMAGINARY:
 61:     case NEP_ALL:
 62:     case NEP_WHICH_USER:
 63:       *sigma = 1.0;   /* arbitrary value */
 64:       break;
 65:     case NEP_SMALLEST_MAGNITUDE:
 66:     case NEP_SMALLEST_IMAGINARY:
 67:       *sigma = 0.0;
 68:       break;
 69:     case NEP_LARGEST_REAL:
 70:       *sigma = PETSC_MAX_REAL;
 71:       break;
 72:     case NEP_SMALLEST_REAL:
 73:       *sigma = PETSC_MIN_REAL;
 74:       break;
 75:     case NEP_TARGET_MAGNITUDE:
 76:     case NEP_TARGET_REAL:
 77:     case NEP_TARGET_IMAGINARY:
 78:       *sigma = nep->target;
 79:       break;
 80:   }
 81:   PetscFunctionReturn(PETSC_SUCCESS);
 82: }

 84: /*
 85:   NEPConvergedRelative - Checks convergence relative to the eigenvalue.
 86: */
 87: PetscErrorCode NEPConvergedRelative(NEP nep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
 88: {
 89:   PetscReal w;

 91:   PetscFunctionBegin;
 92:   w = SlepcAbsEigenvalue(eigr,eigi);
 93:   *errest = (w!=0.0)? res/w: PETSC_MAX_REAL;
 94:   PetscFunctionReturn(PETSC_SUCCESS);
 95: }

 97: /*
 98:   NEPConvergedAbsolute - Checks convergence absolutely.
 99: */
100: PetscErrorCode NEPConvergedAbsolute(NEP nep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
101: {
102:   PetscFunctionBegin;
103:   *errest = res;
104:   PetscFunctionReturn(PETSC_SUCCESS);
105: }

107: /*
108:   NEPConvergedNorm - Checks convergence relative to the matrix norms.
109: */
110: PetscErrorCode NEPConvergedNorm(NEP nep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
111: {
112:   PetscScalar    s;
113:   PetscReal      w=0.0;
114:   PetscInt       j;
115:   PetscBool      flg;

117:   PetscFunctionBegin;
118:   if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
119:     PetscCall(NEPComputeFunction(nep,eigr,nep->function,nep->function));
120:     PetscCall(MatHasOperation(nep->function,MATOP_NORM,&flg));
121:     PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
122:     PetscCall(MatNorm(nep->function,NORM_INFINITY,&w));
123:   } else {
124:     /* initialization of matrix norms */
125:     if (!nep->nrma[0]) {
126:       for (j=0;j<nep->nt;j++) {
127:         PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&flg));
128:         PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The convergence test related to the matrix norms requires a matrix norm operation");
129:         PetscCall(MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]));
130:       }
131:     }
132:     for (j=0;j<nep->nt;j++) {
133:       PetscCall(FNEvaluateFunction(nep->f[j],eigr,&s));
134:       w = w + nep->nrma[j]*PetscAbsScalar(s);
135:     }
136:   }
137:   *errest = res/w;
138:   PetscFunctionReturn(PETSC_SUCCESS);
139: }

141: /*@C
142:    NEPStoppingBasic - Default routine to determine whether the outer eigensolver
143:    iteration must be stopped.

145:    Collective

147:    Input Parameters:
148: +  nep    - nonlinear eigensolver context obtained from NEPCreate()
149: .  its    - current number of iterations
150: .  max_it - maximum number of iterations
151: .  nconv  - number of currently converged eigenpairs
152: .  nev    - number of requested eigenpairs
153: -  ctx    - context (not used here)

155:    Output Parameter:
156: .  reason - result of the stopping test

158:    Notes:
159:    A positive value of reason indicates that the iteration has finished successfully
160:    (converged), and a negative value indicates an error condition (diverged). If
161:    the iteration needs to be continued, reason must be set to NEP_CONVERGED_ITERATING
162:    (zero).

164:    NEPStoppingBasic() will stop if all requested eigenvalues are converged, or if
165:    the maximum number of iterations has been reached.

167:    Use NEPSetStoppingTest() to provide your own test instead of using this one.

169:    Level: advanced

171: .seealso: NEPSetStoppingTest(), NEPConvergedReason, NEPGetConvergedReason()
172: @*/
173: PetscErrorCode NEPStoppingBasic(NEP nep,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,NEPConvergedReason *reason,void *ctx)
174: {
175:   PetscFunctionBegin;
176:   *reason = NEP_CONVERGED_ITERATING;
177:   if (nconv >= nev) {
178:     PetscCall(PetscInfo(nep,"Nonlinear eigensolver finished successfully: %" PetscInt_FMT " eigenpairs converged at iteration %" PetscInt_FMT "\n",nconv,its));
179:     *reason = NEP_CONVERGED_TOL;
180:   } else if (its >= max_it) {
181:     *reason = NEP_DIVERGED_ITS;
182:     PetscCall(PetscInfo(nep,"Nonlinear eigensolver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
183:   }
184:   PetscFunctionReturn(PETSC_SUCCESS);
185: }

187: PetscErrorCode NEPComputeVectors_Schur(NEP nep)
188: {
189:   Mat            Z;

191:   PetscFunctionBegin;
192:   PetscCall(DSVectors(nep->ds,DS_MAT_X,NULL,NULL));
193:   PetscCall(DSGetMat(nep->ds,DS_MAT_X,&Z));
194:   PetscCall(BVMultInPlace(nep->V,Z,0,nep->nconv));
195:   PetscCall(DSRestoreMat(nep->ds,DS_MAT_X,&Z));
196:   PetscCall(BVNormalize(nep->V,nep->eigi));
197:   PetscFunctionReturn(PETSC_SUCCESS);
198: }