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: SVD routines related to the solution process
12: */
14: #include <slepc/private/svdimpl.h> /*I "slepcsvd.h" I*/
16: PetscErrorCode SVDComputeVectors(SVD svd) 17: {
19: Vec tl,uj,vj;
20: PetscInt j,oldsize;
23: SVDCheckSolved(svd,1);
24: if (svd->state==SVD_STATE_SOLVED) {
25: /* generate left singular vectors on U */
26: if (!svd->U) { SVDGetBV(svd,NULL,&svd->U); }
27: BVGetSizes(svd->U,NULL,NULL,&oldsize);
28: if (!oldsize) {
29: if (!((PetscObject)(svd->U))->type_name) {
30: BVSetType(svd->U,BVSVEC);
31: }
32: SVDMatCreateVecsEmpty(svd,NULL,&tl);
33: BVSetSizesFromVec(svd->U,tl,svd->ncv);
34: VecDestroy(&tl);
35: }
36: for (j=0;j<svd->nconv;j++) {
37: BVGetColumn(svd->V,j,&vj);
38: BVGetColumn(svd->U,j,&uj);
39: SVDMatMult(svd,PETSC_FALSE,vj,uj);
40: BVRestoreColumn(svd->V,j,&vj);
41: BVRestoreColumn(svd->U,j,&uj);
42: BVOrthonormalizeColumn(svd->U,j,PETSC_FALSE,NULL,NULL);
43: }
44: }
45: svd->state = SVD_STATE_VECTORS;
46: return(0);
47: }
49: /*@
50: SVDSolve - Solves the singular value problem.
52: Collective on SVD 54: Input Parameter:
55: . svd - singular value solver context obtained from SVDCreate()
57: Options Database Keys:
58: + -svd_view - print information about the solver used
59: . -svd_view_mat binary - save the matrix to the default binary viewer
60: . -svd_view_vectors binary - save the computed singular vectors to the default binary viewer
61: . -svd_view_values - print computed singular values
62: . -svd_converged_reason - print reason for convergence, and number of iterations
63: . -svd_error_absolute - print absolute errors of each singular triplet
64: - -svd_error_relative - print relative errors of each singular triplet
66: Level: beginner
68: .seealso: SVDCreate(), SVDSetUp(), SVDDestroy()
69: @*/
70: PetscErrorCode SVDSolve(SVD svd) 71: {
73: PetscInt i,*workperm;
77: if (svd->state>=SVD_STATE_SOLVED) return(0);
78: PetscLogEventBegin(SVD_Solve,svd,0,0,0);
80: /* call setup */
81: SVDSetUp(svd);
82: svd->its = 0;
83: svd->nconv = 0;
84: for (i=0;i<svd->ncv;i++) {
85: svd->sigma[i] = 0.0;
86: svd->errest[i] = 0.0;
87: svd->perm[i] = i;
88: }
89: SVDViewFromOptions(svd,NULL,"-svd_view_pre");
91: (*svd->ops->solve)(svd);
92: svd->state = (svd->leftbasis)? SVD_STATE_VECTORS: SVD_STATE_SOLVED;
94: /* sort singular triplets */
95: if (svd->which == SVD_SMALLEST) {
96: PetscSortRealWithPermutation(svd->nconv,svd->sigma,svd->perm);
97: } else {
98: PetscMalloc1(svd->nconv,&workperm);
99: for (i=0;i<svd->nconv;i++) workperm[i] = i;
100: PetscSortRealWithPermutation(svd->nconv,svd->sigma,workperm);
101: for (i=0;i<svd->nconv;i++) svd->perm[i] = workperm[svd->nconv-i-1];
102: PetscFree(workperm);
103: }
104: PetscLogEventEnd(SVD_Solve,svd,0,0,0);
106: /* various viewers */
107: SVDViewFromOptions(svd,NULL,"-svd_view");
108: SVDReasonViewFromOptions(svd);
109: SVDErrorViewFromOptions(svd);
110: SVDValuesViewFromOptions(svd);
111: SVDVectorsViewFromOptions(svd);
112: MatViewFromOptions(svd->OP,(PetscObject)svd,"-svd_view_mat");
114: /* Remove the initial subspaces */
115: svd->nini = 0;
116: svd->ninil = 0;
117: return(0);
118: }
120: /*@
121: SVDGetIterationNumber - Gets the current iteration number. If the
122: call to SVDSolve() is complete, then it returns the number of iterations
123: carried out by the solution method.
125: Not Collective
127: Input Parameter:
128: . svd - the singular value solver context
130: Output Parameter:
131: . its - number of iterations
133: Note:
134: During the i-th iteration this call returns i-1. If SVDSolve() is
135: complete, then parameter "its" contains either the iteration number at
136: which convergence was successfully reached, or failure was detected.
137: Call SVDGetConvergedReason() to determine if the solver converged or
138: failed and why.
140: Level: intermediate
142: .seealso: SVDGetConvergedReason(), SVDSetTolerances()
143: @*/
144: PetscErrorCode SVDGetIterationNumber(SVD svd,PetscInt *its)145: {
149: *its = svd->its;
150: return(0);
151: }
153: /*@
154: SVDGetConvergedReason - Gets the reason why the SVDSolve() iteration was
155: stopped.
157: Not Collective
159: Input Parameter:
160: . svd - the singular value solver context
162: Output Parameter:
163: . reason - negative value indicates diverged, positive value converged
164: (see SVDConvergedReason)
166: Notes:
168: Possible values for reason are
169: + SVD_CONVERGED_TOL - converged up to tolerance
170: . SVD_CONVERGED_USER - converged due to a user-defined condition
171: . SVD_DIVERGED_ITS - required more than max_it iterations to reach convergence
172: - SVD_DIVERGED_BREAKDOWN - generic breakdown in method
174: Can only be called after the call to SVDSolve() is complete.
176: Level: intermediate
178: .seealso: SVDSetTolerances(), SVDSolve(), SVDConvergedReason179: @*/
180: PetscErrorCode SVDGetConvergedReason(SVD svd,SVDConvergedReason *reason)181: {
185: SVDCheckSolved(svd,1);
186: *reason = svd->reason;
187: return(0);
188: }
190: /*@
191: SVDGetConverged - Gets the number of converged singular values.
193: Not Collective
195: Input Parameter:
196: . svd - the singular value solver context
198: Output Parameter:
199: . nconv - number of converged singular values
201: Note:
202: This function should be called after SVDSolve() has finished.
204: Level: beginner
206: @*/
207: PetscErrorCode SVDGetConverged(SVD svd,PetscInt *nconv)208: {
212: SVDCheckSolved(svd,1);
213: *nconv = svd->nconv;
214: return(0);
215: }
217: /*@C
218: SVDGetSingularTriplet - Gets the i-th triplet of the singular value decomposition
219: as computed by SVDSolve(). The solution consists in the singular value and its left
220: and right singular vectors.
222: Not Collective, but vectors are shared by all processors that share the SVD224: Input Parameters:
225: + svd - singular value solver context
226: - i - index of the solution
228: Output Parameters:
229: + sigma - singular value
230: . u - left singular vector
231: - v - right singular vector
233: Note:
234: Both U or V can be NULL if singular vectors are not required.
235: Otherwise, the caller must provide valid Vec objects, i.e.,
236: they must be created by the calling program with e.g. MatCreateVecs().
238: The index i should be a value between 0 and nconv-1 (see SVDGetConverged()).
239: Singular triplets are indexed according to the ordering criterion established
240: with SVDSetWhichSingularTriplets().
242: Level: beginner
244: .seealso: SVDSolve(), SVDGetConverged(), SVDSetWhichSingularTriplets()
245: @*/
246: PetscErrorCode SVDGetSingularTriplet(SVD svd,PetscInt i,PetscReal *sigma,Vec u,Vec v)247: {
249: PetscInt M,N;
250: Vec w;
255: SVDCheckSolved(svd,1);
258: if (i<0 || i>=svd->nconv) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
259: if (sigma) *sigma = svd->sigma[svd->perm[i]];
260: MatGetSize(svd->OP,&M,&N);
261: if (M<N) { w = u; u = v; v = w; }
262: if (u) {
263: SVDComputeVectors(svd);
264: BVCopyVec(svd->U,svd->perm[i],u);
265: }
266: if (v) {
267: BVCopyVec(svd->V,svd->perm[i],v);
268: }
269: return(0);
270: }
272: /*
273: SVDComputeResidualNorms_Private - Computes the norms of the left and
274: right residuals associated with the i-th computed singular triplet.
275: @*/
276: static PetscErrorCode SVDComputeResidualNorms_Private(SVD svd,PetscInt i,PetscReal *norm1,PetscReal *norm2)277: {
279: Vec u,v,x = NULL,y = NULL;
280: PetscReal sigma;
281: PetscInt M,N;
284: MatCreateVecs(svd->OP,&v,&u);
285: SVDGetSingularTriplet(svd,i,&sigma,u,v);
286: /* norm1 = ||A*v-sigma*u||_2 */
287: if (norm1) {
288: VecDuplicate(u,&x);
289: MatMult(svd->OP,v,x);
290: VecAXPY(x,-sigma,u);
291: VecNorm(x,NORM_2,norm1);
292: }
293: /* norm2 = ||A^T*u-sigma*v||_2 */
294: if (norm2) {
295: VecDuplicate(v,&y);
296: if (svd->A && svd->AT) {
297: MatGetSize(svd->OP,&M,&N);
298: if (M<N) {
299: MatMult(svd->A,u,y);
300: } else {
301: MatMult(svd->AT,u,y);
302: }
303: } else {
304: #if defined(PETSC_USE_COMPLEX)
305: MatMultHermitianTranspose(svd->OP,u,y);
306: #else
307: MatMultTranspose(svd->OP,u,y);
308: #endif
309: }
310: VecAXPY(y,-sigma,v);
311: VecNorm(y,NORM_2,norm2);
312: }
314: VecDestroy(&v);
315: VecDestroy(&u);
316: VecDestroy(&x);
317: VecDestroy(&y);
318: return(0);
319: }
321: /*@
322: SVDComputeError - Computes the error (based on the residual norm) associated
323: with the i-th singular triplet.
325: Collective on SVD327: Input Parameter:
328: + svd - the singular value solver context
329: . i - the solution index
330: - type - the type of error to compute
332: Output Parameter:
333: . error - the error
335: Notes:
336: The error can be computed in various ways, all of them based on the residual
337: norm obtained as sqrt(n1^2+n2^2) with n1 = ||A*v-sigma*u||_2 and
338: n2 = ||A^T*u-sigma*v||_2, where sigma is the singular value, u is the left
339: singular vector and v is the right singular vector.
341: Level: beginner
343: .seealso: SVDErrorType, SVDSolve()
344: @*/
345: PetscErrorCode SVDComputeError(SVD svd,PetscInt i,SVDErrorType type,PetscReal *error)346: {
348: PetscReal sigma,norm1,norm2;
355: SVDCheckSolved(svd,1);
356: SVDGetSingularTriplet(svd,i,&sigma,NULL,NULL);
357: SVDComputeResidualNorms_Private(svd,i,&norm1,&norm2);
358: *error = PetscSqrtReal(norm1*norm1+norm2*norm2);
359: switch (type) {
360: case SVD_ERROR_ABSOLUTE:
361: break;
362: case SVD_ERROR_RELATIVE:
363: *error /= sigma;
364: break;
365: default:366: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
367: }
368: return(0);
369: }