LCOV - code coverage report
Current view: top level - svd/interface - svdsolve.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 180 182 98.9 %
Date: 2021-08-02 00:32:28 Functions: 10 10 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-2021, 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             : 
      16             : /*
      17             :   SVDComputeVectors_Left - Compute left singular vectors as U=A*V.
      18             :   Only done if the leftbasis flag is false. Assumes V is available.
      19             :  */
      20          45 : PetscErrorCode SVDComputeVectors_Left(SVD svd)
      21             : {
      22          45 :   PetscErrorCode ierr;
      23          45 :   Vec            tl;
      24          45 :   PetscInt       oldsize;
      25             : 
      26          45 :   PetscFunctionBegin;
      27          45 :   if (!svd->leftbasis) {
      28             :     /* generate left singular vectors on U */
      29          30 :     if (!svd->U) { ierr = SVDGetBV(svd,NULL,&svd->U);CHKERRQ(ierr); }
      30          30 :     ierr = BVGetSizes(svd->U,NULL,NULL,&oldsize);CHKERRQ(ierr);
      31          30 :     if (!oldsize) {
      32          25 :       if (!((PetscObject)(svd->U))->type_name) {
      33           1 :         ierr = BVSetType(svd->U,BVSVEC);CHKERRQ(ierr);
      34             :       }
      35          25 :       ierr = MatCreateVecsEmpty(svd->A,NULL,&tl);CHKERRQ(ierr);
      36          25 :       ierr = BVSetSizesFromVec(svd->U,tl,svd->ncv);CHKERRQ(ierr);
      37          25 :       ierr = VecDestroy(&tl);CHKERRQ(ierr);
      38             :     }
      39          30 :     ierr = BVSetActiveColumns(svd->V,0,svd->nconv);CHKERRQ(ierr);
      40          30 :     ierr = BVSetActiveColumns(svd->U,0,svd->nconv);CHKERRQ(ierr);
      41          30 :     ierr = BVMatMult(svd->V,svd->A,svd->U);CHKERRQ(ierr);
      42          30 :     ierr = BVOrthogonalize(svd->U,NULL);CHKERRQ(ierr);
      43             :   }
      44          45 :   PetscFunctionReturn(0);
      45             : }
      46             : 
      47        1022 : PetscErrorCode SVDComputeVectors(SVD svd)
      48             : {
      49        1022 :   PetscErrorCode ierr;
      50             : 
      51        1022 :   PetscFunctionBegin;
      52        1022 :   SVDCheckSolved(svd,1);
      53        1022 :   if (svd->state==SVD_STATE_SOLVED && svd->ops->computevectors) {
      54          86 :     ierr = (*svd->ops->computevectors)(svd);CHKERRQ(ierr);
      55             :   }
      56        1022 :   svd->state = SVD_STATE_VECTORS;
      57        1022 :   PetscFunctionReturn(0);
      58             : }
      59             : 
      60             : /*@
      61             :    SVDSolve - Solves the singular value problem.
      62             : 
      63             :    Collective on svd
      64             : 
      65             :    Input Parameter:
      66             : .  svd - singular value solver context obtained from SVDCreate()
      67             : 
      68             :    Options Database Keys:
      69             : +  -svd_view - print information about the solver used
      70             : .  -svd_view_mat0 binary - save the first matrix (A) to the default binary viewer
      71             : .  -svd_view_mat1 binary - save the second matrix (B) to the default binary viewer
      72             : .  -svd_view_vectors binary - save the computed singular vectors to the default binary viewer
      73             : .  -svd_view_values - print computed singular values
      74             : .  -svd_converged_reason - print reason for convergence, and number of iterations
      75             : .  -svd_error_absolute - print absolute errors of each singular triplet
      76             : -  -svd_error_relative - print relative errors of each singular triplet
      77             : 
      78             :    Level: beginner
      79             : 
      80             : .seealso: SVDCreate(), SVDSetUp(), SVDDestroy()
      81             : @*/
      82         187 : PetscErrorCode SVDSolve(SVD svd)
      83             : {
      84         187 :   PetscErrorCode ierr;
      85         187 :   PetscInt       i,*workperm;
      86             : 
      87         187 :   PetscFunctionBegin;
      88         187 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
      89         187 :   if (svd->state>=SVD_STATE_SOLVED) PetscFunctionReturn(0);
      90         187 :   ierr = PetscLogEventBegin(SVD_Solve,svd,0,0,0);CHKERRQ(ierr);
      91             : 
      92             :   /* call setup */
      93         187 :   ierr = SVDSetUp(svd);CHKERRQ(ierr);
      94         187 :   svd->its = 0;
      95         187 :   svd->nconv = 0;
      96        2789 :   for (i=0;i<svd->ncv;i++) {
      97        2602 :     svd->sigma[i]  = 0.0;
      98        2602 :     svd->errest[i] = 0.0;
      99        2602 :     svd->perm[i]   = i;
     100             :   }
     101         187 :   ierr = SVDViewFromOptions(svd,NULL,"-svd_view_pre");CHKERRQ(ierr);
     102             : 
     103         187 :   switch (svd->problem_type) {
     104         140 :     case SVD_STANDARD:
     105         140 :       ierr = (*svd->ops->solve)(svd);CHKERRQ(ierr);
     106             :       break;
     107          47 :     case SVD_GENERALIZED:
     108          47 :       ierr = (*svd->ops->solveg)(svd);CHKERRQ(ierr);
     109             :       break;
     110             :   }
     111         187 :   svd->state = SVD_STATE_SOLVED;
     112             : 
     113             :   /* sort singular triplets */
     114         187 :   if (svd->which == SVD_SMALLEST) {
     115          21 :     ierr = PetscSortRealWithPermutation(svd->nconv,svd->sigma,svd->perm);CHKERRQ(ierr);
     116             :   } else {
     117         166 :     ierr = PetscMalloc1(svd->nconv,&workperm);CHKERRQ(ierr);
     118        1427 :     for (i=0;i<svd->nconv;i++) workperm[i] = i;
     119         166 :     ierr = PetscSortRealWithPermutation(svd->nconv,svd->sigma,workperm);CHKERRQ(ierr);
     120        1427 :     for (i=0;i<svd->nconv;i++) svd->perm[i] = workperm[svd->nconv-i-1];
     121         166 :     ierr = PetscFree(workperm);CHKERRQ(ierr);
     122             :   }
     123         187 :   ierr = PetscLogEventEnd(SVD_Solve,svd,0,0,0);CHKERRQ(ierr);
     124             : 
     125             :   /* various viewers */
     126         187 :   ierr = SVDViewFromOptions(svd,NULL,"-svd_view");CHKERRQ(ierr);
     127         187 :   ierr = SVDConvergedReasonViewFromOptions(svd);CHKERRQ(ierr);
     128         187 :   ierr = SVDErrorViewFromOptions(svd);CHKERRQ(ierr);
     129         187 :   ierr = SVDValuesViewFromOptions(svd);CHKERRQ(ierr);
     130         187 :   ierr = SVDVectorsViewFromOptions(svd);CHKERRQ(ierr);
     131         187 :   ierr = MatViewFromOptions(svd->OP,(PetscObject)svd,"-svd_view_mat0");CHKERRQ(ierr);
     132         187 :   if (svd->isgeneralized) {
     133          47 :     ierr = MatViewFromOptions(svd->OPb,(PetscObject)svd,"-svd_view_mat1");CHKERRQ(ierr);
     134             :   }
     135             : 
     136             :   /* Remove the initial subspaces */
     137         187 :   svd->nini = 0;
     138         187 :   svd->ninil = 0;
     139         187 :   PetscFunctionReturn(0);
     140             : }
     141             : 
     142             : /*@
     143             :    SVDGetIterationNumber - Gets the current iteration number. If the
     144             :    call to SVDSolve() is complete, then it returns the number of iterations
     145             :    carried out by the solution method.
     146             : 
     147             :    Not Collective
     148             : 
     149             :    Input Parameter:
     150             : .  svd - the singular value solver context
     151             : 
     152             :    Output Parameter:
     153             : .  its - number of iterations
     154             : 
     155             :    Note:
     156             :    During the i-th iteration this call returns i-1. If SVDSolve() is
     157             :    complete, then parameter "its" contains either the iteration number at
     158             :    which convergence was successfully reached, or failure was detected.
     159             :    Call SVDGetConvergedReason() to determine if the solver converged or
     160             :    failed and why.
     161             : 
     162             :    Level: intermediate
     163             : 
     164             : .seealso: SVDGetConvergedReason(), SVDSetTolerances()
     165             : @*/
     166          47 : PetscErrorCode SVDGetIterationNumber(SVD svd,PetscInt *its)
     167             : {
     168          47 :   PetscFunctionBegin;
     169          47 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     170          47 :   PetscValidIntPointer(its,2);
     171          47 :   *its = svd->its;
     172          47 :   PetscFunctionReturn(0);
     173             : }
     174             : 
     175             : /*@
     176             :    SVDGetConvergedReason - Gets the reason why the SVDSolve() iteration was
     177             :    stopped.
     178             : 
     179             :    Not Collective
     180             : 
     181             :    Input Parameter:
     182             : .  svd - the singular value solver context
     183             : 
     184             :    Output Parameter:
     185             : .  reason - negative value indicates diverged, positive value converged
     186             :    (see SVDConvergedReason)
     187             : 
     188             :    Notes:
     189             : 
     190             :    Possible values for reason are
     191             : +  SVD_CONVERGED_TOL - converged up to tolerance
     192             : .  SVD_CONVERGED_USER - converged due to a user-defined condition
     193             : .  SVD_CONVERGED_MAXIT - reached the maximum number of iterations with SVD_CONV_MAXIT criterion
     194             : .  SVD_DIVERGED_ITS - required more than max_it iterations to reach convergence
     195             : -  SVD_DIVERGED_BREAKDOWN - generic breakdown in method
     196             : 
     197             :    Can only be called after the call to SVDSolve() is complete.
     198             : 
     199             :    Level: intermediate
     200             : 
     201             : .seealso: SVDSetTolerances(), SVDSolve(), SVDConvergedReason
     202             : @*/
     203          13 : PetscErrorCode SVDGetConvergedReason(SVD svd,SVDConvergedReason *reason)
     204             : {
     205          13 :   PetscFunctionBegin;
     206          13 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     207          13 :   PetscValidIntPointer(reason,2);
     208          13 :   SVDCheckSolved(svd,1);
     209          13 :   *reason = svd->reason;
     210          13 :   PetscFunctionReturn(0);
     211             : }
     212             : 
     213             : /*@
     214             :    SVDGetConverged - Gets the number of converged singular values.
     215             : 
     216             :    Not Collective
     217             : 
     218             :    Input Parameter:
     219             : .  svd - the singular value solver context
     220             : 
     221             :    Output Parameter:
     222             : .  nconv - number of converged singular values
     223             : 
     224             :    Note:
     225             :    This function should be called after SVDSolve() has finished.
     226             : 
     227             :    Level: beginner
     228             : 
     229             : @*/
     230          78 : PetscErrorCode SVDGetConverged(SVD svd,PetscInt *nconv)
     231             : {
     232          78 :   PetscFunctionBegin;
     233          78 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     234          78 :   PetscValidIntPointer(nconv,2);
     235          78 :   SVDCheckSolved(svd,1);
     236          78 :   *nconv = svd->nconv;
     237          78 :   PetscFunctionReturn(0);
     238             : }
     239             : 
     240             : /*@C
     241             :    SVDGetSingularTriplet - Gets the i-th triplet of the singular value decomposition
     242             :    as computed by SVDSolve(). The solution consists in the singular value and its left
     243             :    and right singular vectors.
     244             : 
     245             :    Not Collective, but vectors are shared by all processors that share the SVD
     246             : 
     247             :    Input Parameters:
     248             : +  svd - singular value solver context
     249             : -  i   - index of the solution
     250             : 
     251             :    Output Parameters:
     252             : +  sigma - singular value
     253             : .  u     - left singular vector
     254             : -  v     - right singular vector
     255             : 
     256             :    Note:
     257             :    Both u or v can be NULL if singular vectors are not required.
     258             :    Otherwise, the caller must provide valid Vec objects, i.e.,
     259             :    they must be created by the calling program with e.g. MatCreateVecs().
     260             : 
     261             :    The index i should be a value between 0 and nconv-1 (see SVDGetConverged()).
     262             :    Singular triplets are indexed according to the ordering criterion established
     263             :    with SVDSetWhichSingularTriplets().
     264             : 
     265             :    In the case of GSVD, the solution consists in three vectors u,v,x that are
     266             :    returned as follows. Vector x is returned in the right singular vector
     267             :    (argument v) and has length equal to the number of columns of A and B.
     268             :    The other two vectors are returned stacked on top of each other [u;v] in
     269             :    the left singular vector argument, with length equal to m+n (number of rows
     270             :    of A plus number of rows of B).
     271             : 
     272             :    Level: beginner
     273             : 
     274             : .seealso: SVDSolve(), SVDGetConverged(), SVDSetWhichSingularTriplets()
     275             : @*/
     276        1507 : PetscErrorCode SVDGetSingularTriplet(SVD svd,PetscInt i,PetscReal *sigma,Vec u,Vec v)
     277             : {
     278        1507 :   PetscErrorCode ierr;
     279        1507 :   PetscInt       M,N;
     280        1507 :   Vec            w;
     281             : 
     282        1507 :   PetscFunctionBegin;
     283        1507 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     284        3014 :   PetscValidLogicalCollectiveInt(svd,i,2);
     285        1507 :   SVDCheckSolved(svd,1);
     286        1507 :   if (u) { PetscValidHeaderSpecific(u,VEC_CLASSID,4); PetscCheckSameComm(svd,1,u,4); }
     287        1507 :   if (v) { PetscValidHeaderSpecific(v,VEC_CLASSID,5); PetscCheckSameComm(svd,1,v,5); }
     288        1507 :   if (i<0) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
     289        1507 :   if (i>=svd->nconv) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see SVDGetConverged()");
     290        1507 :   if (sigma) *sigma = svd->sigma[svd->perm[i]];
     291        1507 :   if (u || v) {
     292        1021 :     if (!svd->isgeneralized) {
     293         691 :       ierr = MatGetSize(svd->OP,&M,&N);CHKERRQ(ierr);
     294         691 :       if (M<N) { w = u; u = v; v = w; }
     295             :     }
     296        1021 :     ierr = SVDComputeVectors(svd);CHKERRQ(ierr);
     297        1021 :     if (u) { ierr = BVCopyVec(svd->U,svd->perm[i],u);CHKERRQ(ierr); }
     298        1021 :     if (v) { ierr = BVCopyVec(svd->V,svd->perm[i],v);CHKERRQ(ierr); }
     299             :   }
     300        1507 :   PetscFunctionReturn(0);
     301             : }
     302             : 
     303             : /*
     304             :    SVDComputeResidualNorms_Standard - Computes the norms of the left and
     305             :    right residuals associated with the i-th computed singular triplet.
     306             : 
     307             :    Input Parameters:
     308             :      sigma - singular value
     309             :      u,v   - singular vectors
     310             :      x,y   - two work vectors with the same dimensions as u,v
     311             : @*/
     312         418 : static PetscErrorCode SVDComputeResidualNorms_Standard(SVD svd,PetscReal sigma,Vec u,Vec v,Vec x,Vec y,PetscReal *norm1,PetscReal *norm2)
     313             : {
     314         418 :   PetscErrorCode ierr;
     315         418 :   PetscInt       M,N;
     316             : 
     317         418 :   PetscFunctionBegin;
     318             :   /* norm1 = ||A*v-sigma*u||_2 */
     319         418 :   if (norm1) {
     320         418 :     ierr = MatMult(svd->OP,v,x);CHKERRQ(ierr);
     321         418 :     ierr = VecAXPY(x,-sigma,u);CHKERRQ(ierr);
     322         418 :     ierr = VecNorm(x,NORM_2,norm1);CHKERRQ(ierr);
     323             :   }
     324             :   /* norm2 = ||A^T*u-sigma*v||_2 */
     325         418 :   if (norm2) {
     326         418 :     ierr = MatGetSize(svd->OP,&M,&N);CHKERRQ(ierr);
     327         418 :     if (M<N) {
     328          67 :       ierr = MatMult(svd->A,u,y);CHKERRQ(ierr);
     329             :     } else {
     330         351 :       ierr = MatMult(svd->AT,u,y);CHKERRQ(ierr);
     331             :     }
     332         418 :     ierr = VecAXPY(y,-sigma,v);CHKERRQ(ierr);
     333         418 :     ierr = VecNorm(y,NORM_2,norm2);CHKERRQ(ierr);
     334             :   }
     335         418 :   PetscFunctionReturn(0);
     336             : }
     337             : 
     338             : /*
     339             :    SVDComputeResidualNorms_Generalized - In GSVD, compute the residual norms
     340             :    norm1 = ||A*x-c*u||_2 and norm2 = ||B*x-s*v||_2.
     341             : 
     342             :    Input Parameters:
     343             :      sigma - singular value
     344             :      x     - right singular vector
     345             : */
     346         197 : static PetscErrorCode SVDComputeResidualNorms_Generalized(SVD svd,PetscReal sigma,Vec uv,Vec x,PetscReal *norm1,PetscReal *norm2)
     347             : {
     348         197 :   PetscErrorCode ierr;
     349         197 :   Vec            u,v,ax,bx,nest,aux[2];
     350         197 :   PetscReal      c,s;
     351             : 
     352         197 :   PetscFunctionBegin;
     353         197 :   ierr = MatCreateVecs(svd->OP,NULL,&u);CHKERRQ(ierr);
     354         197 :   ierr = MatCreateVecs(svd->OPb,NULL,&v);CHKERRQ(ierr);
     355         197 :   aux[0] = u;
     356         197 :   aux[1] = v;
     357         197 :   ierr = VecCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,aux,&nest);CHKERRQ(ierr);
     358         197 :   ierr = VecCopy(uv,nest);CHKERRQ(ierr);
     359         197 :   ierr = VecDuplicate(u,&ax);CHKERRQ(ierr);
     360         197 :   ierr = VecDuplicate(v,&bx);CHKERRQ(ierr);
     361             : 
     362         197 :   s = 1.0/PetscSqrtReal(1.0+sigma*sigma);
     363         197 :   c = sigma*s;
     364             : 
     365             :   /* norm1 = ||A*x-c*u||_2 */
     366         197 :   if (norm1) {
     367         197 :     ierr = MatMult(svd->OP,x,ax);CHKERRQ(ierr);
     368         197 :     ierr = VecAXPY(ax,-c,u);CHKERRQ(ierr);
     369         197 :     ierr = VecNorm(ax,NORM_2,norm1);CHKERRQ(ierr);
     370             :   }
     371             :   /* norm2 = ||B*x-s*v||_2 */
     372         197 :   if (norm2) {
     373         197 :     ierr = MatMult(svd->OPb,x,bx);CHKERRQ(ierr);
     374         197 :     ierr = VecAXPY(bx,-s,v);CHKERRQ(ierr);
     375         197 :     ierr = VecNorm(bx,NORM_2,norm2);CHKERRQ(ierr);
     376             :   }
     377             : 
     378         197 :   ierr = VecDestroy(&v);CHKERRQ(ierr);
     379         197 :   ierr = VecDestroy(&u);CHKERRQ(ierr);
     380         197 :   ierr = VecDestroy(&nest);CHKERRQ(ierr);
     381         197 :   ierr = VecDestroy(&ax);CHKERRQ(ierr);
     382         197 :   ierr = VecDestroy(&bx);CHKERRQ(ierr);
     383         197 :   PetscFunctionReturn(0);
     384             : }
     385             : 
     386             : /*@
     387             :    SVDComputeError - Computes the error (based on the residual norm) associated
     388             :    with the i-th singular triplet.
     389             : 
     390             :    Collective on svd
     391             : 
     392             :    Input Parameters:
     393             : +  svd  - the singular value solver context
     394             : .  i    - the solution index
     395             : -  type - the type of error to compute
     396             : 
     397             :    Output Parameter:
     398             : .  error - the error
     399             : 
     400             :    Notes:
     401             :    The error can be computed in various ways, all of them based on the residual
     402             :    norm obtained as sqrt(n1^2+n2^2) with n1 = ||A*v-sigma*u||_2 and
     403             :    n2 = ||A^T*u-sigma*v||_2, where sigma is the singular value, u is the left
     404             :    singular vector and v is the right singular vector.
     405             : 
     406             :    Level: beginner
     407             : 
     408             : .seealso: SVDErrorType, SVDSolve()
     409             : @*/
     410         615 : PetscErrorCode SVDComputeError(SVD svd,PetscInt i,SVDErrorType type,PetscReal *error)
     411             : {
     412         615 :   PetscErrorCode ierr;
     413         615 :   PetscReal      sigma,norm1,norm2;
     414         615 :   Vec            u,v,x,y;
     415             : 
     416         615 :   PetscFunctionBegin;
     417         615 :   PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
     418        1230 :   PetscValidLogicalCollectiveInt(svd,i,2);
     419        1230 :   PetscValidLogicalCollectiveEnum(svd,type,3);
     420         615 :   PetscValidRealPointer(error,4);
     421         615 :   SVDCheckSolved(svd,1);
     422             : 
     423             :   /* allocate work vectors */
     424         615 :   ierr = SVDSetWorkVecs(svd,2,2);CHKERRQ(ierr);
     425         615 :   u = svd->workl[0];
     426         615 :   v = svd->workr[0];
     427         615 :   x = svd->workl[1];
     428         615 :   y = svd->workr[1];
     429             : 
     430             :   /* compute residual norm and error */
     431         615 :   ierr = SVDGetSingularTriplet(svd,i,&sigma,u,v);CHKERRQ(ierr);
     432         615 :   switch (svd->problem_type) {
     433         418 :     case SVD_STANDARD:
     434         418 :       ierr = SVDComputeResidualNorms_Standard(svd,sigma,u,v,x,y,&norm1,&norm2);CHKERRQ(ierr);
     435             :       break;
     436         197 :     case SVD_GENERALIZED:
     437         197 :       ierr = SVDComputeResidualNorms_Generalized(svd,sigma,u,v,&norm1,&norm2);CHKERRQ(ierr);
     438             :       break;
     439             :   }
     440         615 :   *error = PetscSqrtReal(norm1*norm1+norm2*norm2);
     441         615 :   switch (type) {
     442             :     case SVD_ERROR_ABSOLUTE:
     443             :       break;
     444         603 :     case SVD_ERROR_RELATIVE:
     445         603 :       *error /= sigma;
     446         603 :       break;
     447           0 :     default:
     448           0 :       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
     449             :   }
     450         615 :   PetscFunctionReturn(0);
     451             : }
     452             : 

Generated by: LCOV version 1.14