Actual source code: ex47.c
petsc-3.4.2 2013-07-02

      2: static char help[] = "Tests the various routines in MatBAIJ format.\n\
      3: Input arguments are:\n\
      4:   -f <input_file> : file to load.  For example see $PETSC_DIR/share/petsc/datafiles/matrices\n\n";

      6: #include <petscmat.h>

     10: int main(int argc,char **args)
     11: {
     12:   Mat               A,B,C;
     13:   PetscViewer       va,vb,vc;
     14:   Vec               x,y;
     15:   PetscErrorCode    ierr;
     16:   PetscInt          i,j,row,m,n,ncols1,ncols2,ct,m2,n2;
     17:   const PetscInt    *cols1,*cols2;
     18:   char              file[PETSC_MAX_PATH_LEN];
     19:   PetscBool         tflg;
     20:   PetscScalar       rval;
     21:   const PetscScalar *vals1,*vals2;
     22:   PetscReal         norm1,norm2,rnorm;
     23:   PetscRandom       r;

     26:   PetscInitialize(&argc,&args,(char*)0,help);
     27: #if defined(PETSC_USE_COMPLEX)
     28:   SETERRQ(PETSC_COMM_WORLD,1,"This example does not work with complex numbers");
     29: #else

     31:   PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,NULL);

     33:   /* Load the matrix as AIJ format */
     34:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&va);
     35:   MatCreate(PETSC_COMM_WORLD,&A);
     36:   MatSetType(A,MATSEQAIJ);
     37:   MatLoad(A,va);
     38:   PetscViewerDestroy(&va);

     40:   /* Load the matrix as BAIJ format */
     41:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&vb);
     42:   MatCreate(PETSC_COMM_WORLD,&B);
     43:   MatSetType(B,MATSEQBAIJ);
     44:   MatLoad(B,vb);
     45:   PetscViewerDestroy(&vb);

     47:   /* Load the matrix as BAIJ format */
     48:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&vc);
     49:   MatCreate(PETSC_COMM_WORLD,&C);
     50:   MatSetType(C,MATSEQBAIJ);
     51:   MatLoad(C,vc);
     52:   PetscViewerDestroy(&vc);

     54:   MatGetSize(A,&m,&n);
     55:   MatGetSize(B,&m2,&n2);
     56:   if (m!=m2) SETERRQ(PETSC_COMM_SELF,1,"Matrices are of different size. Cannot run this example");

     58:   /* Test MatEqual() */
     59:   MatEqual(B,C,&tflg);
     60:   if (!tflg) SETERRQ(PETSC_COMM_SELF,1,"MatEqual() failed");

     62:   /* Test MatGetDiagonal() */
     63:   VecCreateSeq(PETSC_COMM_SELF,m,&x);
     64:   VecCreateSeq(PETSC_COMM_SELF,m,&y);

     66:   MatGetDiagonal(A,x);
     67:   MatGetDiagonal(B,y);

     69:   VecEqual(x,y,&tflg);
     70:   if (!tflg) SETERRQ(PETSC_COMM_SELF,1,"MatGetDiagonal() failed");

     72:   /* Test MatDiagonalScale() */
     73:   PetscRandomCreate(PETSC_COMM_SELF,&r);
     74:   PetscRandomSetFromOptions(r);
     75:   VecSetRandom(x,r);
     76:   VecSetRandom(y,r);

     78:   MatDiagonalScale(A,x,y);
     79:   MatDiagonalScale(B,x,y);
     80:   MatMult(A,x,y);
     81:   VecNorm(y,NORM_2,&norm1);
     82:   MatMult(B,x,y);
     83:   VecNorm(y,NORM_2,&norm2);
     84:   rnorm = ((norm1-norm2)*100)/norm1;
     85:   if (rnorm<-0.1 || rnorm>0.01) {
     86:     PetscPrintf(PETSC_COMM_SELF,"Norm1=%e Norm2=%e\n",norm1,norm2);
     87:     SETERRQ(PETSC_COMM_SELF,1,"MatDiagonalScale() failed");
     88:   }

     90:   /* Test MatGetRow()/ MatRestoreRow() */
     91:   for (ct=0; ct<100; ct++) {
     92:     PetscRandomGetValue(r,&rval);
     93:     row  = (int)(rval*m);
     94:     MatGetRow(A,row,&ncols1,&cols1,&vals1);
     95:     MatGetRow(B,row,&ncols2,&cols2,&vals2);

     97:     for (i=0,j=0; i<ncols1 && j<ncols2; i++) {
     98:       while (cols2[j] != cols1[i]) j++;
     99:       if (vals1[i] != vals2[j]) SETERRQ(PETSC_COMM_SELF,1,"MatGetRow() failed - vals incorrect.");
    100:     }
    101:     if (i<ncols1) SETERRQ(PETSC_COMM_SELF,1,"MatGetRow() failed - cols incorrect");

    103:     MatRestoreRow(A,row,&ncols1,&cols1,&vals1);
    104:     MatRestoreRow(B,row,&ncols2,&cols2,&vals2);
    105:   }

    107:   MatDestroy(&A);
    108:   MatDestroy(&B);
    109:   MatDestroy(&C);
    110:   VecDestroy(&x);
    111:   VecDestroy(&y);
    112:   PetscRandomDestroy(&r);
    113:   PetscFinalize();
    114: #endif
    115:   return 0;
    116: }