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 a user-defined convergence test (based on ex8.c).\n\n"
      12             :   "The command line options are:\n"
      13             :   "  -n <n>, where <n> = matrix dimension.\n\n";
      14             : 
      15             : #include <slepcsvd.h>
      16             : 
      17             : /*
      18             :    This example computes the singular values of an nxn Grcar matrix,
      19             :    which is a nonsymmetric Toeplitz matrix:
      20             : 
      21             :               |  1  1  1  1               |
      22             :               | -1  1  1  1  1            |
      23             :               |    -1  1  1  1  1         |
      24             :               |       .  .  .  .  .       |
      25             :           A = |          .  .  .  .  .    |
      26             :               |            -1  1  1  1  1 |
      27             :               |               -1  1  1  1 |
      28             :               |                  -1  1  1 |
      29             :               |                     -1  1 |
      30             : 
      31             :  */
      32             : 
      33             : /*
      34             :   MyConvergedRel - Convergence test relative to the norm of A (given in ctx).
      35             : */
      36          14 : PetscErrorCode MyConvergedRel(SVD svd,PetscReal sigma,PetscReal res,PetscReal *errest,void *ctx)
      37             : {
      38          14 :   PetscReal norm = *(PetscReal*)ctx;
      39             : 
      40          14 :   PetscFunctionBegin;
      41          14 :   *errest = res/norm;
      42          14 :   PetscFunctionReturn(PETSC_SUCCESS);
      43             : }
      44             : 
      45           1 : int main(int argc,char **argv)
      46             : {
      47           1 :   Mat            A;               /* Grcar matrix */
      48           1 :   SVD            svd;             /* singular value solver context */
      49           1 :   PetscInt       N=30,Istart,Iend,i,col[5],nconv1,nconv2;
      50           1 :   PetscScalar    value[] = { -1, 1, 1, 1, 1 };
      51           1 :   PetscReal      sigma_1,sigma_n;
      52             : 
      53           1 :   PetscFunctionBeginUser;
      54           1 :   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
      55             : 
      56           1 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL));
      57           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nEstimate the condition number of a Grcar matrix, n=%" PetscInt_FMT "\n\n",N));
      58             : 
      59             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      60             :         Generate the matrix
      61             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      62             : 
      63           1 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
      64           1 :   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
      65           1 :   PetscCall(MatSetFromOptions(A));
      66           1 :   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
      67          31 :   for (i=Istart;i<Iend;i++) {
      68          30 :     col[0]=i-1; col[1]=i; col[2]=i+1; col[3]=i+2; col[4]=i+3;
      69          30 :     if (i==0) PetscCall(MatSetValues(A,1,&i,PetscMin(4,N-i),col+1,value+1,INSERT_VALUES));
      70          30 :     else PetscCall(MatSetValues(A,1,&i,PetscMin(5,N-i+1),col,value,INSERT_VALUES));
      71             :   }
      72           1 :   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
      73           1 :   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
      74             : 
      75             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      76             :              Create the SVD solver and set the solution method
      77             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      78             : 
      79           1 :   PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
      80           1 :   PetscCall(SVDSetOperators(svd,A,NULL));
      81           1 :   PetscCall(SVDSetType(svd,SVDTRLANCZOS));
      82           1 :   PetscCall(SVDSetFromOptions(svd));
      83             : 
      84             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      85             :                       Solve the singular value problem
      86             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      87             : 
      88           1 :   PetscCall(SVDSetWhichSingularTriplets(svd,SVD_LARGEST));
      89           1 :   PetscCall(SVDSolve(svd));
      90           1 :   PetscCall(SVDGetConverged(svd,&nconv1));
      91           1 :   if (nconv1 > 0) PetscCall(SVDGetSingularTriplet(svd,0,&sigma_1,NULL,NULL));
      92           0 :   else PetscCall(PetscPrintf(PETSC_COMM_WORLD," Unable to compute large singular value!\n\n"));
      93             : 
      94             :   /* compute smallest singular value relative to the matrix norm */
      95           1 :   PetscCall(SVDSetConvergenceTestFunction(svd,MyConvergedRel,&sigma_1,NULL));
      96           1 :   PetscCall(SVDSetWhichSingularTriplets(svd,SVD_SMALLEST));
      97           1 :   PetscCall(SVDSolve(svd));
      98           1 :   PetscCall(SVDGetConverged(svd,&nconv2));
      99           1 :   if (nconv2 > 0) PetscCall(SVDGetSingularTriplet(svd,0,&sigma_n,NULL,NULL));
     100           0 :   else PetscCall(PetscPrintf(PETSC_COMM_WORLD," Unable to compute small singular value!\n\n"));
     101             : 
     102             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     103             :                     Display solution and clean up
     104             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     105           1 :   if (nconv1 > 0 && nconv2 > 0) {
     106           1 :     PetscCall(PetscPrintf(PETSC_COMM_WORLD," Computed singular values: sigma_1=%.4f, sigma_n=%.4f\n",(double)sigma_1,(double)sigma_n));
     107           1 :     PetscCall(PetscPrintf(PETSC_COMM_WORLD," Estimated condition number: sigma_1/sigma_n=%.4f\n\n",(double)(sigma_1/sigma_n)));
     108             :   }
     109             : 
     110           1 :   PetscCall(SVDDestroy(&svd));
     111           1 :   PetscCall(MatDestroy(&A));
     112           1 :   PetscCall(SlepcFinalize());
     113             :   return 0;
     114             : }
     115             : 
     116             : /*TEST
     117             : 
     118             :    test:
     119             :       suffix: 1
     120             : 
     121             : TEST*/