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] = {
 11:     {0, 1, 2},
 12:     {3, 4, 5},
 13:     {6, 7, 8}
 14:   };
 15:   Mat         A, mC, C;
 16:   PetscScalar one = 1.;
 17:   PetscMPIInt size;
 18:   PetscBool   flg;

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

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

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

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

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

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

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

 72:   /* User API */
 73:   MatPtAP(A, P, MAT_INITIAL_MATRIX, 1., &mC);
 74:   MatPtAP(A, P, MAT_REUSE_MATRIX, 1., &mC);

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

 81:   /* Cleanup */
 82:   MatDestroy(&P);
 83:   MatDestroy(&aijP);
 84:   MatDestroy(&A);
 85:   MatDestroy(&C);
 86:   PetscFinalize();
 87:   return 0;
 88: }

 90: /*TEST

 92:    test:
 93:       output_file: output/ex101.out

 95: TEST*/