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 - singular value solver context obtained from SVDCreate()

 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, 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:    All the command-line options listed above admit an optional argument specifying
107:    the viewer type and options. For instance, use '-svd_view_mat0 binary:amatrix.bin'
108:    to save the A matrix to a binary file, '-svd_view_values draw' to draw the computed
109:    singular values graphically, or '-svd_error_relative :myerr.m:ascii_matlab' to save
110:    the errors in a file that can be executed in Matlab.

112:    Level: beginner

114: .seealso: `SVDCreate()`, `SVDSetUp()`, `SVDDestroy()`
115: @*/
116: PetscErrorCode SVDSolve(SVD svd)
117: {
118:   PetscInt       i,m,n,*workperm;

120:   PetscFunctionBegin;
122:   if (svd->state>=SVD_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
123:   PetscCall(PetscLogEventBegin(SVD_Solve,svd,0,0,0));

125:   /* call setup */
126:   PetscCall(SVDSetUp(svd));

128:   /* safeguard for matrices with zero rows or columns */
129:   PetscCall(MatGetSize(svd->OP,&m,&n));
130:   if (m == 0 || n == 0) {
131:     svd->nconv  = 0;
132:     svd->reason = SVD_CONVERGED_TOL;
133:     svd->state  = SVD_STATE_SOLVED;
134:     PetscFunctionReturn(PETSC_SUCCESS);
135:   }

137:   svd->its = 0;
138:   svd->nconv = 0;
139:   for (i=0;i<svd->ncv;i++) {
140:     svd->sigma[i]  = 0.0;
141:     svd->errest[i] = 0.0;
142:     svd->perm[i]   = i;
143:   }
144:   PetscCall(SVDViewFromOptions(svd,NULL,"-svd_view_pre"));

146:   switch (svd->problem_type) {
147:     case SVD_STANDARD:
148:       PetscUseTypeMethod(svd,solve);
149:       break;
150:     case SVD_GENERALIZED:
151:       PetscUseTypeMethod(svd,solveg);
152:       break;
153:     case SVD_HYPERBOLIC:
154:       PetscUseTypeMethod(svd,solveh);
155:       break;
156:   }
157:   svd->state = SVD_STATE_SOLVED;

159:   /* sort singular triplets */
160:   if (svd->which == SVD_SMALLEST) PetscCall(PetscSortRealWithPermutation(svd->nconv,svd->sigma,svd->perm));
161:   else {
162:     PetscCall(PetscMalloc1(svd->nconv,&workperm));
163:     for (i=0;i<svd->nconv;i++) workperm[i] = i;
164:     PetscCall(PetscSortRealWithPermutation(svd->nconv,svd->sigma,workperm));
165:     for (i=0;i<svd->nconv;i++) svd->perm[i] = workperm[svd->nconv-i-1];
166:     PetscCall(PetscFree(workperm));
167:   }
168:   PetscCall(PetscLogEventEnd(SVD_Solve,svd,0,0,0));

170:   /* various viewers */
171:   PetscCall(SVDViewFromOptions(svd,NULL,"-svd_view"));
172:   PetscCall(SVDConvergedReasonViewFromOptions(svd));
173:   PetscCall(SVDErrorViewFromOptions(svd));
174:   PetscCall(SVDValuesViewFromOptions(svd));
175:   PetscCall(SVDVectorsViewFromOptions(svd));
176:   PetscCall(MatViewFromOptions(svd->OP,(PetscObject)svd,"-svd_view_mat0"));
177:   if (svd->isgeneralized) PetscCall(MatViewFromOptions(svd->OPb,(PetscObject)svd,"-svd_view_mat1"));
178:   if (svd->ishyperbolic) PetscCall(VecViewFromOptions(svd->omega,(PetscObject)svd,"-svd_view_signature"));

180:   /* Remove the initial subspaces */
181:   svd->nini = 0;
182:   svd->ninil = 0;
183:   PetscFunctionReturn(PETSC_SUCCESS);
184: }

186: /*@
187:    SVDGetIterationNumber - Gets the current iteration number. If the
188:    call to SVDSolve() is complete, then it returns the number of iterations
189:    carried out by the solution method.

191:    Not Collective

193:    Input Parameter:
194: .  svd - the singular value solver context

196:    Output Parameter:
197: .  its - number of iterations

199:    Note:
200:    During the i-th iteration this call returns i-1. If SVDSolve() is
201:    complete, then parameter "its" contains either the iteration number at
202:    which convergence was successfully reached, or failure was detected.
203:    Call SVDGetConvergedReason() to determine if the solver converged or
204:    failed and why.

206:    Level: intermediate

208: .seealso: `SVDGetConvergedReason()`, `SVDSetTolerances()`
209: @*/
210: PetscErrorCode SVDGetIterationNumber(SVD svd,PetscInt *its)
211: {
212:   PetscFunctionBegin;
214:   PetscAssertPointer(its,2);
215:   *its = svd->its;
216:   PetscFunctionReturn(PETSC_SUCCESS);
217: }

219: /*@
220:    SVDGetConvergedReason - Gets the reason why the SVDSolve() iteration was
221:    stopped.

223:    Not Collective

225:    Input Parameter:
226: .  svd - the singular value solver context

228:    Output Parameter:
229: .  reason - negative value indicates diverged, positive value converged
230:    (see SVDConvergedReason)

232:    Options Database Key:
233: .  -svd_converged_reason - print the reason to a viewer

235:    Notes:
236:    Possible values for reason are
237: +  SVD_CONVERGED_TOL - converged up to tolerance
238: .  SVD_CONVERGED_USER - converged due to a user-defined condition
239: .  SVD_CONVERGED_MAXIT - reached the maximum number of iterations with SVD_CONV_MAXIT criterion
240: .  SVD_DIVERGED_ITS - required more than max_it iterations to reach convergence
241: .  SVD_DIVERGED_BREAKDOWN - generic breakdown in method
242: -  SVD_DIVERGED_SYMMETRY_LOST - underlying indefinite eigensolver was not able to keep symmetry

244:    Can only be called after the call to SVDSolve() is complete.

246:    Level: intermediate

248: .seealso: `SVDSetTolerances()`, `SVDSolve()`, `SVDConvergedReason`
249: @*/
250: PetscErrorCode SVDGetConvergedReason(SVD svd,SVDConvergedReason *reason)
251: {
252:   PetscFunctionBegin;
254:   PetscAssertPointer(reason,2);
255:   SVDCheckSolved(svd,1);
256:   *reason = svd->reason;
257:   PetscFunctionReturn(PETSC_SUCCESS);
258: }

260: /*@
261:    SVDGetConverged - Gets the number of converged singular values.

263:    Not Collective

265:    Input Parameter:
266: .  svd - the singular value solver context

268:    Output Parameter:
269: .  nconv - number of converged singular values

271:    Note:
272:    This function should be called after SVDSolve() has finished.

274:    Level: beginner

276: .seealso: `SVDSetDimensions()`, `SVDSolve()`, `SVDGetSingularTriplet()`
277: @*/
278: PetscErrorCode SVDGetConverged(SVD svd,PetscInt *nconv)
279: {
280:   PetscFunctionBegin;
282:   PetscAssertPointer(nconv,2);
283:   SVDCheckSolved(svd,1);
284:   *nconv = svd->nconv;
285:   PetscFunctionReturn(PETSC_SUCCESS);
286: }

288: /*@
289:    SVDGetSingularTriplet - Gets the i-th triplet of the singular value decomposition
290:    as computed by SVDSolve(). The solution consists in the singular value and its left
291:    and right singular vectors.

293:    Collective

295:    Input Parameters:
296: +  svd - singular value solver context
297: -  i   - index of the solution

299:    Output Parameters:
300: +  sigma - singular value
301: .  u     - left singular vector
302: -  v     - right singular vector

304:    Note:
305:    Both u or v can be NULL if singular vectors are not required.
306:    Otherwise, the caller must provide valid Vec objects, i.e.,
307:    they must be created by the calling program with e.g. MatCreateVecs().

309:    The index i should be a value between 0 and nconv-1 (see SVDGetConverged()).
310:    Singular triplets are indexed according to the ordering criterion established
311:    with SVDSetWhichSingularTriplets().

313:    In the case of GSVD, the solution consists in three vectors u,v,x that are
314:    returned as follows. Vector x is returned in the right singular vector
315:    (argument v) and has length equal to the number of columns of A and B.
316:    The other two vectors are returned stacked on top of each other [u;v] in
317:    the left singular vector argument, with length equal to m+n (number of rows
318:    of A plus number of rows of B).

320:    Level: beginner

322: .seealso: `SVDSolve()`, `SVDGetConverged()`, `SVDSetWhichSingularTriplets()`
323: @*/
324: PetscErrorCode SVDGetSingularTriplet(SVD svd,PetscInt i,PetscReal *sigma,Vec u,Vec v)
325: {
326:   PetscInt       M,N;
327:   Vec            w;

329:   PetscFunctionBegin;
332:   SVDCheckSolved(svd,1);
335:   PetscCheck(i>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
336:   PetscCheck(i<svd->nconv,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see SVDGetConverged()");
337:   if (sigma) *sigma = svd->sigma[svd->perm[i]];
338:   if (u || v) {
339:     if (!svd->isgeneralized) {
340:       PetscCall(MatGetSize(svd->OP,&M,&N));
341:       if (M<N) { w = u; u = v; v = w; }
342:     }
343:     PetscCall(SVDComputeVectors(svd));
344:     if (u) PetscCall(BVCopyVec(svd->U,svd->perm[i],u));
345:     if (v) PetscCall(BVCopyVec(svd->V,svd->perm[i],v));
346:   }
347:   PetscFunctionReturn(PETSC_SUCCESS);
348: }

350: /*
351:    SVDComputeResidualNorms_Standard - Computes the norms of the left and
352:    right residuals associated with the i-th computed singular triplet.

354:    Input Parameters:
355:      sigma - singular value
356:      u,v   - singular vectors
357:      x,y   - two work vectors with the same dimensions as u,v
358: */
359: static PetscErrorCode SVDComputeResidualNorms_Standard(SVD svd,PetscReal sigma,Vec u,Vec v,Vec x,Vec y,PetscReal *norm1,PetscReal *norm2)
360: {
361:   PetscInt       M,N;

363:   PetscFunctionBegin;
364:   /* norm1 = ||A*v-sigma*u||_2 */
365:   if (norm1) {
366:     PetscCall(MatMult(svd->OP,v,x));
367:     PetscCall(VecAXPY(x,-sigma,u));
368:     PetscCall(VecNorm(x,NORM_2,norm1));
369:   }
370:   /* norm2 = ||A^T*u-sigma*v||_2 */
371:   if (norm2) {
372:     PetscCall(MatGetSize(svd->OP,&M,&N));
373:     if (M<N) PetscCall(MatMult(svd->A,u,y));
374:     else PetscCall(MatMult(svd->AT,u,y));
375:     PetscCall(VecAXPY(y,-sigma,v));
376:     PetscCall(VecNorm(y,NORM_2,norm2));
377:   }
378:   PetscFunctionReturn(PETSC_SUCCESS);
379: }

381: /*
382:    SVDComputeResidualNorms_Generalized - In GSVD, compute the residual norms
383:    norm1 = ||s^2*A'*u-c*B'*B*x||_2 and norm2 = ||c^2*B'*v-s*A'*A*x||_2.

385:    Input Parameters:
386:      sigma - singular value
387:      uv    - left singular vectors [u;v]
388:      x     - right singular vector
389:      y,z   - two work vectors with the same dimension as x
390: */
391: static PetscErrorCode SVDComputeResidualNorms_Generalized(SVD svd,PetscReal sigma,Vec uv,Vec x,Vec y,Vec z,PetscReal *norm1,PetscReal *norm2)
392: {
393:   Vec            u,v,ax,bx,nest,aux[2];
394:   PetscReal      c,s;

396:   PetscFunctionBegin;
397:   PetscCall(MatCreateVecs(svd->OP,NULL,&u));
398:   PetscCall(MatCreateVecs(svd->OPb,NULL,&v));
399:   aux[0] = u;
400:   aux[1] = v;
401:   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,aux,&nest));
402:   PetscCall(VecCopy(uv,nest));

404:   s = 1.0/PetscSqrtReal(1.0+sigma*sigma);
405:   c = sigma*s;

407:   /* norm1 = ||s^2*A'*u-c*B'*B*x||_2 */
408:   if (norm1) {
409:     PetscCall(VecDuplicate(v,&bx));
410:     PetscCall(MatMultHermitianTranspose(svd->OP,u,z));
411:     PetscCall(MatMult(svd->OPb,x,bx));
412:     PetscCall(MatMultHermitianTranspose(svd->OPb,bx,y));
413:     PetscCall(VecAXPBY(y,s*s,-c,z));
414:     PetscCall(VecNorm(y,NORM_2,norm1));
415:     PetscCall(VecDestroy(&bx));
416:   }
417:   /* norm2 = ||c^2*B'*v-s*A'*A*x||_2 */
418:   if (norm2) {
419:     PetscCall(VecDuplicate(u,&ax));
420:     PetscCall(MatMultHermitianTranspose(svd->OPb,v,z));
421:     PetscCall(MatMult(svd->OP,x,ax));
422:     PetscCall(MatMultHermitianTranspose(svd->OP,ax,y));
423:     PetscCall(VecAXPBY(y,c*c,-s,z));
424:     PetscCall(VecNorm(y,NORM_2,norm2));
425:     PetscCall(VecDestroy(&ax));
426:   }

428:   PetscCall(VecDestroy(&v));
429:   PetscCall(VecDestroy(&u));
430:   PetscCall(VecDestroy(&nest));
431:   PetscFunctionReturn(PETSC_SUCCESS);
432: }

434: /*
435:    SVDComputeResidualNorms_Hyperbolic - Computes the norms of the left and
436:    right residuals associated with the i-th computed singular triplet.

438:    Input Parameters:
439:      sigma - singular value
440:      sign  - corresponding element of the signature Omega2
441:      u,v   - singular vectors
442:      x,y,z - three work vectors with the same dimensions as u,v,u
443: */
444: static PetscErrorCode SVDComputeResidualNorms_Hyperbolic(SVD svd,PetscReal sigma,PetscReal sign,Vec u,Vec v,Vec x,Vec y,Vec z,PetscReal *norm1,PetscReal *norm2)
445: {
446:   PetscInt       M,N;

448:   PetscFunctionBegin;
449:   /* norm1 = ||A*v-sigma*u||_2 */
450:   if (norm1) {
451:     PetscCall(MatMult(svd->OP,v,x));
452:     PetscCall(VecAXPY(x,-sigma,u));
453:     PetscCall(VecNorm(x,NORM_2,norm1));
454:   }
455:   /* norm2 = ||A^T*Omega*u-sigma*sign*v||_2 */
456:   if (norm2) {
457:     PetscCall(MatGetSize(svd->OP,&M,&N));
458:     PetscCall(VecPointwiseMult(z,u,svd->omega));
459:     if (M<N) PetscCall(MatMult(svd->A,z,y));
460:     else PetscCall(MatMult(svd->AT,z,y));
461:     PetscCall(VecAXPY(y,-sigma*sign,v));
462:     PetscCall(VecNorm(y,NORM_2,norm2));
463:   }
464:   PetscFunctionReturn(PETSC_SUCCESS);
465: }

467: /*@
468:    SVDComputeError - Computes the error (based on the residual norm) associated
469:    with the i-th singular triplet.

471:    Collective

473:    Input Parameters:
474: +  svd  - the singular value solver context
475: .  i    - the solution index
476: -  type - the type of error to compute

478:    Output Parameter:
479: .  error - the error

481:    Notes:
482:    The error can be computed in various ways, all of them based on the residual
483:    norm obtained as sqrt(n1^2+n2^2) with n1 = ||A*v-sigma*u||_2 and
484:    n2 = ||A^T*u-sigma*v||_2, where sigma is the singular value, u is the left
485:    singular vector and v is the right singular vector.

487:    In the case of the GSVD, the two components of the residual norm are
488:    n1 = ||s^2*A'*u-c*B'*B*x||_2 and n2 = ||c^2*B'*v-s*A'*A*x||_2, where [u;v]
489:    are the left singular vectors and x is the right singular vector, with
490:    sigma=c/s.

492:    Level: beginner

494: .seealso: `SVDErrorType`, `SVDSolve()`
495: @*/
496: PetscErrorCode SVDComputeError(SVD svd,PetscInt i,SVDErrorType type,PetscReal *error)
497: {
498:   PetscReal      sigma,norm1,norm2,c,s;
499:   Vec            u=NULL,v=NULL,x=NULL,y=NULL,z=NULL;
500:   PetscReal      vecnorm=1.0;

502:   PetscFunctionBegin;
506:   PetscAssertPointer(error,4);
507:   SVDCheckSolved(svd,1);

509:   /* allocate work vectors */
510:   switch (svd->problem_type) {
511:     case SVD_STANDARD:
512:       PetscCall(SVDSetWorkVecs(svd,2,2));
513:       u = svd->workl[0];
514:       v = svd->workr[0];
515:       x = svd->workl[1];
516:       y = svd->workr[1];
517:       break;
518:     case SVD_GENERALIZED:
519:       PetscCall(SVDSetWorkVecs(svd,1,3));
520:       u = svd->workl[0];
521:       v = svd->workr[0];
522:       x = svd->workr[1];
523:       y = svd->workr[2];
524:       break;
525:     case SVD_HYPERBOLIC:
526:       PetscCall(SVDSetWorkVecs(svd,3,2));
527:       u = svd->workl[0];
528:       v = svd->workr[0];
529:       x = svd->workl[1];
530:       y = svd->workr[1];
531:       z = svd->workl[2];
532:       break;
533:   }

535:   /* compute residual norm */
536:   PetscCall(SVDGetSingularTriplet(svd,i,&sigma,u,v));
537:   switch (svd->problem_type) {
538:     case SVD_STANDARD:
539:       PetscCall(SVDComputeResidualNorms_Standard(svd,sigma,u,v,x,y,&norm1,&norm2));
540:       break;
541:     case SVD_GENERALIZED:
542:       PetscCall(SVDComputeResidualNorms_Generalized(svd,sigma,u,v,x,y,&norm1,&norm2));
543:       break;
544:     case SVD_HYPERBOLIC:
545:       PetscCall(SVDComputeResidualNorms_Hyperbolic(svd,sigma,svd->sign[svd->perm[i]],u,v,x,y,z,&norm1,&norm2));
546:       break;
547:   }
548:   *error = SlepcAbs(norm1,norm2);

550:   /* compute 2-norm of eigenvector of the cyclic form */
551:   if (type!=SVD_ERROR_ABSOLUTE) {
552:     switch (svd->problem_type) {
553:       case SVD_STANDARD:
554:         vecnorm = PETSC_SQRT2;
555:         break;
556:       case SVD_GENERALIZED:
557:         PetscCall(VecNorm(v,NORM_2,&vecnorm));
558:         vecnorm = PetscSqrtReal(1.0+vecnorm*vecnorm);
559:         break;
560:       case SVD_HYPERBOLIC:
561:         PetscCall(VecNorm(u,NORM_2,&vecnorm));
562:         vecnorm = PetscSqrtReal(1.0+vecnorm*vecnorm);
563:         break;
564:     }
565:   }

567:   /* compute error */
568:   switch (type) {
569:     case SVD_ERROR_ABSOLUTE:
570:       break;
571:     case SVD_ERROR_RELATIVE:
572:       if (svd->isgeneralized) {
573:         s = 1.0/PetscSqrtReal(1.0+sigma*sigma);
574:         c = sigma*s;
575:         norm1 /= c*vecnorm;
576:         norm2 /= s*vecnorm;
577:         *error = PetscMax(norm1,norm2);
578:       } else *error /= sigma*vecnorm;
579:       break;
580:     case SVD_ERROR_NORM:
581:       if (!svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
582:       if (svd->isgeneralized && !svd->nrmb) PetscCall(MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb));
583:       *error /= PetscMax(svd->nrma,svd->nrmb)*vecnorm;
584:       break;
585:     default:
586:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
587:   }
588:   PetscFunctionReturn(PETSC_SUCCESS);
589: }