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: 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: {
22: Vec tl;
23: PetscInt oldsize;
25: if (!svd->leftbasis) {
26: /* generate left singular vectors on U */
27: if (!svd->U) SVDGetBV(svd,NULL,&svd->U);
28: BVGetSizes(svd->U,NULL,NULL,&oldsize);
29: if (!oldsize) {
30: if (!((PetscObject)(svd->U))->type_name) BVSetType(svd->U,((PetscObject)(svd->V))->type_name);
31: MatCreateVecsEmpty(svd->A,NULL,&tl);
32: BVSetSizesFromVec(svd->U,tl,svd->ncv);
33: VecDestroy(&tl);
34: }
35: BVSetActiveColumns(svd->V,0,svd->nconv);
36: BVSetActiveColumns(svd->U,0,svd->nconv);
37: BVMatMult(svd->V,svd->A,svd->U);
38: BVOrthogonalize(svd->U,NULL);
39: }
40: PetscFunctionReturn(0);
41: }
43: PetscErrorCode SVDComputeVectors(SVD svd) 44: {
45: SVDCheckSolved(svd,1);
46: if (svd->state==SVD_STATE_SOLVED && svd->ops->computevectors) (*svd->ops->computevectors)(svd);
47: svd->state = SVD_STATE_VECTORS;
48: PetscFunctionReturn(0);
49: }
51: /*@
52: SVDSolve - Solves the singular value problem.
54: Collective on svd
56: Input Parameter:
57: . svd - singular value solver context obtained from SVDCreate()
59: Options Database Keys:
60: + -svd_view - print information about the solver used
61: . -svd_view_mat0 - view the first matrix (A)
62: . -svd_view_mat1 - view the second matrix (B)
63: . -svd_view_vectors - view the computed singular vectors
64: . -svd_view_values - view the computed singular values
65: . -svd_converged_reason - print reason for convergence, and number of iterations
66: . -svd_error_absolute - print absolute errors of each singular triplet
67: . -svd_error_relative - print relative errors of each singular triplet
68: - -svd_error_norm - print errors relative to the matrix norms of each singular triplet
70: Notes:
71: All the command-line options listed above admit an optional argument specifying
72: the viewer type and options. For instance, use '-svd_view_mat0 binary:amatrix.bin'
73: to save the A matrix to a binary file, '-svd_view_values draw' to draw the computed
74: singular values graphically, or '-svd_error_relative :myerr.m:ascii_matlab' to save
75: the errors in a file that can be executed in Matlab.
77: Level: beginner
79: .seealso: SVDCreate(), SVDSetUp(), SVDDestroy()
80: @*/
81: PetscErrorCode SVDSolve(SVD svd) 82: {
83: PetscInt i,*workperm;
86: if (svd->state>=SVD_STATE_SOLVED) PetscFunctionReturn(0);
87: PetscLogEventBegin(SVD_Solve,svd,0,0,0);
89: /* call setup */
90: SVDSetUp(svd);
91: svd->its = 0;
92: svd->nconv = 0;
93: for (i=0;i<svd->ncv;i++) {
94: svd->sigma[i] = 0.0;
95: svd->errest[i] = 0.0;
96: svd->perm[i] = i;
97: }
98: SVDViewFromOptions(svd,NULL,"-svd_view_pre");
100: switch (svd->problem_type) {
101: case SVD_STANDARD:
102: (*svd->ops->solve)(svd);
103: break;
104: case SVD_GENERALIZED:
105: (*svd->ops->solveg)(svd);
106: break;
107: }
108: svd->state = SVD_STATE_SOLVED;
110: /* sort singular triplets */
111: if (svd->which == SVD_SMALLEST) PetscSortRealWithPermutation(svd->nconv,svd->sigma,svd->perm);
112: else {
113: PetscMalloc1(svd->nconv,&workperm);
114: for (i=0;i<svd->nconv;i++) workperm[i] = i;
115: PetscSortRealWithPermutation(svd->nconv,svd->sigma,workperm);
116: for (i=0;i<svd->nconv;i++) svd->perm[i] = workperm[svd->nconv-i-1];
117: PetscFree(workperm);
118: }
119: PetscLogEventEnd(SVD_Solve,svd,0,0,0);
121: /* various viewers */
122: SVDViewFromOptions(svd,NULL,"-svd_view");
123: SVDConvergedReasonViewFromOptions(svd);
124: SVDErrorViewFromOptions(svd);
125: SVDValuesViewFromOptions(svd);
126: SVDVectorsViewFromOptions(svd);
127: MatViewFromOptions(svd->OP,(PetscObject)svd,"-svd_view_mat0");
128: if (svd->isgeneralized) MatViewFromOptions(svd->OPb,(PetscObject)svd,"-svd_view_mat1");
130: /* Remove the initial subspaces */
131: svd->nini = 0;
132: svd->ninil = 0;
133: PetscFunctionReturn(0);
134: }
136: /*@
137: SVDGetIterationNumber - Gets the current iteration number. If the
138: call to SVDSolve() is complete, then it returns the number of iterations
139: carried out by the solution method.
141: Not Collective
143: Input Parameter:
144: . svd - the singular value solver context
146: Output Parameter:
147: . its - number of iterations
149: Note:
150: During the i-th iteration this call returns i-1. If SVDSolve() is
151: complete, then parameter "its" contains either the iteration number at
152: which convergence was successfully reached, or failure was detected.
153: Call SVDGetConvergedReason() to determine if the solver converged or
154: failed and why.
156: Level: intermediate
158: .seealso: SVDGetConvergedReason(), SVDSetTolerances()
159: @*/
160: PetscErrorCode SVDGetIterationNumber(SVD svd,PetscInt *its)161: {
164: *its = svd->its;
165: PetscFunctionReturn(0);
166: }
168: /*@
169: SVDGetConvergedReason - Gets the reason why the SVDSolve() iteration was
170: stopped.
172: Not Collective
174: Input Parameter:
175: . svd - the singular value solver context
177: Output Parameter:
178: . reason - negative value indicates diverged, positive value converged
179: (see SVDConvergedReason)
181: Options Database Key:
182: . -svd_converged_reason - print the reason to a viewer
184: Notes:
185: Possible values for reason are
186: + SVD_CONVERGED_TOL - converged up to tolerance
187: . SVD_CONVERGED_USER - converged due to a user-defined condition
188: . SVD_CONVERGED_MAXIT - reached the maximum number of iterations with SVD_CONV_MAXIT criterion
189: . SVD_DIVERGED_ITS - required more than max_it iterations to reach convergence
190: - SVD_DIVERGED_BREAKDOWN - generic breakdown in method
192: Can only be called after the call to SVDSolve() is complete.
194: Level: intermediate
196: .seealso: SVDSetTolerances(), SVDSolve(), SVDConvergedReason197: @*/
198: PetscErrorCode SVDGetConvergedReason(SVD svd,SVDConvergedReason *reason)199: {
202: SVDCheckSolved(svd,1);
203: *reason = svd->reason;
204: PetscFunctionReturn(0);
205: }
207: /*@
208: SVDGetConverged - Gets the number of converged singular values.
210: Not Collective
212: Input Parameter:
213: . svd - the singular value solver context
215: Output Parameter:
216: . nconv - number of converged singular values
218: Note:
219: This function should be called after SVDSolve() has finished.
221: Level: beginner
223: .seealso: SVDSetDimensions(), SVDSolve(), SVDGetSingularTriplet()
224: @*/
225: PetscErrorCode SVDGetConverged(SVD svd,PetscInt *nconv)226: {
229: SVDCheckSolved(svd,1);
230: *nconv = svd->nconv;
231: PetscFunctionReturn(0);
232: }
234: /*@C
235: SVDGetSingularTriplet - Gets the i-th triplet of the singular value decomposition
236: as computed by SVDSolve(). The solution consists in the singular value and its left
237: and right singular vectors.
239: Not Collective, but vectors are shared by all processors that share the SVD241: Input Parameters:
242: + svd - singular value solver context
243: - i - index of the solution
245: Output Parameters:
246: + sigma - singular value
247: . u - left singular vector
248: - v - right singular vector
250: Note:
251: Both u or v can be NULL if singular vectors are not required.
252: Otherwise, the caller must provide valid Vec objects, i.e.,
253: they must be created by the calling program with e.g. MatCreateVecs().
255: The index i should be a value between 0 and nconv-1 (see SVDGetConverged()).
256: Singular triplets are indexed according to the ordering criterion established
257: with SVDSetWhichSingularTriplets().
259: In the case of GSVD, the solution consists in three vectors u,v,x that are
260: returned as follows. Vector x is returned in the right singular vector
261: (argument v) and has length equal to the number of columns of A and B.
262: The other two vectors are returned stacked on top of each other [u;v] in
263: the left singular vector argument, with length equal to m+n (number of rows
264: of A plus number of rows of B).
266: Level: beginner
268: .seealso: SVDSolve(), SVDGetConverged(), SVDSetWhichSingularTriplets()
269: @*/
270: PetscErrorCode SVDGetSingularTriplet(SVD svd,PetscInt i,PetscReal *sigma,Vec u,Vec v)271: {
272: PetscInt M,N;
273: Vec w;
277: SVDCheckSolved(svd,1);
282: if (sigma) *sigma = svd->sigma[svd->perm[i]];
283: if (u || v) {
284: if (!svd->isgeneralized) {
285: MatGetSize(svd->OP,&M,&N);
286: if (M<N) { w = u; u = v; v = w; }
287: }
288: SVDComputeVectors(svd);
289: if (u) BVCopyVec(svd->U,svd->perm[i],u);
290: if (v) BVCopyVec(svd->V,svd->perm[i],v);
291: }
292: PetscFunctionReturn(0);
293: }
295: /*
296: SVDComputeResidualNorms_Standard - Computes the norms of the left and
297: right residuals associated with the i-th computed singular triplet.
299: Input Parameters:
300: sigma - singular value
301: u,v - singular vectors
302: x,y - two work vectors with the same dimensions as u,v
303: */
304: static PetscErrorCode SVDComputeResidualNorms_Standard(SVD svd,PetscReal sigma,Vec u,Vec v,Vec x,Vec y,PetscReal *norm1,PetscReal *norm2)305: {
306: PetscInt M,N;
308: /* norm1 = ||A*v-sigma*u||_2 */
309: if (norm1) {
310: MatMult(svd->OP,v,x);
311: VecAXPY(x,-sigma,u);
312: VecNorm(x,NORM_2,norm1);
313: }
314: /* norm2 = ||A^T*u-sigma*v||_2 */
315: if (norm2) {
316: MatGetSize(svd->OP,&M,&N);
317: if (M<N) MatMult(svd->A,u,y);
318: else MatMult(svd->AT,u,y);
319: VecAXPY(y,-sigma,v);
320: VecNorm(y,NORM_2,norm2);
321: }
322: PetscFunctionReturn(0);
323: }
325: /*
326: SVDComputeResidualNorms_Generalized - In GSVD, compute the residual norms
327: norm1 = ||s^2*A'*u-c*B'*B*x||_2 and norm2 = ||c^2*B'*v-s*A'*A*x||_2.
329: Input Parameters:
330: sigma - singular value
331: uv - left singular vectors [u;v]
332: x - right singular vector
333: y,z - two work vectors with the same dimension as x
334: */
335: static PetscErrorCode SVDComputeResidualNorms_Generalized(SVD svd,PetscReal sigma,Vec uv,Vec x,Vec y,Vec z,PetscReal *norm1,PetscReal *norm2)336: {
337: Vec u,v,ax,bx,nest,aux[2];
338: PetscReal c,s;
340: MatCreateVecs(svd->OP,NULL,&u);
341: MatCreateVecs(svd->OPb,NULL,&v);
342: aux[0] = u;
343: aux[1] = v;
344: VecCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,aux,&nest);
345: VecCopy(uv,nest);
347: s = 1.0/PetscSqrtReal(1.0+sigma*sigma);
348: c = sigma*s;
350: /* norm1 = ||s^2*A'*u-c*B'*B*x||_2 */
351: if (norm1) {
352: VecDuplicate(v,&bx);
353: MatMultHermitianTranspose(svd->OP,u,z);
354: MatMult(svd->OPb,x,bx);
355: MatMultHermitianTranspose(svd->OPb,bx,y);
356: VecAXPBY(y,s*s,-c,z);
357: VecNorm(y,NORM_2,norm1);
358: VecDestroy(&bx);
359: }
360: /* norm2 = ||c^2*B'*v-s*A'*A*x||_2 */
361: if (norm2) {
362: VecDuplicate(u,&ax);
363: MatMultHermitianTranspose(svd->OPb,v,z);
364: MatMult(svd->OP,x,ax);
365: MatMultHermitianTranspose(svd->OP,ax,y);
366: VecAXPBY(y,c*c,-s,z);
367: VecNorm(y,NORM_2,norm2);
368: VecDestroy(&ax);
369: }
371: VecDestroy(&v);
372: VecDestroy(&u);
373: VecDestroy(&nest);
374: PetscFunctionReturn(0);
375: }
377: /*@
378: SVDComputeError - Computes the error (based on the residual norm) associated
379: with the i-th singular triplet.
381: Collective on svd
383: Input Parameters:
384: + svd - the singular value solver context
385: . i - the solution index
386: - type - the type of error to compute
388: Output Parameter:
389: . error - the error
391: Notes:
392: The error can be computed in various ways, all of them based on the residual
393: norm obtained as sqrt(n1^2+n2^2) with n1 = ||A*v-sigma*u||_2 and
394: n2 = ||A^T*u-sigma*v||_2, where sigma is the singular value, u is the left
395: singular vector and v is the right singular vector.
397: In the case of the GSVD, the two components of the residual norm are
398: n1 = ||s^2*A'*u-c*B'*B*x||_2 and n2 = ||c^2*B'*v-s*A'*A*x||_2, where [u;v]
399: are the left singular vectors and x is the right singular vector, with
400: sigma=c/s.
402: Level: beginner
404: .seealso: SVDErrorType, SVDSolve()
405: @*/
406: PetscErrorCode SVDComputeError(SVD svd,PetscInt i,SVDErrorType type,PetscReal *error)407: {
408: PetscReal sigma,norm1,norm2;
409: Vec u=NULL,v=NULL,x=NULL,y=NULL;
415: SVDCheckSolved(svd,1);
417: /* allocate work vectors */
418: switch (svd->problem_type) {
419: case SVD_STANDARD:
420: SVDSetWorkVecs(svd,2,2);
421: u = svd->workl[0];
422: v = svd->workr[0];
423: x = svd->workl[1];
424: y = svd->workr[1];
425: break;
426: case SVD_GENERALIZED:
428: SVDSetWorkVecs(svd,1,3);
429: u = svd->workl[0];
430: v = svd->workr[0];
431: x = svd->workr[1];
432: y = svd->workr[2];
433: break;
434: }
436: /* compute residual norm and error */
437: SVDGetSingularTriplet(svd,i,&sigma,u,v);
438: switch (svd->problem_type) {
439: case SVD_STANDARD:
440: SVDComputeResidualNorms_Standard(svd,sigma,u,v,x,y,&norm1,&norm2);
441: break;
442: case SVD_GENERALIZED:
443: SVDComputeResidualNorms_Generalized(svd,sigma,u,v,x,y,&norm1,&norm2);
444: break;
445: }
446: *error = PetscSqrtReal(norm1*norm1+norm2*norm2);
447: switch (type) {
448: case SVD_ERROR_ABSOLUTE:
449: break;
450: case SVD_ERROR_RELATIVE:
451: *error /= sigma;
452: break;
453: case SVD_ERROR_NORM:
454: if (!svd->nrma) MatNorm(svd->OP,NORM_INFINITY,&svd->nrma);
455: if (svd->isgeneralized && !svd->nrmb) MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb);
456: *error /= PetscMax(svd->nrma,svd->nrmb);
457: break;
458: default:459: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
460: }
461: PetscFunctionReturn(0);
462: }