1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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> 16: /*
17: SVDComputeVectors_Left - Compute left singular vectors as U=A*V.
18: Only done if the leftbasis flag is false. Assumes V is available.
19: */
20: PetscErrorCode SVDComputeVectors_Left(SVD svd) 21: {
23: Vec tl;
24: PetscInt oldsize;
27: if (!svd->leftbasis) {
28: /* generate left singular vectors on U */
29: if (!svd->U) { SVDGetBV(svd,NULL,&svd->U); }
30: BVGetSizes(svd->U,NULL,NULL,&oldsize);
31: if (!oldsize) {
32: if (!((PetscObject)(svd->U))->type_name) {
33: BVSetType(svd->U,BVSVEC);
34: }
35: MatCreateVecsEmpty(svd->A,NULL,&tl);
36: BVSetSizesFromVec(svd->U,tl,svd->ncv);
37: VecDestroy(&tl);
38: }
39: BVSetActiveColumns(svd->V,0,svd->nconv);
40: BVSetActiveColumns(svd->U,0,svd->nconv);
41: BVMatMult(svd->V,svd->A,svd->U);
42: BVOrthogonalize(svd->U,NULL);
43: }
44: return(0);
45: }
47: PetscErrorCode SVDComputeVectors(SVD svd) 48: {
52: SVDCheckSolved(svd,1);
53: if (svd->state==SVD_STATE_SOLVED && svd->ops->computevectors) {
54: (*svd->ops->computevectors)(svd);
55: }
56: svd->state = SVD_STATE_VECTORS;
57: return(0);
58: }
60: /*@
61: SVDSolve - Solves the singular value problem.
63: Collective on svd
65: Input Parameter:
66: . svd - singular value solver context obtained from SVDCreate()
68: Options Database Keys:
69: + -svd_view - print information about the solver used
70: . -svd_view_mat0 binary - save the first matrix (A) to the default binary viewer
71: . -svd_view_mat1 binary - save the second matrix (B) to the default binary viewer
72: . -svd_view_vectors binary - save the computed singular vectors to the default binary viewer
73: . -svd_view_values - print computed singular values
74: . -svd_converged_reason - print reason for convergence, and number of iterations
75: . -svd_error_absolute - print absolute errors of each singular triplet
76: - -svd_error_relative - print relative errors of each singular triplet
78: Level: beginner
80: .seealso: SVDCreate(), SVDSetUp(), SVDDestroy()
81: @*/
82: PetscErrorCode SVDSolve(SVD svd) 83: {
85: PetscInt i,*workperm;
89: if (svd->state>=SVD_STATE_SOLVED) return(0);
90: PetscLogEventBegin(SVD_Solve,svd,0,0,0);
92: /* call setup */
93: SVDSetUp(svd);
94: svd->its = 0;
95: svd->nconv = 0;
96: for (i=0;i<svd->ncv;i++) {
97: svd->sigma[i] = 0.0;
98: svd->errest[i] = 0.0;
99: svd->perm[i] = i;
100: }
101: SVDViewFromOptions(svd,NULL,"-svd_view_pre");
103: switch (svd->problem_type) {
104: case SVD_STANDARD:
105: (*svd->ops->solve)(svd);
106: break;
107: case SVD_GENERALIZED:
108: (*svd->ops->solveg)(svd);
109: break;
110: }
111: svd->state = SVD_STATE_SOLVED;
113: /* sort singular triplets */
114: if (svd->which == SVD_SMALLEST) {
115: PetscSortRealWithPermutation(svd->nconv,svd->sigma,svd->perm);
116: } else {
117: PetscMalloc1(svd->nconv,&workperm);
118: for (i=0;i<svd->nconv;i++) workperm[i] = i;
119: PetscSortRealWithPermutation(svd->nconv,svd->sigma,workperm);
120: for (i=0;i<svd->nconv;i++) svd->perm[i] = workperm[svd->nconv-i-1];
121: PetscFree(workperm);
122: }
123: PetscLogEventEnd(SVD_Solve,svd,0,0,0);
125: /* various viewers */
126: SVDViewFromOptions(svd,NULL,"-svd_view");
127: SVDConvergedReasonViewFromOptions(svd);
128: SVDErrorViewFromOptions(svd);
129: SVDValuesViewFromOptions(svd);
130: SVDVectorsViewFromOptions(svd);
131: MatViewFromOptions(svd->OP,(PetscObject)svd,"-svd_view_mat0");
132: if (svd->isgeneralized) {
133: MatViewFromOptions(svd->OPb,(PetscObject)svd,"-svd_view_mat1");
134: }
136: /* Remove the initial subspaces */
137: svd->nini = 0;
138: svd->ninil = 0;
139: return(0);
140: }
142: /*@
143: SVDGetIterationNumber - Gets the current iteration number. If the
144: call to SVDSolve() is complete, then it returns the number of iterations
145: carried out by the solution method.
147: Not Collective
149: Input Parameter:
150: . svd - the singular value solver context
152: Output Parameter:
153: . its - number of iterations
155: Note:
156: During the i-th iteration this call returns i-1. If SVDSolve() is
157: complete, then parameter "its" contains either the iteration number at
158: which convergence was successfully reached, or failure was detected.
159: Call SVDGetConvergedReason() to determine if the solver converged or
160: failed and why.
162: Level: intermediate
164: .seealso: SVDGetConvergedReason(), SVDSetTolerances()
165: @*/
166: PetscErrorCode SVDGetIterationNumber(SVD svd,PetscInt *its)167: {
171: *its = svd->its;
172: return(0);
173: }
175: /*@
176: SVDGetConvergedReason - Gets the reason why the SVDSolve() iteration was
177: stopped.
179: Not Collective
181: Input Parameter:
182: . svd - the singular value solver context
184: Output Parameter:
185: . reason - negative value indicates diverged, positive value converged
186: (see SVDConvergedReason)
188: Notes:
190: Possible values for reason are
191: + SVD_CONVERGED_TOL - converged up to tolerance
192: . SVD_CONVERGED_USER - converged due to a user-defined condition
193: . SVD_CONVERGED_MAXIT - reached the maximum number of iterations with SVD_CONV_MAXIT criterion
194: . SVD_DIVERGED_ITS - required more than max_it iterations to reach convergence
195: - SVD_DIVERGED_BREAKDOWN - generic breakdown in method
197: Can only be called after the call to SVDSolve() is complete.
199: Level: intermediate
201: .seealso: SVDSetTolerances(), SVDSolve(), SVDConvergedReason202: @*/
203: PetscErrorCode SVDGetConvergedReason(SVD svd,SVDConvergedReason *reason)204: {
208: SVDCheckSolved(svd,1);
209: *reason = svd->reason;
210: return(0);
211: }
213: /*@
214: SVDGetConverged - Gets the number of converged singular values.
216: Not Collective
218: Input Parameter:
219: . svd - the singular value solver context
221: Output Parameter:
222: . nconv - number of converged singular values
224: Note:
225: This function should be called after SVDSolve() has finished.
227: Level: beginner
229: @*/
230: PetscErrorCode SVDGetConverged(SVD svd,PetscInt *nconv)231: {
235: SVDCheckSolved(svd,1);
236: *nconv = svd->nconv;
237: return(0);
238: }
240: /*@C
241: SVDGetSingularTriplet - Gets the i-th triplet of the singular value decomposition
242: as computed by SVDSolve(). The solution consists in the singular value and its left
243: and right singular vectors.
245: Not Collective, but vectors are shared by all processors that share the SVD247: Input Parameters:
248: + svd - singular value solver context
249: - i - index of the solution
251: Output Parameters:
252: + sigma - singular value
253: . u - left singular vector
254: - v - right singular vector
256: Note:
257: Both u or v can be NULL if singular vectors are not required.
258: Otherwise, the caller must provide valid Vec objects, i.e.,
259: they must be created by the calling program with e.g. MatCreateVecs().
261: The index i should be a value between 0 and nconv-1 (see SVDGetConverged()).
262: Singular triplets are indexed according to the ordering criterion established
263: with SVDSetWhichSingularTriplets().
265: In the case of GSVD, the solution consists in three vectors u,v,x that are
266: returned as follows. Vector x is returned in the right singular vector
267: (argument v) and has length equal to the number of columns of A and B.
268: The other two vectors are returned stacked on top of each other [u;v] in
269: the left singular vector argument, with length equal to m+n (number of rows
270: of A plus number of rows of B).
272: Level: beginner
274: .seealso: SVDSolve(), SVDGetConverged(), SVDSetWhichSingularTriplets()
275: @*/
276: PetscErrorCode SVDGetSingularTriplet(SVD svd,PetscInt i,PetscReal *sigma,Vec u,Vec v)277: {
279: PetscInt M,N;
280: Vec w;
285: SVDCheckSolved(svd,1);
288: if (i<0) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
289: if (i>=svd->nconv) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see SVDGetConverged()");
290: if (sigma) *sigma = svd->sigma[svd->perm[i]];
291: if (u || v) {
292: if (!svd->isgeneralized) {
293: MatGetSize(svd->OP,&M,&N);
294: if (M<N) { w = u; u = v; v = w; }
295: SVDComputeVectors(svd);
296: }
297: if (u) { BVCopyVec(svd->U,svd->perm[i],u); }
298: if (v) { BVCopyVec(svd->V,svd->perm[i],v); }
299: }
300: return(0);
301: }
303: /*
304: SVDComputeResidualNorms_Standard - Computes the norms of the left and
305: right residuals associated with the i-th computed singular triplet.
307: Input Parameters:
308: sigma - singular value
309: u,v - singular vectors
310: x,y - two work vectors with the same dimensions as u,v
311: @*/
312: static PetscErrorCode SVDComputeResidualNorms_Standard(SVD svd,PetscReal sigma,Vec u,Vec v,Vec x,Vec y,PetscReal *norm1,PetscReal *norm2)313: {
315: PetscInt M,N;
318: /* norm1 = ||A*v-sigma*u||_2 */
319: if (norm1) {
320: MatMult(svd->OP,v,x);
321: VecAXPY(x,-sigma,u);
322: VecNorm(x,NORM_2,norm1);
323: }
324: /* norm2 = ||A^T*u-sigma*v||_2 */
325: if (norm2) {
326: MatGetSize(svd->OP,&M,&N);
327: if (M<N) {
328: MatMult(svd->A,u,y);
329: } else {
330: MatMult(svd->AT,u,y);
331: }
332: VecAXPY(y,-sigma,v);
333: VecNorm(y,NORM_2,norm2);
334: }
335: return(0);
336: }
338: /*
339: SVDComputeResidualNorms_Generalized - In GSVD, compute the residual norms
340: norm1 = ||A*x-c*u||_2 and norm2 = ||B*x-s*v||_2.
342: Input Parameters:
343: sigma - singular value
344: x - right singular vector
345: */
346: static PetscErrorCode SVDComputeResidualNorms_Generalized(SVD svd,PetscReal sigma,Vec uv,Vec x,PetscReal *norm1,PetscReal *norm2)347: {
349: Vec u,v,ax,bx,nest,aux[2];
350: PetscReal c,s;
353: MatCreateVecs(svd->OP,NULL,&u);
354: MatCreateVecs(svd->OPb,NULL,&v);
355: aux[0] = u;
356: aux[1] = v;
357: VecCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,aux,&nest);
358: VecCopy(uv,nest);
359: VecDuplicate(u,&ax);
360: VecDuplicate(v,&bx);
362: s = 1.0/PetscSqrtReal(1.0+sigma*sigma);
363: c = sigma*s;
365: /* norm1 = ||A*x-c*u||_2 */
366: if (norm1) {
367: MatMult(svd->OP,x,ax);
368: VecAXPY(ax,-c,u);
369: VecNorm(ax,NORM_2,norm1);
370: }
371: /* norm2 = ||B*x-s*v||_2 */
372: if (norm2) {
373: MatMult(svd->OPb,x,bx);
374: VecAXPY(bx,-s,v);
375: VecNorm(bx,NORM_2,norm2);
376: }
378: VecDestroy(&v);
379: VecDestroy(&u);
380: VecDestroy(&nest);
381: VecDestroy(&ax);
382: VecDestroy(&bx);
383: return(0);
384: }
386: /*@
387: SVDComputeError - Computes the error (based on the residual norm) associated
388: with the i-th singular triplet.
390: Collective on svd
392: Input Parameter:
393: + svd - the singular value solver context
394: . i - the solution index
395: - type - the type of error to compute
397: Output Parameter:
398: . error - the error
400: Notes:
401: The error can be computed in various ways, all of them based on the residual
402: norm obtained as sqrt(n1^2+n2^2) with n1 = ||A*v-sigma*u||_2 and
403: n2 = ||A^T*u-sigma*v||_2, where sigma is the singular value, u is the left
404: singular vector and v is the right singular vector.
406: Level: beginner
408: .seealso: SVDErrorType, SVDSolve()
409: @*/
410: PetscErrorCode SVDComputeError(SVD svd,PetscInt i,SVDErrorType type,PetscReal *error)411: {
413: PetscReal sigma,norm1,norm2;
414: Vec u,v,x,y;
421: SVDCheckSolved(svd,1);
423: /* allocate work vectors */
424: SVDSetWorkVecs(svd,2,2);
425: u = svd->workl[0];
426: v = svd->workr[0];
427: x = svd->workl[1];
428: y = svd->workr[1];
430: /* compute residual norm and error */
431: SVDGetSingularTriplet(svd,i,&sigma,u,v);
432: switch (svd->problem_type) {
433: case SVD_STANDARD:
434: SVDComputeResidualNorms_Standard(svd,sigma,u,v,x,y,&norm1,&norm2);
435: break;
436: case SVD_GENERALIZED:
437: SVDComputeResidualNorms_Generalized(svd,sigma,u,v,&norm1,&norm2);
438: break;
439: }
440: *error = PetscSqrtReal(norm1*norm1+norm2*norm2);
441: switch (type) {
442: case SVD_ERROR_ABSOLUTE:
443: break;
444: case SVD_ERROR_RELATIVE:
445: *error /= sigma;
446: break;
447: default:448: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
449: }
450: return(0);
451: }