Actual source code: nepdefault.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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>     /*I "slepcnep.h" I*/

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

 19:    Collective on NEP

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

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

 29:    Level: developer
 30: @*/
 31: PetscErrorCode NEPSetWorkVecs(NEP nep,PetscInt nw)
 32: {
 34:   Vec            t;

 37:   if (nep->nwork < nw) {
 38:     VecDestroyVecs(nep->nwork,&nep->work);
 39:     nep->nwork = nw;
 40:     BVGetColumn(nep->V,0,&t);
 41:     VecDuplicateVecs(t,nw,&nep->work);
 42:     BVRestoreColumn(nep->V,0,&t);
 43:     PetscLogObjectParents(nep,nw,nep->work);
 44:   }
 45:   return(0);
 46: }

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

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

 89:   w = SlepcAbsEigenvalue(eigr,eigi);
 90:   *errest = res/w;
 91:   return(0);
 92: }

 94: /*
 95:   NEPConvergedAbsolute - Checks convergence absolutely.
 96: */
 97: PetscErrorCode NEPConvergedAbsolute(NEP nep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
 98: {
100:   *errest = res;
101:   return(0);
102: }

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

116:   if (nep->fui!=NEP_USER_INTERFACE_SPLIT) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Backward error only available in split form");
117:   /* initialization of matrix norms */
118:   if (!nep->nrma[0]) {
119:     for (j=0;j<nep->nt;j++) {
120:       MatHasOperation(nep->A[j],MATOP_NORM,&flg);
121:       if (!flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The convergence test related to the matrix norms requires a matrix norm operation");
122:       MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]);
123:     }
124:   }
125:   for (j=0;j<nep->nt;j++) {
126:     FNEvaluateFunction(nep->f[j],eigr,&s);
127:     w = w + nep->nrma[j]*PetscAbsScalar(s);
128:   }
129:   *errest = res/w;
130:   return(0);
131: }

133: /*@C
134:    NEPStoppingBasic - Default routine to determine whether the outer eigensolver
135:    iteration must be stopped.

137:    Collective on NEP

139:    Input Parameters:
140: +  nep    - nonlinear eigensolver context obtained from NEPCreate()
141: .  its    - current number of iterations
142: .  max_it - maximum number of iterations
143: .  nconv  - number of currently converged eigenpairs
144: .  nev    - number of requested eigenpairs
145: -  ctx    - context (not used here)

147:    Output Parameter:
148: .  reason - result of the stopping test

150:    Notes:
151:    A positive value of reason indicates that the iteration has finished successfully
152:    (converged), and a negative value indicates an error condition (diverged). If
153:    the iteration needs to be continued, reason must be set to NEP_CONVERGED_ITERATING
154:    (zero).

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

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

161:    Level: advanced

163: .seealso: NEPSetStoppingTest(), NEPConvergedReason, NEPGetConvergedReason()
164: @*/
165: PetscErrorCode NEPStoppingBasic(NEP nep,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,NEPConvergedReason *reason,void *ctx)
166: {

170:   *reason = NEP_CONVERGED_ITERATING;
171:   if (nconv >= nev) {
172:     PetscInfo2(nep,"Nonlinear eigensolver finished successfully: %D eigenpairs converged at iteration %D\n",nconv,its);
173:     *reason = NEP_CONVERGED_TOL;
174:   } else if (its >= max_it) {
175:     *reason = NEP_DIVERGED_ITS;
176:     PetscInfo1(nep,"Nonlinear eigensolver iteration reached maximum number of iterations (%D)\n",its);
177:   }
178:   return(0);
179: }

181: PetscErrorCode NEPComputeVectors_Schur(NEP nep)
182: {
184:   PetscInt       n,i;
185:   Mat            Z;
186:   Vec            v;

189:   DSGetDimensions(nep->ds,&n,NULL,NULL,NULL,NULL);
190:   DSVectors(nep->ds,DS_MAT_X,NULL,NULL);
191:   DSGetMat(nep->ds,DS_MAT_X,&Z);
192:   BVSetActiveColumns(nep->V,0,n);
193:   BVMultInPlace(nep->V,Z,0,n);
194:   MatDestroy(&Z);

196:   /* normalization */
197:   for (i=0;i<n;i++) {
198:     BVGetColumn(nep->V,i,&v);
199:     VecNormalize(v,NULL);
200:     BVRestoreColumn(nep->V,i,&v);
201:   }
202:   return(0);
203: }