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[] = "SVD via the cyclic matrix with a user-provided EPS.\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           1 : int main(int argc,char **argv)
      31             : {
      32           1 :   Mat                  A;
      33           1 :   SVD                  svd;
      34           1 :   EPS                  eps;
      35           1 :   ST                   st;
      36           1 :   KSP                  ksp;
      37           1 :   PC                   pc;
      38           1 :   PetscInt             m=20,n,Istart,Iend,i,col[2];
      39           1 :   PetscScalar          value[] = { 1, 2 };
      40           1 :   PetscBool            flg,expmat;
      41             : 
      42           1 :   PetscFunctionBeginUser;
      43           1 :   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
      44             : 
      45           1 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
      46           1 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,&flg));
      47           1 :   if (!flg) n=m+2;
      48           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nRectangular bidiagonal matrix, m=%" PetscInt_FMT " n=%" PetscInt_FMT "\n\n",m,n));
      49             : 
      50             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      51             :         Generate the matrix
      52             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      53             : 
      54           1 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
      55           1 :   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n));
      56           1 :   PetscCall(MatSetFromOptions(A));
      57           1 :   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
      58          21 :   for (i=Istart;i<Iend;i++) {
      59          20 :     col[0]=i; col[1]=i+1;
      60          20 :     if (i<n-1) PetscCall(MatSetValues(A,1,&i,2,col,value,INSERT_VALUES));
      61          20 :     else if (i==n-1) PetscCall(MatSetValue(A,i,col[0],value[0],INSERT_VALUES));
      62             :   }
      63           1 :   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
      64           1 :   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
      65             : 
      66             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      67             :         Create a standalone EPS with appropriate settings
      68             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      69             : 
      70           1 :   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
      71           1 :   PetscCall(EPSSetWhichEigenpairs(eps,EPS_TARGET_MAGNITUDE));
      72           1 :   PetscCall(EPSSetTarget(eps,1.0));
      73           1 :   PetscCall(EPSGetST(eps,&st));
      74           1 :   PetscCall(STSetType(st,STSINVERT));
      75           1 :   PetscCall(STSetShift(st,1.01));
      76           1 :   PetscCall(STGetKSP(st,&ksp));
      77           1 :   PetscCall(KSPSetType(ksp,KSPPREONLY));
      78           1 :   PetscCall(KSPGetPC(ksp,&pc));
      79           1 :   PetscCall(PCSetType(pc,PCLU));
      80           1 :   PetscCall(EPSSetFromOptions(eps));
      81             : 
      82             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      83             :         Compute singular values
      84             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      85             : 
      86           1 :   PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
      87           1 :   PetscCall(SVDSetOperators(svd,A,NULL));
      88           1 :   PetscCall(SVDSetType(svd,SVDCYCLIC));
      89           1 :   PetscCall(SVDCyclicSetEPS(svd,eps));
      90           1 :   PetscCall(SVDCyclicSetExplicitMatrix(svd,PETSC_TRUE));
      91           1 :   PetscCall(SVDSetWhichSingularTriplets(svd,SVD_SMALLEST));
      92           1 :   PetscCall(SVDSetFromOptions(svd));
      93           1 :   PetscCall(PetscObjectTypeCompare((PetscObject)svd,SVDCYCLIC,&flg));
      94           1 :   if (flg) {
      95           1 :     PetscCall(SVDCyclicGetExplicitMatrix(svd,&expmat));
      96           1 :     if (expmat) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Using explicit matrix with cyclic solver\n"));
      97             :   }
      98           1 :   PetscCall(SVDSolve(svd));
      99             : 
     100             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     101             :                     Display solution and clean up
     102             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     103           1 :   PetscCall(SVDErrorView(svd,SVD_ERROR_ABSOLUTE,PETSC_VIEWER_STDOUT_WORLD));
     104           1 :   PetscCall(SVDDestroy(&svd));
     105           1 :   PetscCall(EPSDestroy(&eps));
     106           1 :   PetscCall(MatDestroy(&A));
     107           1 :   PetscCall(SlepcFinalize());
     108             :   return 0;
     109             : }
     110             : 
     111             : /*TEST
     112             : 
     113             :    test:
     114             :       suffix: 1
     115             :       args: -log_exclude svd
     116             : 
     117             :    test:
     118             :       suffix: 2_cuda
     119             :       args: -log_exclude svd -mat_type aijcusparse
     120             :       requires: cuda
     121             :       output_file: output/test7_1.out
     122             : 
     123             : TEST*/