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: /*@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 - 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: /*@C
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,lastsv;
134: PetscBool rel;
135: SVDWhich which;
137: PetscFunctionBegin;
138: *reason = SVD_CONVERGED_ITERATING;
139: firstsv = ((SVDStoppingCtx)ctx)->firstsv;
140: lastsv = ((SVDStoppingCtx)ctx)->lastsv;
141: thres = ((SVDStoppingCtx)ctx)->thres;
142: rel = ((SVDStoppingCtx)ctx)->threlative;
143: which = ((SVDStoppingCtx)ctx)->which;
144: if (nconv && ((which==SVD_LARGEST && ((rel && lastsv<thres*firstsv) || (!rel && lastsv<thres))) || (which==SVD_SMALLEST && lastsv>thres))) {
145: 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));
146: else if (!rel) PetscCall(PetscInfo(svd,"Singular value solver finished successfully: singular value %g is below the threshold %g\n",(double)lastsv,(double)thres));
147: 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));
148: *reason = SVD_CONVERGED_TOL;
149: } else if (nsv && nconv >= nsv) {
150: PetscCall(PetscInfo(svd,"Singular value solver finished successfully: %" PetscInt_FMT " singular triplets converged at iteration %" PetscInt_FMT "\n",nconv,its));
151: *reason = SVD_CONVERGED_TOL;
152: } else if (its >= max_it) {
153: if (svd->conv == SVD_CONV_MAXIT) *reason = SVD_CONVERGED_MAXIT;
154: else {
155: *reason = SVD_DIVERGED_ITS;
156: PetscCall(PetscInfo(svd,"Singular value solver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
157: }
158: }
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: /*@
163: SVDSetWorkVecs - Sets a number of work vectors into an `SVD` object.
165: Collective
167: Input Parameters:
168: + svd - the singular value solver context
169: . nleft - number of work vectors of dimension equal to left singular vector
170: - nright - number of work vectors of dimension equal to right singular vector
172: Developer Note:
173: This is `SLEPC_EXTERN` because it may be required by user plugin `SVD`
174: implementations.
176: Level: developer
178: .seealso: [](ch:svd), `SVDSetUp()`
179: @*/
180: PetscErrorCode SVDSetWorkVecs(SVD svd,PetscInt nleft,PetscInt nright)
181: {
182: Vec t;
184: PetscFunctionBegin;
188: PetscCheck(nleft>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nleft must be >= 0: nleft = %" PetscInt_FMT,nleft);
189: PetscCheck(nright>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nright must be >= 0: nright = %" PetscInt_FMT,nright);
190: PetscCheck(nleft>0 || nright>0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"nleft and nright cannot be both zero");
191: if (svd->nworkl < nleft) {
192: PetscCall(VecDestroyVecs(svd->nworkl,&svd->workl));
193: svd->nworkl = nleft;
194: if (svd->isgeneralized) PetscCall(SVDCreateLeftTemplate(svd,&t));
195: else PetscCall(MatCreateVecsEmpty(svd->OP,NULL,&t));
196: PetscCall(VecDuplicateVecs(t,nleft,&svd->workl));
197: PetscCall(VecDestroy(&t));
198: }
199: if (svd->nworkr < nright) {
200: PetscCall(VecDestroyVecs(svd->nworkr,&svd->workr));
201: svd->nworkr = nright;
202: PetscCall(MatCreateVecsEmpty(svd->OP,&t,NULL));
203: PetscCall(VecDuplicateVecs(t,nright,&svd->workr));
204: PetscCall(VecDestroy(&t));
205: }
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }