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

      1: static char help[] = "Testing PtAP for SeqMAIJ matrix, P, with SeqAIJ matrix, A.\n\n";

      3: #include <petscmat.h>

      7: int main(int argc,char **argv)
      8: {
      9:   Mat            pA,P,aijP;
     10:   PetscScalar    pa[]={1.,-1.,0.,0.,1.,-1.,0.,0.,1.};
     11:   PetscInt       i,pij[]={0,1,2};
     12:   PetscInt       aij[3][3]={{0,1,2},{3,4,5},{6,7,8}};
     13:   Mat            A,mC,C;
     14:   PetscScalar    one=1.;
     16:   PetscMPIInt    size;

     18:   PetscInitialize(&argc,&argv,(char*)0,help);
     19:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     20:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");

     22:   /* Create MAIJ matrix, P */
     23:   MatCreate(PETSC_COMM_SELF,&pA);
     24:   MatSetSizes(pA,3,3,3,3);
     25:   MatSetType(pA,MATSEQAIJ);
     26:   MatSetUp(pA);
     27:   MatSetOption(pA,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
     28:   MatSetValues(pA,3,pij,3,pij,pa,ADD_VALUES);
     29:   MatAssemblyBegin(pA,MAT_FINAL_ASSEMBLY);
     30:   MatAssemblyEnd(pA,MAT_FINAL_ASSEMBLY);
     31:   MatCreateMAIJ(pA,3,&P);
     32:   MatDestroy(&pA);

     34:   /* Create AIJ equivalent matrix, aijP, for comparison testing */
     35:   MatConvert(P,MATSEQAIJ,MAT_INITIAL_MATRIX,&aijP);

     37:   /* Create AIJ matrix, A */
     38:   MatCreate(PETSC_COMM_SELF,&A);
     39:   MatSetSizes(A,9,9,9,9);
     40:   MatSetType(A,MATSEQAIJ);
     41:   MatSetUp(A);
     42:   MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
     43:   MatSetValues(A,3,aij[0],3,aij[0],pa,ADD_VALUES);
     44:   MatSetValues(A,3,aij[1],3,aij[1],pa,ADD_VALUES);
     45:   MatSetValues(A,3,aij[2],3,aij[2],pa,ADD_VALUES);
     46:   for (i=0; i<9; i++) {
     47:     MatSetValue(A,i,i,one,ADD_VALUES);
     48:   }
     49:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
     50:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

     52:   /* Perform PtAP_SeqAIJ_SeqMAIJ */
     53:   MatPtAP(A,P,MAT_INITIAL_MATRIX,1.,&mC);
     54:   MatPtAP(A,P,MAT_REUSE_MATRIX,1.,&mC);
     55:   MatView(mC,PETSC_VIEWER_STDOUT_SELF);

     57:   /* Perform PtAP_SeqAIJ_SeqAIJ for comparison testing */
     58:   MatPtAP(A,aijP,MAT_INITIAL_MATRIX,1.,&C);
     59:   MatView(C,PETSC_VIEWER_STDOUT_SELF);

     61:   /* Check mC = C */
     62:   MatAXPY(C,-1.0,mC,DIFFERENT_NONZERO_PATTERN);
     63:   /* Note: We should be able to use SAME_NONZERO_PATTERN on the line above, */
     64:   /*       but don't because this flag doesn't assist testing. */
     65:   MatView(C,PETSC_VIEWER_STDOUT_SELF);

     67:   /* Cleanup */
     68:   MatDestroy(&P);
     69:   MatDestroy(&aijP);
     70:   MatDestroy(&A);
     71:   MatDestroy(&C);
     72:   MatDestroy(&mC);
     73:   PetscFinalize();
     74:   return(0);
     75: }