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[] = "Test BV created from a Mat.\n\n";
      12             : 
      13             : #include <slepcbv.h>
      14             : 
      15           9 : int main(int argc,char **argv)
      16             : {
      17           9 :   BV             X;
      18           9 :   Mat            A,B,M;
      19           9 :   PetscInt       i,j,n=20,k=8,Istart,Iend;
      20           9 :   PetscViewer    view;
      21           9 :   PetscBool      sparse=PETSC_FALSE,verbose;
      22           9 :   PetscReal      norm;
      23           9 :   PetscScalar    alpha;
      24             : 
      25           9 :   PetscFunctionBeginUser;
      26           9 :   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
      27           9 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
      28           9 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
      29           9 :   PetscCall(PetscOptionsHasName(NULL,NULL,"-sparse",&sparse));
      30           9 :   PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
      31          17 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test BV created from a %s Mat (length %" PetscInt_FMT ", k=%" PetscInt_FMT ").\n",sparse?"sparse":"dense",n,k));
      32             : 
      33             :   /* Create matrix */
      34           9 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
      35           9 :   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,k));
      36           9 :   if (!sparse) PetscCall(MatSetType(A,MATDENSE));
      37           1 :   else PetscCall(MatSetType(A,MATAIJ));
      38           9 :   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
      39          81 :   for (j=0;j<k;j++) {
      40         864 :     for (i=0;i<=n/2;i++) {
      41         792 :       if (i+j<n && i>=Istart && i<Iend) {
      42         440 :         alpha = (3.0*i+j-2)/(2*(i+j+1));
      43         792 :         PetscCall(MatSetValue(A,i+j,j,alpha,INSERT_VALUES));
      44             :       }
      45             :     }
      46             :   }
      47           9 :   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
      48           9 :   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
      49             : 
      50             :   /* Create BV object X */
      51           9 :   PetscCall(BVCreateFromMat(A,&X));
      52           9 :   PetscCall(BVSetFromOptions(X));
      53           9 :   PetscCall(PetscObjectSetName((PetscObject)X,"X"));
      54             : 
      55             :   /* Set up viewer */
      56           9 :   PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view));
      57           9 :   if (verbose) {
      58           0 :     PetscCall(PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB));
      59           0 :     PetscCall(BVView(X,view));
      60             :   }
      61             : 
      62             :   /* Test BVCreateMat */
      63           9 :   PetscCall(BVCreateMat(X,&B));
      64           9 :   PetscCall(MatAXPY(B,-1.0,A,SAME_NONZERO_PATTERN));
      65           9 :   PetscCall(MatNorm(B,NORM_1,&norm));
      66           9 :   if (norm<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of difference < 100*eps\n"));
      67           0 :   else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of difference: %g\n",(double)norm));
      68             : 
      69             :   /* Test BVOrthogonalize */
      70           9 :   PetscCall(BVOrthogonalize(X,NULL));
      71           9 :   if (verbose) PetscCall(BVView(X,view));
      72             : 
      73             :   /* Check orthogonality */
      74           9 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&M));
      75           9 :   PetscCall(MatShift(M,1.0));   /* set leading part to identity */
      76           9 :   PetscCall(BVDot(X,X,M));
      77           9 :   PetscCall(MatShift(M,-1.0));
      78           9 :   PetscCall(MatNorm(M,NORM_1,&norm));
      79           9 :   if (norm<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality < 100*eps\n"));
      80           0 :   else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g\n",(double)norm));
      81             : 
      82           9 :   PetscCall(MatDestroy(&M));
      83           9 :   PetscCall(MatDestroy(&A));
      84           9 :   PetscCall(MatDestroy(&B));
      85           9 :   PetscCall(BVDestroy(&X));
      86           9 :   PetscCall(SlepcFinalize());
      87             :   return 0;
      88             : }
      89             : 
      90             : /*TEST
      91             : 
      92             :    test:
      93             :       suffix: 1
      94             :       nsize: 2
      95             :       args: -bv_type {{vecs contiguous svec mat}shared output}
      96             : 
      97             :    test:
      98             :       suffix: 2
      99             :       args: -sparse
     100             : 
     101             : TEST*/