Actual source code: svddefault.c
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: /*@
57: SVDStoppingBasic - Default routine to determine whether the outer singular value
58: solver iteration must be stopped.
60: Collective
62: Input Parameters:
63: + svd - the singular value solver context
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: `SVDStoppingBasic()` will stop if all requested singular values are converged, or if
75: the maximum number of iterations has been reached.
77: This is the default stopping test.
78: Use `SVDSetStoppingTest()` to provide your own test instead of using this one.
80: Level: advanced
82: .seealso: [](ch:svd), `SVDSetStoppingTest()`, `SVDStoppingThreshold()`, `SVDConvergedReason`, `SVDGetConvergedReason()`
83: @*/
84: PetscErrorCode SVDStoppingBasic(SVD svd,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nsv,SVDConvergedReason *reason,void *ctx)
85: {
86: PetscFunctionBegin;
87: *reason = SVD_CONVERGED_ITERATING;
88: if (nconv >= nsv && svd->conv != SVD_CONV_MAXIT) {
89: PetscCall(PetscInfo(svd,"Singular value solver finished successfully: %" PetscInt_FMT " singular triplets converged at iteration %" PetscInt_FMT "\n",nconv,its));
90: *reason = SVD_CONVERGED_TOL;
91: } else if (its >= max_it) {
92: if (svd->conv == SVD_CONV_MAXIT) *reason = SVD_CONVERGED_MAXIT;
93: else {
94: *reason = SVD_DIVERGED_ITS;
95: PetscCall(PetscInfo(svd,"Singular value solver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
96: }
97: }
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: /*@
102: SVDStoppingThreshold - Routine to determine whether the outer singular value
103: solver iteration must be stopped, according to some threshold for the computed values.
105: Collective
107: Input Parameters:
108: + svd - the singular value solver context
109: . its - current number of iterations
110: . max_it - maximum number of iterations
111: . nconv - number of currently converged singular triplets (ignored here)
112: . nsv - number of requested singular triplets (ignored here)
113: - ctx - context containing additional data (`SVDStoppingCtx`)
115: Output Parameter:
116: . reason - result of the stopping test
118: Notes:
119: `SVDStoppingThreshold()` will stop when one of the computed singular values is not
120: above/below the threshold given at `SVDSetThreshold()`. If a number of wanted singular
121: values has been specified via `SVDSetDimensions()` then it is also taken into account,
122: and the solver will stop when one of the two conditions (threshold or number of
123: converged values) is met.
125: Use `SVDSetStoppingTest()` to provide your own test instead of using this one.
127: Level: advanced
129: .seealso: [](ch:svd), `SVDSetStoppingTest()`, `SVDStoppingBasic()`, `SVDSetThreshold()`, `SVDSetDimensions()`, `SVDConvergedReason`, `SVDGetConvergedReason()`
130: @*/
131: PetscErrorCode SVDStoppingThreshold(SVD svd,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nsv,SVDConvergedReason *reason,void *ctx)
132: {
133: PetscReal thres,firstsv,firstnc,errest;
134: PetscBool rel;
135: PetscInt napprox;
136: SVDWhich which;
138: PetscFunctionBegin;
139: *reason = SVD_CONVERGED_ITERATING;
140: firstsv = ((SVDStoppingCtx)ctx)->firstsv;
141: firstnc = ((SVDStoppingCtx)ctx)->firstnc;
142: errest = ((SVDStoppingCtx)ctx)->errest;
143: thres = ((SVDStoppingCtx)ctx)->thres;
144: rel = ((SVDStoppingCtx)ctx)->threlative;
145: napprox = ((SVDStoppingCtx)ctx)->napprox;
146: which = ((SVDStoppingCtx)ctx)->which;
148: if (!nconv) { /* no converged singular values yet */
149: if (its >= max_it) {
150: if (svd->conv == SVD_CONV_MAXIT) *reason = SVD_CONVERGED_MAXIT;
151: else {
152: PetscCall(PetscInfo(svd,"Singular value solver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
153: *reason = SVD_DIVERGED_ITS;
154: }
155: }
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
158: if (nconv==napprox) { /* all approximations converged, very unusual */
159: PetscCall(PetscInfo(svd,"Singular value solver finished successfully: all available singular value approximations have converged\n"));
160: *reason = SVD_CONVERGED_TOL;
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
163: if ((which==SVD_LARGEST && ((rel && firstnc+errest<thres*firstsv) || (!rel && firstnc+errest<thres))) || (which==SVD_SMALLEST && firstnc-errest>thres)) {
164: if (its==((SVDStoppingCtx)ctx)->its+1) {
165: if (which==SVD_SMALLEST) PetscCall(PetscInfo(svd,"Singular value solver finished successfully: the approximation %g (minus error estimate %g) is above the threshold %g\n",(double)firstnc,(double)errest,(double)thres));
166: else if (!rel) PetscCall(PetscInfo(svd,"Singular value solver finished successfully: the approximation %g (plus error estimate %g) is below the threshold %g\n",(double)firstnc,(double)errest,(double)thres));
167: else PetscCall(PetscInfo(svd,"Singular value solver finished successfully: the ratio %g/%g is below the threshold %g\n",(double)(firstnc+errest),(double)firstsv,(double)thres));
168: *reason = SVD_CONVERGED_TOL;
169: } else ((SVDStoppingCtx)ctx)->its = its; /* wait until next iteration */
170: } else if (nsv && nconv >= nsv) {
171: PetscCall(PetscInfo(svd,"Singular value solver finished successfully: %" PetscInt_FMT " singular triplets converged at iteration %" PetscInt_FMT "\n",nconv,its));
172: *reason = SVD_CONVERGED_TOL;
173: } else if (its >= max_it) {
174: if (svd->conv == SVD_CONV_MAXIT) *reason = SVD_CONVERGED_MAXIT;
175: else {
176: PetscCall(PetscInfo(svd,"Singular value solver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
177: *reason = SVD_DIVERGED_ITS;
178: }
179: }
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: /*@
184: SVDSetWorkVecs - Sets a number of work vectors into an `SVD` object.
186: Collective
188: Input Parameters:
189: + svd - the singular value solver context
190: . nleft - number of work vectors of dimension equal to left singular vector
191: - nright - number of work vectors of dimension equal to right singular vector
193: Developer Note:
194: This is `SLEPC_EXTERN` because it may be required by user plugin `SVD`
195: implementations.
197: Level: developer
199: .seealso: [](ch:svd), `SVDSetUp()`
200: @*/
201: PetscErrorCode SVDSetWorkVecs(SVD svd,PetscInt nleft,PetscInt nright)
202: {
203: Vec t;
205: PetscFunctionBegin;
209: PetscCheck(nleft>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nleft must be >= 0: nleft = %" PetscInt_FMT,nleft);
210: PetscCheck(nright>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nright must be >= 0: nright = %" PetscInt_FMT,nright);
211: PetscCheck(nleft>0 || nright>0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nleft and nright cannot be both zero");
212: if (svd->nworkl < nleft) {
213: PetscCall(VecDestroyVecs(svd->nworkl,&svd->workl));
214: svd->nworkl = nleft;
215: if (svd->isgeneralized) PetscCall(SVDCreateLeftTemplate(svd,&t));
216: else PetscCall(MatCreateVecsEmpty(svd->OP,NULL,&t));
217: PetscCall(VecDuplicateVecs(t,nleft,&svd->workl));
218: PetscCall(VecDestroy(&t));
219: }
220: if (svd->nworkr < nright) {
221: PetscCall(VecDestroyVecs(svd->nworkr,&svd->workr));
222: svd->nworkr = nright;
223: PetscCall(MatCreateVecsEmpty(svd->OP,&t,NULL));
224: PetscCall(VecDuplicateVecs(t,nright,&svd->workr));
225: PetscCall(VecDestroy(&t));
226: }
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }