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             : static char help[] = "Lanczos SVD. Also illustrates the use of SVDSetBV().\n\n"
      12             :   "The command line options are:\n"
      13             :   "  -m <m>, where <m> = matrix rows.\n"
      14             :   "  -n <n>, where <n> = matrix columns (defaults to m+2).\n\n";
      15             : 
      16             : #include <slepcsvd.h>
      17             : 
      18             : /*
      19             :    This example computes the singular values of a rectangular bidiagonal matrix
      20             : 
      21             :               |  1  2                     |
      22             :               |     1  2                  |
      23             :               |        1  2               |
      24             :           A = |          .  .             |
      25             :               |             .  .          |
      26             :               |                1  2       |
      27             :               |                   1  2    |
      28             :  */
      29             : 
      30           3 : int main(int argc,char **argv)
      31             : {
      32           3 :   Mat            A;
      33           3 :   SVD            svd;
      34           3 :   PetscInt       m=20,n,Istart,Iend,i,k=6,col[2];
      35           3 :   PetscScalar    value[] = { 1, 2 };
      36           3 :   PetscBool      flg,oneside=PETSC_FALSE;
      37           3 :   const char     *prefix;
      38           3 :   BV             U,V;
      39           3 :   Vec            u,v;
      40             : 
      41           3 :   PetscFunctionBeginUser;
      42           3 :   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
      43           3 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
      44           3 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,&flg));
      45           3 :   if (!flg) n=m+2;
      46           3 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
      47           3 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nRectangular bidiagonal matrix, m=%" PetscInt_FMT " n=%" PetscInt_FMT "\n\n",m,n));
      48             : 
      49             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      50             :         Generate the matrix
      51             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      52             : 
      53           3 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
      54           3 :   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n));
      55           3 :   PetscCall(MatSetFromOptions(A));
      56           3 :   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
      57          65 :   for (i=Istart;i<Iend;i++) {
      58          62 :     col[0]=i; col[1]=i+1;
      59          62 :     if (i<n-1) PetscCall(MatSetValues(A,1,&i,2,col,value,INSERT_VALUES));
      60          62 :     else if (i==n-1) PetscCall(MatSetValue(A,i,col[0],value[0],INSERT_VALUES));
      61             :   }
      62           3 :   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
      63           3 :   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
      64           3 :   PetscCall(MatCreateVecs(A,&v,&u));
      65             : 
      66             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      67             :          Create standalone BV objects to illustrate use of SVDSetBV()
      68             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      69             : 
      70           3 :   PetscCall(BVCreate(PETSC_COMM_WORLD,&U));
      71           3 :   PetscCall(PetscObjectSetName((PetscObject)U,"U"));
      72           3 :   PetscCall(BVSetSizesFromVec(U,u,k));
      73           3 :   PetscCall(BVSetFromOptions(U));
      74           3 :   PetscCall(BVCreate(PETSC_COMM_WORLD,&V));
      75           3 :   PetscCall(PetscObjectSetName((PetscObject)V,"V"));
      76           3 :   PetscCall(BVSetSizesFromVec(V,v,k));
      77           3 :   PetscCall(BVSetFromOptions(V));
      78             : 
      79             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      80             :         Compute singular values
      81             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      82             : 
      83           3 :   PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
      84           3 :   PetscCall(SVDSetBV(svd,V,U));
      85           3 :   PetscCall(SVDSetOptionsPrefix(svd,"check_"));
      86           3 :   PetscCall(SVDAppendOptionsPrefix(svd,"myprefix_"));
      87           3 :   PetscCall(SVDGetOptionsPrefix(svd,&prefix));
      88           3 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"SVD prefix is currently: %s\n\n",prefix));
      89           3 :   PetscCall(PetscObjectSetName((PetscObject)svd,"SVD_solver"));
      90             : 
      91           3 :   PetscCall(SVDSetOperators(svd,A,NULL));
      92           3 :   PetscCall(SVDSetType(svd,SVDLANCZOS));
      93           3 :   PetscCall(SVDSetFromOptions(svd));
      94             : 
      95           3 :   PetscCall(PetscObjectTypeCompare((PetscObject)svd,SVDLANCZOS,&flg));
      96           3 :   if (flg) {
      97           2 :     PetscCall(SVDLanczosGetOneSide(svd,&oneside));
      98           4 :     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Running Lanczos %s\n\n",oneside?"(onesided)":""));
      99             :   }
     100           3 :   PetscCall(PetscObjectTypeCompare((PetscObject)svd,SVDTRLANCZOS,&flg));
     101           3 :   if (flg) {
     102           1 :     PetscCall(SVDTRLanczosGetOneSide(svd,&oneside));
     103           2 :     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Running thick-restart Lanczos %s\n\n",oneside?"(onesided)":""));
     104             :   }
     105             : 
     106           3 :   PetscCall(SVDSolve(svd));
     107             : 
     108             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     109             :                     Display solution and clean up
     110             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     111           3 :   PetscCall(SVDErrorView(svd,SVD_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
     112           3 :   PetscCall(BVDestroy(&U));
     113           3 :   PetscCall(BVDestroy(&V));
     114           3 :   PetscCall(VecDestroy(&u));
     115           3 :   PetscCall(VecDestroy(&v));
     116           3 :   PetscCall(SVDDestroy(&svd));
     117           3 :   PetscCall(MatDestroy(&A));
     118           3 :   PetscCall(SlepcFinalize());
     119             :   return 0;
     120             : }
     121             : 
     122             : /*TEST
     123             : 
     124             :    testset:
     125             :       args: -check_myprefix_svd_nsv 3
     126             :       requires: double
     127             :       test:
     128             :          suffix: 1
     129             :          args: -check_myprefix_svd_view_vectors ::ascii_info
     130             :       test:
     131             :          suffix: 2
     132             :          args: -check_myprefix_svd_type trlanczos -check_myprefix_svd_monitor -check_myprefix_svd_view_values ::ascii_matlab
     133             :          filter: sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
     134             :       test:
     135             :          suffix: 3
     136             :          args: -m 22 -n 20
     137             : 
     138             : TEST*/