Actual source code: ex195.c
1: /*
2: * ex195.c
3: *
4: * Created on: Aug 24, 2015
5: * Author: Fande Kong <fdkong.jd@gmail.com>
6: */
8: static char help[] = " Demonstrate the use of MatConvert_Nest_AIJ\n";
10: #include <petscmat.h>
12: int main(int argc, char **args)
13: {
14: Mat A1, A2, A3, A4, A5, B, C, C1, nest;
15: Mat mata[4];
16: Mat aij;
17: MPI_Comm comm;
18: PetscInt m, M, n, istart, iend, ii, i, J, j, K = 10;
19: PetscScalar v;
20: PetscMPIInt size;
21: PetscBool equal;
24: PetscInitialize(&argc, &args, (char *)0, help);
25: comm = PETSC_COMM_WORLD;
26: MPI_Comm_size(comm, &size);
28: /*
29: Assemble the matrix for the five point stencil, YET AGAIN
30: */
31: MatCreate(comm, &A1);
32: m = 2, n = 2;
33: MatSetSizes(A1, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n);
34: MatSetFromOptions(A1);
35: MatSetUp(A1);
36: MatGetOwnershipRange(A1, &istart, &iend);
37: for (ii = istart; ii < iend; ii++) {
38: v = -1.0;
39: i = ii / n;
40: j = ii - i * n;
41: if (i > 0) {
42: J = ii - n;
43: MatSetValues(A1, 1, &ii, 1, &J, &v, INSERT_VALUES);
44: }
45: if (i < m - 1) {
46: J = ii + n;
47: MatSetValues(A1, 1, &ii, 1, &J, &v, INSERT_VALUES);
48: }
49: if (j > 0) {
50: J = ii - 1;
51: MatSetValues(A1, 1, &ii, 1, &J, &v, INSERT_VALUES);
52: }
53: if (j < n - 1) {
54: J = ii + 1;
55: MatSetValues(A1, 1, &ii, 1, &J, &v, INSERT_VALUES);
56: }
57: v = 4.0;
58: MatSetValues(A1, 1, &ii, 1, &ii, &v, INSERT_VALUES);
59: }
60: MatAssemblyBegin(A1, MAT_FINAL_ASSEMBLY);
61: MatAssemblyEnd(A1, MAT_FINAL_ASSEMBLY);
62: MatView(A1, PETSC_VIEWER_STDOUT_WORLD);
64: MatDuplicate(A1, MAT_COPY_VALUES, &A2);
65: MatDuplicate(A1, MAT_COPY_VALUES, &A3);
66: MatDuplicate(A1, MAT_COPY_VALUES, &A4);
68: /*create a nest matrix */
69: MatCreate(comm, &nest);
70: MatSetType(nest, MATNEST);
71: mata[0] = A1, mata[1] = A2, mata[2] = A3, mata[3] = A4;
72: MatNestSetSubMats(nest, 2, NULL, 2, NULL, mata);
73: MatSetUp(nest);
74: MatConvert(nest, MATAIJ, MAT_INITIAL_MATRIX, &aij);
75: MatView(aij, PETSC_VIEWER_STDOUT_WORLD);
77: /* create a dense matrix */
78: MatGetSize(nest, &M, NULL);
79: MatGetLocalSize(nest, &m, NULL);
80: MatCreateDense(comm, m, PETSC_DECIDE, M, K, NULL, &B);
81: MatSetRandom(B, PETSC_NULL);
83: /* C = nest*B_dense */
84: MatMatMult(nest, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C);
85: MatMatMult(nest, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C);
86: MatMatMultEqual(nest, B, C, 10, &equal);
89: /* Test B = nest*C, reuse C and B with MatProductCreateWithMat() */
90: /* C has been obtained from nest*B. Clear internal data structures related to factors to prevent circular references */
91: MatProductClear(C);
92: MatProductCreateWithMat(nest, C, NULL, B);
93: MatProductSetType(B, MATPRODUCT_AB);
94: MatProductSetFromOptions(B);
95: MatProductSymbolic(B);
96: MatProductNumeric(B);
97: MatMatMultEqual(nest, C, B, 10, &equal);
99: MatConvert(nest, MATAIJ, MAT_INPLACE_MATRIX, &nest);
100: MatEqual(nest, aij, &equal);
102: MatDestroy(&nest);
104: if (size > 1) { /* Do not know why this test fails for size = 1 */
105: MatCreateTranspose(A1, &A5); /* A1 is symmetric */
106: mata[0] = A5;
107: MatCreate(comm, &nest);
108: MatSetType(nest, MATNEST);
109: MatNestSetSubMats(nest, 2, NULL, 2, NULL, mata);
110: MatSetUp(nest);
111: MatMatMult(nest, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C1);
112: MatMatMult(nest, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C1);
114: MatMatMultEqual(nest, B, C1, 10, &equal);
116: MatDestroy(&C1);
117: MatDestroy(&A5);
118: MatDestroy(&nest);
119: }
121: MatDestroy(&C);
122: MatDestroy(&B);
123: MatDestroy(&aij);
124: MatDestroy(&A1);
125: MatDestroy(&A2);
126: MatDestroy(&A3);
127: MatDestroy(&A4);
128: PetscFinalize();
129: return 0;
130: }
132: /*TEST
134: test:
135: nsize: 2
137: test:
138: suffix: 2
140: TEST*/