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[] = "Tests multiple calls to SVDSolve with equal matrix size.\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 two rectangular bidiagonal matrices
      20             : 
      21             :               |  1  2                     |       |  1                        |
      22             :               |     1  2                  |       |  2  1                     |
      23             :               |        1  2               |       |     2  1                  |
      24             :           A = |          .  .             |   B = |       .  .                |
      25             :               |             .  .          |       |          .  .             |
      26             :               |                1  2       |       |             2  1          |
      27             :               |                   1  2    |       |                2  1       |
      28             :  */
      29             : 
      30          14 : int main(int argc,char **argv)
      31             : {
      32          14 :   Mat            A,B;
      33          14 :   SVD            svd;
      34          14 :   PetscInt       m=20,n,Istart,Iend,i,col[2];
      35          14 :   PetscScalar    valsa[] = { 1, 2 }, valsb[] = { 2, 1 };
      36          14 :   PetscBool      flg;
      37             : 
      38          14 :   PetscFunctionBeginUser;
      39          14 :   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
      40          14 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
      41          14 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,&flg));
      42          14 :   if (!flg) n=m+2;
      43          14 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nRectangular bidiagonal matrix, m=%" PetscInt_FMT " n=%" PetscInt_FMT "\n\n",m,n));
      44             : 
      45             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      46             :                      Generate the matrices
      47             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      48             : 
      49          14 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
      50          14 :   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n));
      51          14 :   PetscCall(MatSetFromOptions(A));
      52          14 :   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
      53         294 :   for (i=Istart;i<Iend;i++) {
      54         280 :     col[0]=i; col[1]=i+1;
      55         280 :     if (i<n-1) PetscCall(MatSetValues(A,1,&i,2,col,valsa,INSERT_VALUES));
      56         280 :     else if (i==n-1) PetscCall(MatSetValue(A,i,col[0],valsa[0],INSERT_VALUES));
      57             :   }
      58          14 :   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
      59          14 :   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
      60             : 
      61          14 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
      62          14 :   PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n));
      63          14 :   PetscCall(MatSetFromOptions(B));
      64          14 :   PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
      65         294 :   for (i=Istart;i<Iend;i++) {
      66         280 :     col[0]=i-1; col[1]=i;
      67         280 :     if (i==0) PetscCall(MatSetValue(B,i,col[1],valsb[1],INSERT_VALUES));
      68         266 :     else if (i<n) PetscCall(MatSetValues(B,1,&i,2,col,valsb,INSERT_VALUES));
      69         280 :     else if (i==n) PetscCall(MatSetValue(B,i,col[0],valsb[0],INSERT_VALUES));
      70             :   }
      71          14 :   PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
      72          14 :   PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
      73             : 
      74             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      75             :          Create the singular value solver, set options and solve
      76             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      77             : 
      78          14 :   PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
      79          14 :   PetscCall(SVDSetOperators(svd,A,NULL));
      80          14 :   PetscCall(SVDSetTolerances(svd,PETSC_DEFAULT,1000));
      81          14 :   PetscCall(SVDSetFromOptions(svd));
      82          14 :   PetscCall(SVDSolve(svd));
      83          14 :   PetscCall(SVDErrorView(svd,SVD_ERROR_RELATIVE,NULL));
      84             : 
      85             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      86             :                        Solve with second matrix
      87             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      88             : 
      89          14 :   PetscCall(SVDSetOperators(svd,B,NULL));
      90          14 :   PetscCall(SVDSolve(svd));
      91          14 :   PetscCall(SVDErrorView(svd,SVD_ERROR_RELATIVE,NULL));
      92             : 
      93             :   /* Free work space */
      94          14 :   PetscCall(SVDDestroy(&svd));
      95          14 :   PetscCall(MatDestroy(&A));
      96          14 :   PetscCall(MatDestroy(&B));
      97          14 :   PetscCall(SlepcFinalize());
      98             :   return 0;
      99             : }
     100             : 
     101             : /*TEST
     102             : 
     103             :    testset:
     104             :       args: -svd_nsv 3
     105             :       requires: !single
     106             :       output_file: output/test14_1.out
     107             :       test:
     108             :          suffix: 1
     109             :          args: -svd_type {{lanczos trlanczos lapack}}
     110             :       test:
     111             :          suffix: 1_cross
     112             :          args: -svd_type cross -svd_cross_explicitmatrix {{0 1}}
     113             :       test:
     114             :          suffix: 1_cyclic
     115             :          args: -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}}
     116             : 
     117             :    testset:
     118             :       args: -n 18 -svd_nsv 3
     119             :       requires: !single
     120             :       output_file: output/test14_2.out
     121             :       test:
     122             :          suffix: 2
     123             :          args: -svd_type {{lanczos trlanczos lapack}}
     124             :       test:
     125             :          suffix: 2_cross
     126             :          args: -svd_type cross -svd_cross_explicitmatrix {{0 1}}
     127             :       test:
     128             :          suffix: 2_cyclic
     129             :          args: -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}}
     130             : 
     131             : TEST*/