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*/