Actual source code: test14.c

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

      6:    This file is part of SLEPc.
      7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: */

     11: static char help[] = "Test BV created from a Mat.\n\n";

     13: #include <slepcbv.h>

     15: int main(int argc,char **argv)
     16: {
     17:   BV             X;
     18:   Mat            A,B,M;
     19:   PetscInt       i,j,n=20,k=8,Istart,Iend;
     20:   PetscViewer    view;
     21:   PetscBool      sparse=PETSC_FALSE,verbose;
     22:   PetscReal      norm;
     23:   PetscScalar    alpha;

     25:   PetscFunctionBeginUser;
     26:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
     27:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
     28:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
     29:   PetscCall(PetscOptionsHasName(NULL,NULL,"-sparse",&sparse));
     30:   PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
     31:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test BV created from a %s Mat (length %" PetscInt_FMT ", k=%" PetscInt_FMT ").\n",sparse?"sparse":"dense",n,k));

     33:   /* Create matrix */
     34:   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
     35:   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,k));
     36:   if (!sparse) PetscCall(MatSetType(A,MATDENSE));
     37:   else PetscCall(MatSetType(A,MATAIJ));
     38:   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
     39:   for (j=0;j<k;j++) {
     40:     for (i=0;i<=n/2;i++) {
     41:       if (i+j<n && i>=Istart && i<Iend) {
     42:         alpha = (3.0*i+j-2)/(2*(i+j+1));
     43:         PetscCall(MatSetValue(A,i+j,j,alpha,INSERT_VALUES));
     44:       }
     45:     }
     46:   }
     47:   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
     48:   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));

     50:   /* Create BV object X */
     51:   PetscCall(BVCreateFromMat(A,&X));
     52:   PetscCall(BVSetFromOptions(X));
     53:   PetscCall(PetscObjectSetName((PetscObject)X,"X"));

     55:   /* Set up viewer */
     56:   PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view));
     57:   if (verbose) {
     58:     PetscCall(PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB));
     59:     PetscCall(BVView(X,view));
     60:   }

     62:   /* Test BVCreateMat */
     63:   PetscCall(BVCreateMat(X,&B));
     64:   PetscCall(MatAXPY(B,-1.0,A,SAME_NONZERO_PATTERN));
     65:   PetscCall(MatNorm(B,NORM_1,&norm));
     66:   if (norm<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of difference < 100*eps\n"));
     67:   else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of difference: %g\n",(double)norm));

     69:   /* Test BVOrthogonalize */
     70:   PetscCall(BVOrthogonalize(X,NULL));
     71:   if (verbose) PetscCall(BVView(X,view));

     73:   /* Check orthogonality */
     74:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&M));
     75:   PetscCall(MatShift(M,1.0));   /* set leading part to identity */
     76:   PetscCall(BVDot(X,X,M));
     77:   PetscCall(MatShift(M,-1.0));
     78:   PetscCall(MatNorm(M,NORM_1,&norm));
     79:   if (norm<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality < 100*eps\n"));
     80:   else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g\n",(double)norm));

     82:   PetscCall(MatDestroy(&M));
     83:   PetscCall(MatDestroy(&A));
     84:   PetscCall(MatDestroy(&B));
     85:   PetscCall(BVDestroy(&X));
     86:   PetscCall(SlepcFinalize());
     87:   return 0;
     88: }

     90: /*TEST

     92:    test:
     93:       suffix: 1
     94:       nsize: 2
     95:       args: -bv_type {{vecs contiguous svec mat}shared output}

     97:    test:
     98:       suffix: 2
     99:       args: -sparse

    101: TEST*/