Actual source code: svdsolve.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: SVD routines related to the solution process
12: */
14: #include <slepc/private/svdimpl.h>
15: #include <slepc/private/bvimpl.h>
17: /*
18: SVDComputeVectors_Left - Compute left singular vectors as U=A*V.
19: Only done if the leftbasis flag is false. Assumes V is available.
20: */
21: PetscErrorCode SVDComputeVectors_Left(SVD svd)
22: {
23: Vec tl,omega2,u,v,w;
24: PetscInt i,oldsize;
25: VecType vtype;
26: const PetscScalar* varray;
28: PetscFunctionBegin;
29: if (!svd->leftbasis) {
30: /* generate left singular vectors on U */
31: if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
32: PetscCall(BVGetSizes(svd->U,NULL,NULL,&oldsize));
33: if (!oldsize) {
34: if (!((PetscObject)svd->U)->type_name) PetscCall(BVSetType(svd->U,((PetscObject)svd->V)->type_name));
35: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&tl));
36: PetscCall(BVSetSizesFromVec(svd->U,tl,svd->ncv));
37: PetscCall(VecDestroy(&tl));
38: }
39: PetscCall(BVSetActiveColumns(svd->V,0,svd->nconv));
40: PetscCall(BVSetActiveColumns(svd->U,0,svd->nconv));
41: if (!svd->ishyperbolic) PetscCall(BVMatMult(svd->V,svd->A,svd->U));
42: else if (svd->swapped) { /* compute right singular vectors as V=A'*Omega*U */
43: PetscCall(MatCreateVecs(svd->A,&w,NULL));
44: for (i=0;i<svd->nconv;i++) {
45: PetscCall(BVGetColumn(svd->V,i,&v));
46: PetscCall(BVGetColumn(svd->U,i,&u));
47: PetscCall(VecPointwiseMult(w,v,svd->omega));
48: PetscCall(MatMult(svd->A,w,u));
49: PetscCall(BVRestoreColumn(svd->V,i,&v));
50: PetscCall(BVRestoreColumn(svd->U,i,&u));
51: }
52: PetscCall(VecDestroy(&w));
53: } else { /* compute left singular vectors as usual U=A*V, and set-up Omega-orthogonalization of U */
54: PetscCall(BV_SetMatrixDiagonal(svd->U,svd->omega,svd->A));
55: PetscCall(BVMatMult(svd->V,svd->A,svd->U));
56: }
57: PetscCall(BVOrthogonalize(svd->U,NULL));
58: if (svd->ishyperbolic && !svd->swapped) { /* store signature after Omega-orthogonalization */
59: PetscCall(MatGetVecType(svd->A,&vtype));
60: PetscCall(VecCreate(PETSC_COMM_SELF,&omega2));
61: PetscCall(VecSetSizes(omega2,svd->nconv,svd->nconv));
62: PetscCall(VecSetType(omega2,vtype));
63: PetscCall(BVGetSignature(svd->U,omega2));
64: PetscCall(VecGetArrayRead(omega2,&varray));
65: for (i=0;i<svd->nconv;i++) {
66: svd->sign[i] = PetscRealPart(varray[i]);
67: if (PetscRealPart(varray[i])<0.0) PetscCall(BVScaleColumn(svd->U,i,-1.0));
68: }
69: PetscCall(VecRestoreArrayRead(omega2,&varray));
70: PetscCall(VecDestroy(&omega2));
71: }
72: }
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: PetscErrorCode SVDComputeVectors(SVD svd)
77: {
78: PetscFunctionBegin;
79: SVDCheckSolved(svd,1);
80: if (svd->state==SVD_STATE_SOLVED) PetscTryTypeMethod(svd,computevectors);
81: svd->state = SVD_STATE_VECTORS;
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: /*@
86: SVDSolve - Solves the singular value problem.
88: Collective
90: Input Parameter:
91: . svd - the singular value solver context
93: Options Database Keys:
94: + -svd_view - print information about the solver used
95: . -svd_view_mat0 - view the first matrix ($A$)
96: . -svd_view_mat1 - view the second matrix ($B$)
97: . -svd_view_signature - view the signature matrix ($\Omega$)
98: . -svd_view_vectors - view the computed singular vectors
99: . -svd_view_values - view the computed singular values
100: . -svd_converged_reason - print reason for convergence/divergence, and number of iterations
101: . -svd_error_absolute - print absolute errors of each singular triplet
102: . -svd_error_relative - print relative errors of each singular triplet
103: - -svd_error_norm - print errors relative to the matrix norms of each singular triplet
105: Notes:
106: The problem matrices are specified with `SVDSetOperators()`.
108: `SVDSolve()` will return without generating an error regardless of whether
109: all requested solutions were computed or not. Call `SVDGetConverged()` to get the
110: actual number of computed solutions, and `SVDGetConvergedReason()` to determine if
111: the solver converged or failed and why.
113: All the command-line options listed above admit an optional argument specifying
114: the viewer type and options. For instance, use `-svd_view_mat0 binary:amatrix.bin`
115: to save the $A$ matrix to a binary file, `-svd_view_values draw` to draw the computed
116: singular values graphically, or `-svd_error_relative :myerr.m:ascii_matlab` to save
117: the errors in a file that can be executed in Matlab.
118: See `PetscObjectViewFromOptions()` for more details.
120: Level: beginner
122: .seealso: [](ch:svd), `SVDCreate()`, `SVDSetUp()`, `SVDDestroy()`, `SVDSetOperators()`, `SVDGetConverged()`, `SVDGetConvergedReason()`
123: @*/
124: PetscErrorCode SVDSolve(SVD svd)
125: {
126: PetscInt i,m,n,*workperm;
128: PetscFunctionBegin;
130: if (svd->state>=SVD_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
131: PetscCall(PetscLogEventBegin(SVD_Solve,svd,0,0,0));
133: /* call setup */
134: PetscCall(SVDSetUp(svd));
136: /* safeguard for matrices with zero rows or columns */
137: PetscCall(MatGetSize(svd->OP,&m,&n));
138: if (m == 0 || n == 0) {
139: svd->nconv = 0;
140: svd->reason = SVD_CONVERGED_TOL;
141: svd->state = SVD_STATE_SOLVED;
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: svd->its = 0;
146: svd->nconv = 0;
147: for (i=0;i<svd->ncv;i++) {
148: svd->sigma[i] = 0.0;
149: svd->errest[i] = 0.0;
150: svd->perm[i] = i;
151: }
152: PetscCall(SVDViewFromOptions(svd,NULL,"-svd_view_pre"));
154: switch (svd->problem_type) {
155: case SVD_STANDARD:
156: PetscUseTypeMethod(svd,solve);
157: break;
158: case SVD_GENERALIZED:
159: PetscUseTypeMethod(svd,solveg);
160: break;
161: case SVD_HYPERBOLIC:
162: PetscUseTypeMethod(svd,solveh);
163: break;
164: }
165: svd->state = SVD_STATE_SOLVED;
167: /* sort singular triplets */
168: if (svd->which == SVD_SMALLEST) PetscCall(PetscSortRealWithPermutation(svd->nconv,svd->sigma,svd->perm));
169: else {
170: PetscCall(PetscMalloc1(svd->nconv,&workperm));
171: for (i=0;i<svd->nconv;i++) workperm[i] = i;
172: PetscCall(PetscSortRealWithPermutation(svd->nconv,svd->sigma,workperm));
173: for (i=0;i<svd->nconv;i++) svd->perm[i] = workperm[svd->nconv-i-1];
174: PetscCall(PetscFree(workperm));
175: }
176: PetscCall(PetscLogEventEnd(SVD_Solve,svd,0,0,0));
178: /* various viewers */
179: PetscCall(SVDViewFromOptions(svd,NULL,"-svd_view"));
180: PetscCall(SVDConvergedReasonViewFromOptions(svd));
181: PetscCall(SVDErrorViewFromOptions(svd));
182: PetscCall(SVDValuesViewFromOptions(svd));
183: PetscCall(SVDVectorsViewFromOptions(svd));
184: PetscCall(MatViewFromOptions(svd->OP,(PetscObject)svd,"-svd_view_mat0"));
185: if (svd->isgeneralized) PetscCall(MatViewFromOptions(svd->OPb,(PetscObject)svd,"-svd_view_mat1"));
186: if (svd->ishyperbolic) PetscCall(VecViewFromOptions(svd->omega,(PetscObject)svd,"-svd_view_signature"));
188: /* Remove the initial subspaces */
189: svd->nini = 0;
190: svd->ninil = 0;
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: /*@
195: SVDGetIterationNumber - Gets the current iteration number. If the
196: call to `SVDSolve()` is complete, then it returns the number of iterations
197: carried out by the solution method.
199: Not Collective
201: Input Parameter:
202: . svd - the singular value solver context
204: Output Parameter:
205: . its - number of iterations
207: Note:
208: During the $i$-th iteration this call returns $i-1$. If `SVDSolve()` is
209: complete, then parameter `its` contains either the iteration number at
210: which convergence was successfully reached, or failure was detected.
211: Call `SVDGetConvergedReason()` to determine if the solver converged or
212: failed and why.
214: Level: intermediate
216: .seealso: [](ch:svd), `SVDGetConvergedReason()`, `SVDSetTolerances()`
217: @*/
218: PetscErrorCode SVDGetIterationNumber(SVD svd,PetscInt *its)
219: {
220: PetscFunctionBegin;
222: PetscAssertPointer(its,2);
223: *its = svd->its;
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: /*@
228: SVDGetConvergedReason - Gets the reason why the `SVDSolve()` iteration was
229: stopped.
231: Not Collective
233: Input Parameter:
234: . svd - the singular value solver context
236: Output Parameter:
237: . reason - negative value indicates diverged, positive value converged, see
238: `SVDConvergedReason` for the possible values
240: Options Database Key:
241: . -svd_converged_reason - print reason for convergence/divergence, and number of iterations
243: Note:
244: If this routine is called before or doing the `SVDSolve()` the value of
245: `SVD_CONVERGED_ITERATING` is returned.
247: Level: intermediate
249: .seealso: [](ch:svd), `SVDSetTolerances()`, `SVDSolve()`, `SVDConvergedReason`
250: @*/
251: PetscErrorCode SVDGetConvergedReason(SVD svd,SVDConvergedReason *reason)
252: {
253: PetscFunctionBegin;
255: PetscAssertPointer(reason,2);
256: SVDCheckSolved(svd,1);
257: *reason = svd->reason;
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: /*@
262: SVDGetConverged - Gets the number of converged singular values.
264: Not Collective
266: Input Parameter:
267: . svd - the singular value solver context
269: Output Parameter:
270: . nconv - number of converged singular values
272: Notes:
273: This function should be called after `SVDSolve()` has finished.
275: The value `nconv` may be different from the number of requested solutions
276: `nsv`, but not larger than `ncv`, see `SVDSetDimensions()`.
278: Level: beginner
280: .seealso: [](ch:svd), `SVDSetDimensions()`, `SVDSolve()`, `SVDGetSingularTriplet()`
281: @*/
282: PetscErrorCode SVDGetConverged(SVD svd,PetscInt *nconv)
283: {
284: PetscFunctionBegin;
286: PetscAssertPointer(nconv,2);
287: SVDCheckSolved(svd,1);
288: *nconv = svd->nconv;
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: /*@
293: SVDGetSingularTriplet - Gets the `i`-th triplet of the singular value decomposition
294: as computed by `SVDSolve()`. The solution consists in the singular value and its left
295: and right singular vectors.
297: Collective
299: Input Parameters:
300: + svd - the singular value solver context
301: - i - index of the solution
303: Output Parameters:
304: + sigma - singular value
305: . u - left singular vector
306: - v - right singular vector
308: Note:
309: Both `u` or `v` can be `NULL` if singular vectors are not required.
310: Otherwise, the caller must provide valid `Vec` objects, i.e.,
311: they must be created by the calling program with e.g. `MatCreateVecs()`.
313: The index `i` should be a value between 0 and `nconv`-1 (see `SVDGetConverged()`).
314: Singular triplets are indexed according to the ordering criterion established
315: with `SVDSetWhichSingularTriplets()`.
317: In the case of GSVD, the solution consists in three vectors $u$, $v$, $x$ that are
318: returned as follows. Vector $x$ is returned in the right singular vector
319: (argument `v`) and has length equal to the number of columns of $A$ and $B$.
320: The other two vectors are returned stacked on top of each other $[u^*,v^*]^*$ in
321: the left singular vector argument `u`, with length equal to $m+n$ (number of rows
322: of $A$ plus number of rows of $B$).
324: Level: beginner
326: .seealso: [](ch:svd), `SVDSolve()`, `SVDGetConverged()`, `SVDSetWhichSingularTriplets()`
327: @*/
328: PetscErrorCode SVDGetSingularTriplet(SVD svd,PetscInt i,PetscReal *sigma,Vec u,Vec v)
329: {
330: PetscInt M,N;
331: Vec w;
333: PetscFunctionBegin;
336: SVDCheckSolved(svd,1);
339: PetscCheck(i>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
340: PetscCheck(i<svd->nconv,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see SVDGetConverged()");
341: if (sigma) *sigma = svd->sigma[svd->perm[i]];
342: if (u || v) {
343: if (!svd->isgeneralized) {
344: PetscCall(MatGetSize(svd->OP,&M,&N));
345: if (M<N) { w = u; u = v; v = w; }
346: }
347: PetscCall(SVDComputeVectors(svd));
348: if (u) PetscCall(BVCopyVec(svd->U,svd->perm[i],u));
349: if (v) PetscCall(BVCopyVec(svd->V,svd->perm[i],v));
350: }
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: /*
355: SVDComputeResidualNorms_Standard - Computes the norms of the left and
356: right residuals associated with the i-th computed singular triplet.
358: Input Parameters:
359: sigma - singular value
360: u,v - singular vectors
361: x,y - two work vectors with the same dimensions as u,v
362: */
363: static PetscErrorCode SVDComputeResidualNorms_Standard(SVD svd,PetscReal sigma,Vec u,Vec v,Vec x,Vec y,PetscReal *norm1,PetscReal *norm2)
364: {
365: PetscInt M,N;
367: PetscFunctionBegin;
368: /* norm1 = ||A*v-sigma*u||_2 */
369: if (norm1) {
370: PetscCall(MatMult(svd->OP,v,x));
371: PetscCall(VecAXPY(x,-sigma,u));
372: PetscCall(VecNorm(x,NORM_2,norm1));
373: }
374: /* norm2 = ||A^T*u-sigma*v||_2 */
375: if (norm2) {
376: PetscCall(MatGetSize(svd->OP,&M,&N));
377: if (M<N) PetscCall(MatMult(svd->A,u,y));
378: else PetscCall(MatMult(svd->AT,u,y));
379: PetscCall(VecAXPY(y,-sigma,v));
380: PetscCall(VecNorm(y,NORM_2,norm2));
381: }
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: /*
386: SVDComputeResidualNorms_Generalized - In GSVD, compute the residual norms
387: norm1 = ||s^2*A'*u-c*B'*B*x||_2 and norm2 = ||c^2*B'*v-s*A'*A*x||_2.
389: Input Parameters:
390: sigma - singular value
391: uv - left singular vectors [u;v]
392: x - right singular vector
393: y,z - two work vectors with the same dimension as x
394: */
395: static PetscErrorCode SVDComputeResidualNorms_Generalized(SVD svd,PetscReal sigma,Vec uv,Vec x,Vec y,Vec z,PetscReal *norm1,PetscReal *norm2)
396: {
397: Vec u,v,ax,bx,nest,aux[2];
398: PetscReal c,s;
400: PetscFunctionBegin;
401: PetscCall(MatCreateVecs(svd->OP,NULL,&u));
402: PetscCall(MatCreateVecs(svd->OPb,NULL,&v));
403: aux[0] = u;
404: aux[1] = v;
405: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,aux,&nest));
406: PetscCall(VecCopy(uv,nest));
408: s = 1.0/PetscSqrtReal(1.0+sigma*sigma);
409: c = sigma*s;
411: /* norm1 = ||s^2*A'*u-c*B'*B*x||_2 */
412: if (norm1) {
413: PetscCall(VecDuplicate(v,&bx));
414: PetscCall(MatMultHermitianTranspose(svd->OP,u,z));
415: PetscCall(MatMult(svd->OPb,x,bx));
416: PetscCall(MatMultHermitianTranspose(svd->OPb,bx,y));
417: PetscCall(VecAXPBY(y,s*s,-c,z));
418: PetscCall(VecNorm(y,NORM_2,norm1));
419: PetscCall(VecDestroy(&bx));
420: }
421: /* norm2 = ||c^2*B'*v-s*A'*A*x||_2 */
422: if (norm2) {
423: PetscCall(VecDuplicate(u,&ax));
424: PetscCall(MatMultHermitianTranspose(svd->OPb,v,z));
425: PetscCall(MatMult(svd->OP,x,ax));
426: PetscCall(MatMultHermitianTranspose(svd->OP,ax,y));
427: PetscCall(VecAXPBY(y,c*c,-s,z));
428: PetscCall(VecNorm(y,NORM_2,norm2));
429: PetscCall(VecDestroy(&ax));
430: }
432: PetscCall(VecDestroy(&v));
433: PetscCall(VecDestroy(&u));
434: PetscCall(VecDestroy(&nest));
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: /*
439: SVDComputeResidualNorms_Hyperbolic - Computes the norms of the left and
440: right residuals associated with the i-th computed singular triplet.
442: Input Parameters:
443: sigma - singular value
444: sign - corresponding element of the signature Omega2
445: u,v - singular vectors
446: x,y,z - three work vectors with the same dimensions as u,v,u
447: */
448: static PetscErrorCode SVDComputeResidualNorms_Hyperbolic(SVD svd,PetscReal sigma,PetscReal sign,Vec u,Vec v,Vec x,Vec y,Vec z,PetscReal *norm1,PetscReal *norm2)
449: {
450: PetscInt M,N;
452: PetscFunctionBegin;
453: /* norm1 = ||A*v-sigma*u||_2 */
454: if (norm1) {
455: PetscCall(MatMult(svd->OP,v,x));
456: PetscCall(VecAXPY(x,-sigma,u));
457: PetscCall(VecNorm(x,NORM_2,norm1));
458: }
459: /* norm2 = ||A^T*Omega*u-sigma*sign*v||_2 */
460: if (norm2) {
461: PetscCall(MatGetSize(svd->OP,&M,&N));
462: PetscCall(VecPointwiseMult(z,u,svd->omega));
463: if (M<N) PetscCall(MatMult(svd->A,z,y));
464: else PetscCall(MatMult(svd->AT,z,y));
465: PetscCall(VecAXPY(y,-sigma*sign,v));
466: PetscCall(VecNorm(y,NORM_2,norm2));
467: }
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
471: /*@
472: SVDComputeError - Computes the error (based on the residual norm) associated
473: with the `i`-th singular triplet.
475: Collective
477: Input Parameters:
478: + svd - the singular value solver context
479: . i - the solution index
480: - type - the type of error to compute, see `SVDErrorType`
482: Output Parameter:
483: . error - the error
485: Notes:
486: The error can be computed in various ways, all of them based on the residual
487: norm obtained as $\sqrt{\eta_1^2+\eta_2^2}$ with $\eta_1 = \|Av-\sigma u\|_2$ and
488: $\eta_2 = \|A^*u-\sigma v\|_2$, where $(\sigma,u,v)$ is the approximate singular
489: triplet.
491: In the case of the GSVD, the two components of the residual norm are
492: $\eta_1 = \|s^2 A^*u-cB^*Bx\|_2$ and $\eta_2 = ||c^2 B^*v-sA^*Ax||_2$, where
493: $(\sigma,u,v,x)$ is the approximate generalized singular quadruple, with
494: $\sigma=c/s$.
496: Level: beginner
498: .seealso: [](ch:svd), `SVDErrorType`, `SVDSolve()`
499: @*/
500: PetscErrorCode SVDComputeError(SVD svd,PetscInt i,SVDErrorType type,PetscReal *error)
501: {
502: PetscReal sigma,norm1,norm2,c,s;
503: Vec u=NULL,v=NULL,x=NULL,y=NULL,z=NULL;
504: PetscReal vecnorm=1.0;
506: PetscFunctionBegin;
510: PetscAssertPointer(error,4);
511: SVDCheckSolved(svd,1);
513: /* allocate work vectors */
514: switch (svd->problem_type) {
515: case SVD_STANDARD:
516: PetscCall(SVDSetWorkVecs(svd,2,2));
517: u = svd->workl[0];
518: v = svd->workr[0];
519: x = svd->workl[1];
520: y = svd->workr[1];
521: break;
522: case SVD_GENERALIZED:
523: PetscCall(SVDSetWorkVecs(svd,1,3));
524: u = svd->workl[0];
525: v = svd->workr[0];
526: x = svd->workr[1];
527: y = svd->workr[2];
528: break;
529: case SVD_HYPERBOLIC:
530: PetscCall(SVDSetWorkVecs(svd,3,2));
531: u = svd->workl[0];
532: v = svd->workr[0];
533: x = svd->workl[1];
534: y = svd->workr[1];
535: z = svd->workl[2];
536: break;
537: }
539: /* compute residual norm */
540: PetscCall(SVDGetSingularTriplet(svd,i,&sigma,u,v));
541: switch (svd->problem_type) {
542: case SVD_STANDARD:
543: PetscCall(SVDComputeResidualNorms_Standard(svd,sigma,u,v,x,y,&norm1,&norm2));
544: break;
545: case SVD_GENERALIZED:
546: PetscCall(SVDComputeResidualNorms_Generalized(svd,sigma,u,v,x,y,&norm1,&norm2));
547: break;
548: case SVD_HYPERBOLIC:
549: PetscCall(SVDComputeResidualNorms_Hyperbolic(svd,sigma,svd->sign[svd->perm[i]],u,v,x,y,z,&norm1,&norm2));
550: break;
551: }
552: *error = SlepcAbs(norm1,norm2);
554: /* compute 2-norm of eigenvector of the cyclic form */
555: if (type!=SVD_ERROR_ABSOLUTE) {
556: switch (svd->problem_type) {
557: case SVD_STANDARD:
558: vecnorm = PETSC_SQRT2;
559: break;
560: case SVD_GENERALIZED:
561: PetscCall(VecNorm(v,NORM_2,&vecnorm));
562: vecnorm = PetscSqrtReal(1.0+vecnorm*vecnorm);
563: break;
564: case SVD_HYPERBOLIC:
565: PetscCall(VecNorm(u,NORM_2,&vecnorm));
566: vecnorm = PetscSqrtReal(1.0+vecnorm*vecnorm);
567: break;
568: }
569: }
571: /* compute error */
572: switch (type) {
573: case SVD_ERROR_ABSOLUTE:
574: break;
575: case SVD_ERROR_RELATIVE:
576: if (svd->isgeneralized) {
577: s = 1.0/PetscSqrtReal(1.0+sigma*sigma);
578: c = sigma*s;
579: norm1 /= c*vecnorm;
580: norm2 /= s*vecnorm;
581: *error = PetscMax(norm1,norm2);
582: } else *error /= sigma*vecnorm;
583: break;
584: case SVD_ERROR_NORM:
585: if (!svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
586: if (svd->isgeneralized && !svd->nrmb) PetscCall(MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb));
587: *error /= PetscMax(svd->nrma,svd->nrmb)*vecnorm;
588: break;
589: default:
590: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
591: }
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }