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 SVD operations
 12: */

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

 16: /*
 17:   SVDConvergedAbsolute - Checks convergence absolutely.
 18: */
 19: PetscErrorCode SVDConvergedAbsolute(SVD svd,PetscReal sigma,PetscReal res,PetscReal *errest,void *ctx)
 20: {
 21:   PetscFunctionBegin;
 22:   *errest = res;
 23:   PetscFunctionReturn(PETSC_SUCCESS);
 24: }

 26: /*
 27:   SVDConvergedRelative - Checks convergence relative to the singular value.
 28: */
 29: PetscErrorCode SVDConvergedRelative(SVD svd,PetscReal sigma,PetscReal res,PetscReal *errest,void *ctx)
 30: {
 31:   PetscFunctionBegin;
 32:   *errest = (sigma!=0.0)? res/sigma: PETSC_MAX_REAL;
 33:   PetscFunctionReturn(PETSC_SUCCESS);
 34: }

 36: /*
 37:   SVDConvergedNorm - Checks convergence relative to the matrix norms.
 38: */
 39: PetscErrorCode SVDConvergedNorm(SVD svd,PetscReal sigma,PetscReal res,PetscReal *errest,void *ctx)
 40: {
 41:   PetscFunctionBegin;
 42:   *errest = res/PetscMax(svd->nrma,svd->nrmb);
 43:   PetscFunctionReturn(PETSC_SUCCESS);
 44: }

 46: /*
 47:   SVDConvergedMaxIt - Always returns Inf to force reaching the maximum number of iterations.
 48: */
 49: PetscErrorCode SVDConvergedMaxIt(SVD svd,PetscReal sigma,PetscReal res,PetscReal *errest,void *ctx)
 50: {
 51:   PetscFunctionBegin;
 52:   *errest = PETSC_MAX_REAL;
 53:   PetscFunctionReturn(PETSC_SUCCESS);
 54: }

 56: /*@C
 57:    SVDStoppingBasic - Default routine to determine whether the outer singular value
 58:    solver iteration must be stopped.

 60:    Collective

 62:    Input Parameters:
 63: +  svd    - singular value solver context obtained from SVDCreate()
 64: .  its    - current number of iterations
 65: .  max_it - maximum number of iterations
 66: .  nconv  - number of currently converged singular triplets
 67: .  nsv    - number of requested singular triplets
 68: -  ctx    - context (not used here)

 70:    Output Parameter:
 71: .  reason - result of the stopping test

 73:    Notes:
 74:    A positive value of reason indicates that the iteration has finished successfully
 75:    (converged), and a negative value indicates an error condition (diverged). If
 76:    the iteration needs to be continued, reason must be set to SVD_CONVERGED_ITERATING
 77:    (zero).

 79:    SVDStoppingBasic() will stop if all requested singular values are converged, or if
 80:    the maximum number of iterations has been reached.

 82:    Use SVDSetStoppingTest() to provide your own test instead of using this one.

 84:    Level: advanced

 86: .seealso: SVDSetStoppingTest(), SVDStoppingThreshold(), SVDConvergedReason, SVDGetConvergedReason()
 87: @*/
 88: PetscErrorCode SVDStoppingBasic(SVD svd,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nsv,SVDConvergedReason *reason,void *ctx)
 89: {
 90:   PetscFunctionBegin;
 91:   *reason = SVD_CONVERGED_ITERATING;
 92:   if (nconv >= nsv) {
 93:     PetscCall(PetscInfo(svd,"Singular value solver finished successfully: %" PetscInt_FMT " singular triplets converged at iteration %" PetscInt_FMT "\n",nconv,its));
 94:     *reason = SVD_CONVERGED_TOL;
 95:   } else if (its >= max_it) {
 96:     if (svd->conv == SVD_CONV_MAXIT) *reason = SVD_CONVERGED_MAXIT;
 97:     else {
 98:       *reason = SVD_DIVERGED_ITS;
 99:       PetscCall(PetscInfo(svd,"Singular value solver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
100:     }
101:   }
102:   PetscFunctionReturn(PETSC_SUCCESS);
103: }

105: /*@C
106:    SVDStoppingThreshold - Routine to determine whether the outer singular value
107:    solver iteration must be stopped, according to some threshold for the computed values.

109:    Collective

111:    Input Parameters:
112: +  svd    - singular value solver context obtained from SVDCreate()
113: .  its    - current number of iterations
114: .  max_it - maximum number of iterations
115: .  nconv  - number of currently converged singular triplets (ignored here)
116: .  nsv    - number of requested singular triplets (ignored here)
117: -  ctx    - context containing additional data (SVDStoppingCtx)

119:    Output Parameter:
120: .  reason - result of the stopping test

122:    Notes:
123:    A positive value of reason indicates that the iteration has finished successfully
124:    (converged), and a negative value indicates an error condition (diverged). If
125:    the iteration needs to be continued, reason must be set to SVD_CONVERGED_ITERATING
126:    (zero).

128:    SVDStoppingThreshold() will stop when one of the computed singular values is not
129:    above/below the threshold given at SVDSetThreshold(). If a number of wanted singular
130:    values has been specified via SVDSetDimensions() then it is also taken into account,
131:    and the solver will stop when one of the two conditions (threshold or number of
132:    converged values) is met.

134:    Use SVDSetStoppingTest() to provide your own test instead of using this one.

136:    Level: advanced

138: .seealso: SVDSetStoppingTest(), SVDStoppingBasic(), SVDSetThreshold(), SVDSetDimensions(), SVDConvergedReason, SVDGetConvergedReason()
139: @*/
140: PetscErrorCode SVDStoppingThreshold(SVD svd,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nsv,SVDConvergedReason *reason,void *ctx)
141: {
142:   PetscReal thres,firstsv,lastsv;
143:   PetscBool rel;
144:   SVDWhich  which;

146:   PetscFunctionBegin;
147:   *reason = SVD_CONVERGED_ITERATING;
148:   firstsv = ((SVDStoppingCtx)ctx)->firstsv;
149:   lastsv  = ((SVDStoppingCtx)ctx)->lastsv;
150:   thres   = ((SVDStoppingCtx)ctx)->thres;
151:   rel     = ((SVDStoppingCtx)ctx)->threlative;
152:   which   = ((SVDStoppingCtx)ctx)->which;
153:   if (nconv && ((which==SVD_LARGEST && ((rel && lastsv<thres*firstsv) || (!rel && lastsv<thres))) || (which==SVD_SMALLEST && lastsv>thres))) {
154:     if (which==SVD_SMALLEST) PetscCall(PetscInfo(svd,"Singular value solver finished successfully: singular value %g is above the threshold %g\n",(double)lastsv,(double)thres));
155:     else if (!rel) PetscCall(PetscInfo(svd,"Singular value solver finished successfully: singular value %g is below the threshold %g\n",(double)lastsv,(double)thres));
156:     else PetscCall(PetscInfo(svd,"Singular value solver finished successfully: the ratio %g/%g is below the threshold %g\n",(double)lastsv,(double)firstsv,(double)thres));
157:     *reason = SVD_CONVERGED_TOL;
158:   } else if (nsv && nconv >= nsv) {
159:     PetscCall(PetscInfo(svd,"Singular value solver finished successfully: %" PetscInt_FMT " singular triplets converged at iteration %" PetscInt_FMT "\n",nconv,its));
160:     *reason = SVD_CONVERGED_TOL;
161:   } else if (its >= max_it) {
162:     if (svd->conv == SVD_CONV_MAXIT) *reason = SVD_CONVERGED_MAXIT;
163:     else {
164:       *reason = SVD_DIVERGED_ITS;
165:       PetscCall(PetscInfo(svd,"Singular value solver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
166:     }
167:   }
168:   PetscFunctionReturn(PETSC_SUCCESS);
169: }

171: /*@
172:    SVDSetWorkVecs - Sets a number of work vectors into an SVD object.

174:    Collective

176:    Input Parameters:
177: +  svd    - singular value solver context
178: .  nleft  - number of work vectors of dimension equal to left singular vector
179: -  nright - number of work vectors of dimension equal to right singular vector

181:    Developer Notes:
182:    This is SLEPC_EXTERN because it may be required by user plugin SVD
183:    implementations.

185:    Level: developer

187: .seealso: SVDSetUp()
188: @*/
189: PetscErrorCode SVDSetWorkVecs(SVD svd,PetscInt nleft,PetscInt nright)
190: {
191:   Vec            t;

193:   PetscFunctionBegin;
197:   PetscCheck(nleft>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nleft must be >= 0: nleft = %" PetscInt_FMT,nleft);
198:   PetscCheck(nright>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nright must be >= 0: nright = %" PetscInt_FMT,nright);
199:   PetscCheck(nleft>0 || nright>0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nleft and nright cannot be both zero");
200:   if (svd->nworkl < nleft) {
201:     PetscCall(VecDestroyVecs(svd->nworkl,&svd->workl));
202:     svd->nworkl = nleft;
203:     if (svd->isgeneralized) PetscCall(SVDCreateLeftTemplate(svd,&t));
204:     else PetscCall(MatCreateVecsEmpty(svd->OP,NULL,&t));
205:     PetscCall(VecDuplicateVecs(t,nleft,&svd->workl));
206:     PetscCall(VecDestroy(&t));
207:   }
208:   if (svd->nworkr < nright) {
209:     PetscCall(VecDestroyVecs(svd->nworkr,&svd->workr));
210:     svd->nworkr = nright;
211:     PetscCall(MatCreateVecsEmpty(svd->OP,&t,NULL));
212:     PetscCall(VecDuplicateVecs(t,nright,&svd->workr));
213:     PetscCall(VecDestroy(&t));
214:   }
215:   PetscFunctionReturn(PETSC_SUCCESS);
216: }