Actual source code: ex101.c

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

  3: #include <petscmat.h>

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

 16:   PetscInitialize(&argc,&argv,(char*)0,help);
 17:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

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

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

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

 50:   /* Perform PtAP_SeqAIJ_SeqAIJ for comparison testing */
 51:   MatPtAP(A,aijP,MAT_INITIAL_MATRIX,1.,&C);

 53:   /* Perform PtAP_SeqAIJ_SeqMAIJ */
 54:   /* Developer API */
 55:   MatProductCreate(A,P,NULL,&mC);
 56:   MatProductSetType(mC,MATPRODUCT_PtAP);
 57:   MatProductSetAlgorithm(mC,"default");
 58:   MatProductSetFill(mC,PETSC_DEFAULT);
 59:   MatProductSetFromOptions(mC);
 60:   MatProductSymbolic(mC);
 61:   MatProductNumeric(mC);
 62:   MatProductNumeric(mC);

 64:   /* Check mC = C */
 65:   MatEqual(C,mC,&flg);
 67:   MatDestroy(&mC);

 69:   /* User API */
 70:   MatPtAP(A,P,MAT_INITIAL_MATRIX,1.,&mC);
 71:   MatPtAP(A,P,MAT_REUSE_MATRIX,1.,&mC);

 73:   /* Check mC = C */
 74:   MatEqual(C,mC,&flg);
 76:   MatDestroy(&mC);

 78:   /* Cleanup */
 79:   MatDestroy(&P);
 80:   MatDestroy(&aijP);
 81:   MatDestroy(&A);
 82:   MatDestroy(&C);
 83:   PetscFinalize();
 84:   return 0;
 85: }

 87: /*TEST

 89:    test:
 90:       output_file: output/ex101.out

 92: TEST*/