LCOV - code coverage report
Current view: top level - svd/interface - svdsolve.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 255 262 97.3 %
Date: 2024-03-29 00:34:52 Functions: 11 11 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       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             : */
      13             : 
      14             : #include <slepc/private/svdimpl.h>   /*I "slepcsvd.h" I*/
      15             : #include <slepc/private/bvimpl.h>
      16             : 
      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          48 : PetscErrorCode SVDComputeVectors_Left(SVD svd)
      22             : {
      23          48 :   Vec                tl,omega2,u,v,w;
      24          48 :   PetscInt           i,oldsize;
      25          48 :   VecType            vtype;
      26          48 :   const PetscScalar* varray;
      27             : 
      28          48 :   PetscFunctionBegin;
      29          48 :   if (!svd->leftbasis) {
      30             :     /* generate left singular vectors on U */
      31          34 :     if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
      32          34 :     PetscCall(BVGetSizes(svd->U,NULL,NULL,&oldsize));
      33          34 :     if (!oldsize) {
      34          29 :       if (!((PetscObject)svd->U)->type_name) PetscCall(BVSetType(svd->U,((PetscObject)svd->V)->type_name));
      35          29 :       PetscCall(MatCreateVecsEmpty(svd->A,NULL,&tl));
      36          29 :       PetscCall(BVSetSizesFromVec(svd->U,tl,svd->ncv));
      37          29 :       PetscCall(VecDestroy(&tl));
      38             :     }
      39          34 :     PetscCall(BVSetActiveColumns(svd->V,0,svd->nconv));
      40          34 :     PetscCall(BVSetActiveColumns(svd->U,0,svd->nconv));
      41          34 :     if (!svd->ishyperbolic) PetscCall(BVMatMult(svd->V,svd->A,svd->U));
      42           7 :     else if (svd->swapped) {  /* compute right singular vectors as V=A'*Omega*U */
      43           2 :       PetscCall(MatCreateVecs(svd->A,&w,NULL));
      44          58 :       for (i=0;i<svd->nconv;i++) {
      45          56 :         PetscCall(BVGetColumn(svd->V,i,&v));
      46          56 :         PetscCall(BVGetColumn(svd->U,i,&u));
      47          56 :         PetscCall(VecPointwiseMult(w,v,svd->omega));
      48          56 :         PetscCall(MatMult(svd->A,w,u));
      49          56 :         PetscCall(BVRestoreColumn(svd->V,i,&v));
      50          56 :         PetscCall(BVRestoreColumn(svd->U,i,&u));
      51             :       }
      52           2 :       PetscCall(VecDestroy(&w));
      53             :     } else {  /* compute left singular vectors as usual U=A*V, and set-up Omega-orthogonalization of U */
      54           5 :       PetscCall(BV_SetMatrixDiagonal(svd->U,svd->omega,svd->A));
      55           5 :       PetscCall(BVMatMult(svd->V,svd->A,svd->U));
      56             :     }
      57          34 :     PetscCall(BVOrthogonalize(svd->U,NULL));
      58          34 :     if (svd->ishyperbolic && !svd->swapped) {  /* store signature after Omega-orthogonalization */
      59           5 :       PetscCall(MatGetVecType(svd->A,&vtype));
      60           5 :       PetscCall(VecCreate(PETSC_COMM_SELF,&omega2));
      61           5 :       PetscCall(VecSetSizes(omega2,svd->nconv,svd->nconv));
      62           5 :       PetscCall(VecSetType(omega2,vtype));
      63           5 :       PetscCall(BVGetSignature(svd->U,omega2));
      64           5 :       PetscCall(VecGetArrayRead(omega2,&varray));
      65         187 :       for (i=0;i<svd->nconv;i++) {
      66         182 :         svd->sign[i] = PetscRealPart(varray[i]);
      67         182 :         if (PetscRealPart(varray[i])<0.0) PetscCall(BVScaleColumn(svd->U,i,-1.0));
      68             :       }
      69           5 :       PetscCall(VecRestoreArrayRead(omega2,&varray));
      70           5 :       PetscCall(VecDestroy(&omega2));
      71             :     }
      72             :   }
      73          48 :   PetscFunctionReturn(PETSC_SUCCESS);
      74             : }
      75             : 
      76        2286 : PetscErrorCode SVDComputeVectors(SVD svd)
      77             : {
      78        2286 :   PetscFunctionBegin;
      79        2286 :   SVDCheckSolved(svd,1);
      80        2286 :   if (svd->state==SVD_STATE_SOLVED) PetscTryTypeMethod(svd,computevectors);
      81        2286 :   svd->state = SVD_STATE_VECTORS;
      82        2286 :   PetscFunctionReturn(PETSC_SUCCESS);
      83             : }
      84             : 
      85             : /*@
      86             :    SVDSolve - Solves the singular value problem.
      87             : 
      88             :    Collective
      89             : 
      90             :    Input Parameter:
      91             : .  svd - singular value solver context obtained from SVDCreate()
      92             : 
      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
     104             : 
     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.
     111             : 
     112             :    Level: beginner
     113             : 
     114             : .seealso: SVDCreate(), SVDSetUp(), SVDDestroy()
     115             : @*/
     116         247 : PetscErrorCode SVDSolve(SVD svd)
     117             : {
     118         247 :   PetscInt       i,*workperm;
     119             : 
     120         247 :   PetscFunctionBegin;
     121         247 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     122         247 :   if (svd->state>=SVD_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
     123         247 :   PetscCall(PetscLogEventBegin(SVD_Solve,svd,0,0,0));
     124             : 
     125             :   /* call setup */
     126         247 :   PetscCall(SVDSetUp(svd));
     127         247 :   svd->its = 0;
     128         247 :   svd->nconv = 0;
     129        5044 :   for (i=0;i<svd->ncv;i++) {
     130        4797 :     svd->sigma[i]  = 0.0;
     131        4797 :     svd->errest[i] = 0.0;
     132        4797 :     svd->perm[i]   = i;
     133             :   }
     134         247 :   PetscCall(SVDViewFromOptions(svd,NULL,"-svd_view_pre"));
     135             : 
     136         247 :   switch (svd->problem_type) {
     137         136 :     case SVD_STANDARD:
     138         136 :       PetscUseTypeMethod(svd,solve);
     139             :       break;
     140          85 :     case SVD_GENERALIZED:
     141          85 :       PetscUseTypeMethod(svd,solveg);
     142             :       break;
     143          26 :     case SVD_HYPERBOLIC:
     144          26 :       PetscUseTypeMethod(svd,solveh);
     145             :       break;
     146             :   }
     147         247 :   svd->state = SVD_STATE_SOLVED;
     148             : 
     149             :   /* sort singular triplets */
     150         247 :   if (svd->which == SVD_SMALLEST) PetscCall(PetscSortRealWithPermutation(svd->nconv,svd->sigma,svd->perm));
     151             :   else {
     152         222 :     PetscCall(PetscMalloc1(svd->nconv,&workperm));
     153        2675 :     for (i=0;i<svd->nconv;i++) workperm[i] = i;
     154         222 :     PetscCall(PetscSortRealWithPermutation(svd->nconv,svd->sigma,workperm));
     155        2675 :     for (i=0;i<svd->nconv;i++) svd->perm[i] = workperm[svd->nconv-i-1];
     156         222 :     PetscCall(PetscFree(workperm));
     157             :   }
     158         247 :   PetscCall(PetscLogEventEnd(SVD_Solve,svd,0,0,0));
     159             : 
     160             :   /* various viewers */
     161         247 :   PetscCall(SVDViewFromOptions(svd,NULL,"-svd_view"));
     162         247 :   PetscCall(SVDConvergedReasonViewFromOptions(svd));
     163         247 :   PetscCall(SVDErrorViewFromOptions(svd));
     164         247 :   PetscCall(SVDValuesViewFromOptions(svd));
     165         247 :   PetscCall(SVDVectorsViewFromOptions(svd));
     166         247 :   PetscCall(MatViewFromOptions(svd->OP,(PetscObject)svd,"-svd_view_mat0"));
     167         247 :   if (svd->isgeneralized) PetscCall(MatViewFromOptions(svd->OPb,(PetscObject)svd,"-svd_view_mat1"));
     168         247 :   if (svd->ishyperbolic) PetscCall(VecViewFromOptions(svd->omega,(PetscObject)svd,"-svd_view_signature"));
     169             : 
     170             :   /* Remove the initial subspaces */
     171         247 :   svd->nini = 0;
     172         247 :   svd->ninil = 0;
     173         247 :   PetscFunctionReturn(PETSC_SUCCESS);
     174             : }
     175             : 
     176             : /*@
     177             :    SVDGetIterationNumber - Gets the current iteration number. If the
     178             :    call to SVDSolve() is complete, then it returns the number of iterations
     179             :    carried out by the solution method.
     180             : 
     181             :    Not Collective
     182             : 
     183             :    Input Parameter:
     184             : .  svd - the singular value solver context
     185             : 
     186             :    Output Parameter:
     187             : .  its - number of iterations
     188             : 
     189             :    Note:
     190             :    During the i-th iteration this call returns i-1. If SVDSolve() is
     191             :    complete, then parameter "its" contains either the iteration number at
     192             :    which convergence was successfully reached, or failure was detected.
     193             :    Call SVDGetConvergedReason() to determine if the solver converged or
     194             :    failed and why.
     195             : 
     196             :    Level: intermediate
     197             : 
     198             : .seealso: SVDGetConvergedReason(), SVDSetTolerances()
     199             : @*/
     200          55 : PetscErrorCode SVDGetIterationNumber(SVD svd,PetscInt *its)
     201             : {
     202          55 :   PetscFunctionBegin;
     203          55 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     204          55 :   PetscAssertPointer(its,2);
     205          55 :   *its = svd->its;
     206          55 :   PetscFunctionReturn(PETSC_SUCCESS);
     207             : }
     208             : 
     209             : /*@
     210             :    SVDGetConvergedReason - Gets the reason why the SVDSolve() iteration was
     211             :    stopped.
     212             : 
     213             :    Not Collective
     214             : 
     215             :    Input Parameter:
     216             : .  svd - the singular value solver context
     217             : 
     218             :    Output Parameter:
     219             : .  reason - negative value indicates diverged, positive value converged
     220             :    (see SVDConvergedReason)
     221             : 
     222             :    Options Database Key:
     223             : .  -svd_converged_reason - print the reason to a viewer
     224             : 
     225             :    Notes:
     226             :    Possible values for reason are
     227             : +  SVD_CONVERGED_TOL - converged up to tolerance
     228             : .  SVD_CONVERGED_USER - converged due to a user-defined condition
     229             : .  SVD_CONVERGED_MAXIT - reached the maximum number of iterations with SVD_CONV_MAXIT criterion
     230             : .  SVD_DIVERGED_ITS - required more than max_it iterations to reach convergence
     231             : .  SVD_DIVERGED_BREAKDOWN - generic breakdown in method
     232             : -  SVD_DIVERGED_SYMMETRY_LOST - underlying indefinite eigensolver was not able to keep symmetry
     233             : 
     234             :    Can only be called after the call to SVDSolve() is complete.
     235             : 
     236             :    Level: intermediate
     237             : 
     238             : .seealso: SVDSetTolerances(), SVDSolve(), SVDConvergedReason
     239             : @*/
     240          12 : PetscErrorCode SVDGetConvergedReason(SVD svd,SVDConvergedReason *reason)
     241             : {
     242          12 :   PetscFunctionBegin;
     243          12 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     244          12 :   PetscAssertPointer(reason,2);
     245          12 :   SVDCheckSolved(svd,1);
     246          12 :   *reason = svd->reason;
     247          12 :   PetscFunctionReturn(PETSC_SUCCESS);
     248             : }
     249             : 
     250             : /*@
     251             :    SVDGetConverged - Gets the number of converged singular values.
     252             : 
     253             :    Not Collective
     254             : 
     255             :    Input Parameter:
     256             : .  svd - the singular value solver context
     257             : 
     258             :    Output Parameter:
     259             : .  nconv - number of converged singular values
     260             : 
     261             :    Note:
     262             :    This function should be called after SVDSolve() has finished.
     263             : 
     264             :    Level: beginner
     265             : 
     266             : .seealso: SVDSetDimensions(), SVDSolve(), SVDGetSingularTriplet()
     267             : @*/
     268         115 : PetscErrorCode SVDGetConverged(SVD svd,PetscInt *nconv)
     269             : {
     270         115 :   PetscFunctionBegin;
     271         115 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     272         115 :   PetscAssertPointer(nconv,2);
     273         115 :   SVDCheckSolved(svd,1);
     274         115 :   *nconv = svd->nconv;
     275         115 :   PetscFunctionReturn(PETSC_SUCCESS);
     276             : }
     277             : 
     278             : /*@C
     279             :    SVDGetSingularTriplet - Gets the i-th triplet of the singular value decomposition
     280             :    as computed by SVDSolve(). The solution consists in the singular value and its left
     281             :    and right singular vectors.
     282             : 
     283             :    Collective
     284             : 
     285             :    Input Parameters:
     286             : +  svd - singular value solver context
     287             : -  i   - index of the solution
     288             : 
     289             :    Output Parameters:
     290             : +  sigma - singular value
     291             : .  u     - left singular vector
     292             : -  v     - right singular vector
     293             : 
     294             :    Note:
     295             :    Both u or v can be NULL if singular vectors are not required.
     296             :    Otherwise, the caller must provide valid Vec objects, i.e.,
     297             :    they must be created by the calling program with e.g. MatCreateVecs().
     298             : 
     299             :    The index i should be a value between 0 and nconv-1 (see SVDGetConverged()).
     300             :    Singular triplets are indexed according to the ordering criterion established
     301             :    with SVDSetWhichSingularTriplets().
     302             : 
     303             :    In the case of GSVD, the solution consists in three vectors u,v,x that are
     304             :    returned as follows. Vector x is returned in the right singular vector
     305             :    (argument v) and has length equal to the number of columns of A and B.
     306             :    The other two vectors are returned stacked on top of each other [u;v] in
     307             :    the left singular vector argument, with length equal to m+n (number of rows
     308             :    of A plus number of rows of B).
     309             : 
     310             :    Level: beginner
     311             : 
     312             : .seealso: SVDSolve(), SVDGetConverged(), SVDSetWhichSingularTriplets()
     313             : @*/
     314        3385 : PetscErrorCode SVDGetSingularTriplet(SVD svd,PetscInt i,PetscReal *sigma,Vec u,Vec v)
     315             : {
     316        3385 :   PetscInt       M,N;
     317        3385 :   Vec            w;
     318             : 
     319        3385 :   PetscFunctionBegin;
     320        3385 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     321       13540 :   PetscValidLogicalCollectiveInt(svd,i,2);
     322        3385 :   SVDCheckSolved(svd,1);
     323        3385 :   if (u) { PetscValidHeaderSpecific(u,VEC_CLASSID,4); PetscCheckSameComm(svd,1,u,4); }
     324        3385 :   if (v) { PetscValidHeaderSpecific(v,VEC_CLASSID,5); PetscCheckSameComm(svd,1,v,5); }
     325        3385 :   PetscCheck(i>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
     326        3385 :   PetscCheck(i<svd->nconv,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see SVDGetConverged()");
     327        3385 :   if (sigma) *sigma = svd->sigma[svd->perm[i]];
     328        3385 :   if (u || v) {
     329        2285 :     if (!svd->isgeneralized) {
     330        1745 :       PetscCall(MatGetSize(svd->OP,&M,&N));
     331        1745 :       if (M<N) { w = u; u = v; v = w; }
     332             :     }
     333        2285 :     PetscCall(SVDComputeVectors(svd));
     334        2285 :     if (u) PetscCall(BVCopyVec(svd->U,svd->perm[i],u));
     335        2285 :     if (v) PetscCall(BVCopyVec(svd->V,svd->perm[i],v));
     336             :   }
     337        3385 :   PetscFunctionReturn(PETSC_SUCCESS);
     338             : }
     339             : 
     340             : /*
     341             :    SVDComputeResidualNorms_Standard - Computes the norms of the left and
     342             :    right residuals associated with the i-th computed singular triplet.
     343             : 
     344             :    Input Parameters:
     345             :      sigma - singular value
     346             :      u,v   - singular vectors
     347             :      x,y   - two work vectors with the same dimensions as u,v
     348             : */
     349         411 : static PetscErrorCode SVDComputeResidualNorms_Standard(SVD svd,PetscReal sigma,Vec u,Vec v,Vec x,Vec y,PetscReal *norm1,PetscReal *norm2)
     350             : {
     351         411 :   PetscInt       M,N;
     352             : 
     353         411 :   PetscFunctionBegin;
     354             :   /* norm1 = ||A*v-sigma*u||_2 */
     355         411 :   if (norm1) {
     356         411 :     PetscCall(MatMult(svd->OP,v,x));
     357         411 :     PetscCall(VecAXPY(x,-sigma,u));
     358         411 :     PetscCall(VecNorm(x,NORM_2,norm1));
     359             :   }
     360             :   /* norm2 = ||A^T*u-sigma*v||_2 */
     361         411 :   if (norm2) {
     362         411 :     PetscCall(MatGetSize(svd->OP,&M,&N));
     363         411 :     if (M<N) PetscCall(MatMult(svd->A,u,y));
     364         341 :     else PetscCall(MatMult(svd->AT,u,y));
     365         411 :     PetscCall(VecAXPY(y,-sigma,v));
     366         411 :     PetscCall(VecNorm(y,NORM_2,norm2));
     367             :   }
     368         411 :   PetscFunctionReturn(PETSC_SUCCESS);
     369             : }
     370             : 
     371             : /*
     372             :    SVDComputeResidualNorms_Generalized - In GSVD, compute the residual norms
     373             :    norm1 = ||s^2*A'*u-c*B'*B*x||_2 and norm2 = ||c^2*B'*v-s*A'*A*x||_2.
     374             : 
     375             :    Input Parameters:
     376             :      sigma - singular value
     377             :      uv    - left singular vectors [u;v]
     378             :      x     - right singular vector
     379             :      y,z   - two work vectors with the same dimension as x
     380             : */
     381         347 : static PetscErrorCode SVDComputeResidualNorms_Generalized(SVD svd,PetscReal sigma,Vec uv,Vec x,Vec y,Vec z,PetscReal *norm1,PetscReal *norm2)
     382             : {
     383         347 :   Vec            u,v,ax,bx,nest,aux[2];
     384         347 :   PetscReal      c,s;
     385             : 
     386         347 :   PetscFunctionBegin;
     387         347 :   PetscCall(MatCreateVecs(svd->OP,NULL,&u));
     388         347 :   PetscCall(MatCreateVecs(svd->OPb,NULL,&v));
     389         347 :   aux[0] = u;
     390         347 :   aux[1] = v;
     391         347 :   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,aux,&nest));
     392         347 :   PetscCall(VecCopy(uv,nest));
     393             : 
     394         347 :   s = 1.0/PetscSqrtReal(1.0+sigma*sigma);
     395         347 :   c = sigma*s;
     396             : 
     397             :   /* norm1 = ||s^2*A'*u-c*B'*B*x||_2 */
     398         347 :   if (norm1) {
     399         347 :     PetscCall(VecDuplicate(v,&bx));
     400         347 :     PetscCall(MatMultHermitianTranspose(svd->OP,u,z));
     401         347 :     PetscCall(MatMult(svd->OPb,x,bx));
     402         347 :     PetscCall(MatMultHermitianTranspose(svd->OPb,bx,y));
     403         347 :     PetscCall(VecAXPBY(y,s*s,-c,z));
     404         347 :     PetscCall(VecNorm(y,NORM_2,norm1));
     405         347 :     PetscCall(VecDestroy(&bx));
     406             :   }
     407             :   /* norm2 = ||c^2*B'*v-s*A'*A*x||_2 */
     408         347 :   if (norm2) {
     409         347 :     PetscCall(VecDuplicate(u,&ax));
     410         347 :     PetscCall(MatMultHermitianTranspose(svd->OPb,v,z));
     411         347 :     PetscCall(MatMult(svd->OP,x,ax));
     412         347 :     PetscCall(MatMultHermitianTranspose(svd->OP,ax,y));
     413         347 :     PetscCall(VecAXPBY(y,c*c,-s,z));
     414         347 :     PetscCall(VecNorm(y,NORM_2,norm2));
     415         347 :     PetscCall(VecDestroy(&ax));
     416             :   }
     417             : 
     418         347 :   PetscCall(VecDestroy(&v));
     419         347 :   PetscCall(VecDestroy(&u));
     420         347 :   PetscCall(VecDestroy(&nest));
     421         347 :   PetscFunctionReturn(PETSC_SUCCESS);
     422             : }
     423             : 
     424             : /*
     425             :    SVDComputeResidualNorms_Hyperbolic - Computes the norms of the left and
     426             :    right residuals associated with the i-th computed singular triplet.
     427             : 
     428             :    Input Parameters:
     429             :      sigma - singular value
     430             :      sign  - corresponding element of the signature Omega2
     431             :      u,v   - singular vectors
     432             :      x,y,z - three work vectors with the same dimensions as u,v,u
     433             : */
     434         531 : static PetscErrorCode SVDComputeResidualNorms_Hyperbolic(SVD svd,PetscReal sigma,PetscReal sign,Vec u,Vec v,Vec x,Vec y,Vec z,PetscReal *norm1,PetscReal *norm2)
     435             : {
     436         531 :   PetscInt       M,N;
     437             : 
     438         531 :   PetscFunctionBegin;
     439             :   /* norm1 = ||A*v-sigma*u||_2 */
     440         531 :   if (norm1) {
     441         531 :     PetscCall(MatMult(svd->OP,v,x));
     442         531 :     PetscCall(VecAXPY(x,-sigma,u));
     443         531 :     PetscCall(VecNorm(x,NORM_2,norm1));
     444             :   }
     445             :   /* norm2 = ||A^T*Omega*u-sigma*sign*v||_2 */
     446         531 :   if (norm2) {
     447         531 :     PetscCall(MatGetSize(svd->OP,&M,&N));
     448         531 :     PetscCall(VecPointwiseMult(z,u,svd->omega));
     449         531 :     if (M<N) PetscCall(MatMult(svd->A,z,y));
     450         495 :     else PetscCall(MatMult(svd->AT,z,y));
     451         531 :     PetscCall(VecAXPY(y,-sigma*sign,v));
     452         531 :     PetscCall(VecNorm(y,NORM_2,norm2));
     453             :   }
     454         531 :   PetscFunctionReturn(PETSC_SUCCESS);
     455             : }
     456             : 
     457             : /*@
     458             :    SVDComputeError - Computes the error (based on the residual norm) associated
     459             :    with the i-th singular triplet.
     460             : 
     461             :    Collective
     462             : 
     463             :    Input Parameters:
     464             : +  svd  - the singular value solver context
     465             : .  i    - the solution index
     466             : -  type - the type of error to compute
     467             : 
     468             :    Output Parameter:
     469             : .  error - the error
     470             : 
     471             :    Notes:
     472             :    The error can be computed in various ways, all of them based on the residual
     473             :    norm obtained as sqrt(n1^2+n2^2) with n1 = ||A*v-sigma*u||_2 and
     474             :    n2 = ||A^T*u-sigma*v||_2, where sigma is the singular value, u is the left
     475             :    singular vector and v is the right singular vector.
     476             : 
     477             :    In the case of the GSVD, the two components of the residual norm are
     478             :    n1 = ||s^2*A'*u-c*B'*B*x||_2 and n2 = ||c^2*B'*v-s*A'*A*x||_2, where [u;v]
     479             :    are the left singular vectors and x is the right singular vector, with
     480             :    sigma=c/s.
     481             : 
     482             :    Level: beginner
     483             : 
     484             : .seealso: SVDErrorType, SVDSolve()
     485             : @*/
     486        1289 : PetscErrorCode SVDComputeError(SVD svd,PetscInt i,SVDErrorType type,PetscReal *error)
     487             : {
     488        1289 :   PetscReal      sigma,norm1,norm2,c,s;
     489        1289 :   Vec            u=NULL,v=NULL,x=NULL,y=NULL,z=NULL;
     490        1289 :   PetscReal      vecnorm=1.0;
     491             : 
     492        1289 :   PetscFunctionBegin;
     493        1289 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     494        5156 :   PetscValidLogicalCollectiveInt(svd,i,2);
     495        5156 :   PetscValidLogicalCollectiveEnum(svd,type,3);
     496        1289 :   PetscAssertPointer(error,4);
     497        1289 :   SVDCheckSolved(svd,1);
     498             : 
     499             :   /* allocate work vectors */
     500        1289 :   switch (svd->problem_type) {
     501         411 :     case SVD_STANDARD:
     502         411 :       PetscCall(SVDSetWorkVecs(svd,2,2));
     503         411 :       u = svd->workl[0];
     504         411 :       v = svd->workr[0];
     505         411 :       x = svd->workl[1];
     506         411 :       y = svd->workr[1];
     507         411 :       break;
     508         347 :     case SVD_GENERALIZED:
     509         347 :       PetscCall(SVDSetWorkVecs(svd,1,3));
     510         347 :       u = svd->workl[0];
     511         347 :       v = svd->workr[0];
     512         347 :       x = svd->workr[1];
     513         347 :       y = svd->workr[2];
     514         347 :       break;
     515         531 :     case SVD_HYPERBOLIC:
     516         531 :       PetscCall(SVDSetWorkVecs(svd,3,2));
     517         531 :       u = svd->workl[0];
     518         531 :       v = svd->workr[0];
     519         531 :       x = svd->workl[1];
     520         531 :       y = svd->workr[1];
     521         531 :       z = svd->workl[2];
     522         531 :       break;
     523             :   }
     524             : 
     525             :   /* compute residual norm */
     526        1289 :   PetscCall(SVDGetSingularTriplet(svd,i,&sigma,u,v));
     527        1289 :   switch (svd->problem_type) {
     528         411 :     case SVD_STANDARD:
     529         411 :       PetscCall(SVDComputeResidualNorms_Standard(svd,sigma,u,v,x,y,&norm1,&norm2));
     530             :       break;
     531         347 :     case SVD_GENERALIZED:
     532         347 :       PetscCall(SVDComputeResidualNorms_Generalized(svd,sigma,u,v,x,y,&norm1,&norm2));
     533             :       break;
     534         531 :     case SVD_HYPERBOLIC:
     535         531 :       PetscCall(SVDComputeResidualNorms_Hyperbolic(svd,sigma,svd->sign[svd->perm[i]],u,v,x,y,z,&norm1,&norm2));
     536             :       break;
     537             :   }
     538        1289 :   *error = SlepcAbs(norm1,norm2);
     539             : 
     540             :   /* compute 2-norm of eigenvector of the cyclic form */
     541        1289 :   if (type!=SVD_ERROR_ABSOLUTE) {
     542        1276 :     switch (svd->problem_type) {
     543         398 :       case SVD_STANDARD:
     544         398 :         vecnorm = PETSC_SQRT2;
     545         398 :         break;
     546         347 :       case SVD_GENERALIZED:
     547         347 :         PetscCall(VecNorm(v,NORM_2,&vecnorm));
     548         347 :         vecnorm = PetscSqrtReal(1.0+vecnorm*vecnorm);
     549         347 :         break;
     550         531 :       case SVD_HYPERBOLIC:
     551         531 :         PetscCall(VecNorm(u,NORM_2,&vecnorm));
     552         531 :         vecnorm = PetscSqrtReal(1.0+vecnorm*vecnorm);
     553         531 :         break;
     554             :     }
     555             :   }
     556             : 
     557             :   /* compute error */
     558        1276 :   switch (type) {
     559             :     case SVD_ERROR_ABSOLUTE:
     560             :       break;
     561         391 :     case SVD_ERROR_RELATIVE:
     562         391 :       if (svd->isgeneralized) {
     563           0 :         s = 1.0/PetscSqrtReal(1.0+sigma*sigma);
     564           0 :         c = sigma*s;
     565           0 :         norm1 /= c*vecnorm;
     566           0 :         norm2 /= s*vecnorm;
     567           0 :         *error = PetscMax(norm1,norm2);
     568         391 :       } else *error /= sigma*vecnorm;
     569             :       break;
     570         885 :     case SVD_ERROR_NORM:
     571         885 :       if (!svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
     572         885 :       if (svd->isgeneralized && !svd->nrmb) PetscCall(MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb));
     573         885 :       *error /= PetscMax(svd->nrma,svd->nrmb)*vecnorm;
     574         885 :       break;
     575           0 :     default:
     576           0 :       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
     577             :   }
     578        1289 :   PetscFunctionReturn(PETSC_SUCCESS);
     579             : }

Generated by: LCOV version 1.14