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: }