Actual source code: ex109.c
1: static char help[] = "Test MatMatMult() for AIJ and Dense matrices.\n\n";
3: #include <petscmat.h>
5: int main(int argc, char **argv)
6: {
7: Mat A, B, C, D, AT;
8: PetscInt i, M, N, Istart, Iend, n = 7, j, J, Ii, m = 8, am, an;
9: PetscScalar v;
10: PetscRandom r;
11: PetscBool equal = PETSC_FALSE, flg;
12: PetscReal fill = 1.0, norm;
13: PetscMPIInt size;
16: PetscInitialize(&argc, &argv, (char *)0, help);
17: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
18: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
19: PetscOptionsGetReal(NULL, NULL, "-fill", &fill, NULL);
21: PetscRandomCreate(PETSC_COMM_WORLD, &r);
22: PetscRandomSetFromOptions(r);
24: /* Create a aij matrix A */
25: M = N = m * n;
26: MatCreate(PETSC_COMM_WORLD, &A);
27: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N);
28: MatSetType(A, MATAIJ);
29: MatSetFromOptions(A);
30: MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL);
31: MatSeqAIJSetPreallocation(A, 5, NULL);
33: MatGetOwnershipRange(A, &Istart, &Iend);
34: am = Iend - Istart;
35: for (Ii = Istart; Ii < Iend; Ii++) {
36: v = -1.0;
37: i = Ii / n;
38: j = Ii - i * n;
39: if (i > 0) {
40: J = Ii - n;
41: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
42: }
43: if (i < m - 1) {
44: J = Ii + n;
45: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
46: }
47: if (j > 0) {
48: J = Ii - 1;
49: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
50: }
51: if (j < n - 1) {
52: J = Ii + 1;
53: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
54: }
55: v = 4.0;
56: MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
57: }
58: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
59: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
61: /* Create a dense matrix B */
62: MatGetLocalSize(A, &am, &an);
63: MatCreate(PETSC_COMM_WORLD, &B);
64: MatSetSizes(B, an, PETSC_DECIDE, PETSC_DECIDE, M);
65: MatSetType(B, MATDENSE);
66: MatSeqDenseSetPreallocation(B, NULL);
67: MatMPIDenseSetPreallocation(B, NULL);
68: MatSetFromOptions(B);
69: MatSetRandom(B, r);
70: PetscRandomDestroy(&r);
72: /* Test reuse of user-provided dense C (unassembled) -- not recommended usage */
73: MatCreate(PETSC_COMM_WORLD, &C);
74: MatSetType(C, MATDENSE);
75: MatSetSizes(C, am, PETSC_DECIDE, PETSC_DECIDE, N);
76: MatSetFromOptions(C);
77: MatSetUp(C);
78: MatZeroEntries(C);
79: MatMatMult(A, B, MAT_REUSE_MATRIX, fill, &C);
80: MatNorm(C, NORM_INFINITY, &norm);
81: MatDestroy(&C);
83: /* Test C = A*B (aij*dense) */
84: MatMatMult(A, B, MAT_INITIAL_MATRIX, fill, &C);
85: MatMatMult(A, B, MAT_REUSE_MATRIX, fill, &C);
87: /* Test developer API */
88: MatProductCreate(A, B, NULL, &D);
89: MatProductSetType(D, MATPRODUCT_AB);
90: MatProductSetAlgorithm(D, "default");
91: MatProductSetFill(D, fill);
92: MatProductSetFromOptions(D);
93: MatProductSymbolic(D);
94: for (i = 0; i < 2; i++) MatProductNumeric(D);
95: MatEqual(C, D, &equal);
97: MatDestroy(&D);
99: /* Test D = AT*B (transpose(aij)*dense) */
100: MatCreateTranspose(A, &AT);
101: MatMatMult(AT, B, MAT_INITIAL_MATRIX, fill, &D);
102: MatMatMultEqual(AT, B, D, 10, &equal);
104: MatDestroy(&D);
105: MatDestroy(&AT);
107: /* Test D = C*A (dense*aij) */
108: MatMatMult(C, A, MAT_INITIAL_MATRIX, fill, &D);
109: MatMatMult(C, A, MAT_REUSE_MATRIX, fill, &D);
110: MatMatMultEqual(C, A, D, 10, &equal);
112: MatDestroy(&D);
114: /* Test D = A*C (aij*dense) */
115: MatMatMult(A, C, MAT_INITIAL_MATRIX, fill, &D);
116: MatMatMult(A, C, MAT_REUSE_MATRIX, fill, &D);
117: MatMatMultEqual(A, C, D, 10, &equal);
119: MatDestroy(&D);
121: /* Test D = B*C (dense*dense) */
122: MPI_Comm_size(PETSC_COMM_WORLD, &size);
123: if (size == 1) {
124: MatMatMult(B, C, MAT_INITIAL_MATRIX, fill, &D);
125: MatMatMult(B, C, MAT_REUSE_MATRIX, fill, &D);
126: MatMatMultEqual(B, C, D, 10, &equal);
128: MatDestroy(&D);
129: }
131: /* Test D = B*C^T (dense*dense) */
132: MatMatTransposeMult(B, C, MAT_INITIAL_MATRIX, fill, &D);
133: MatMatTransposeMult(B, C, MAT_REUSE_MATRIX, fill, &D);
134: MatMatTransposeMultEqual(B, C, D, 10, &equal);
136: MatDestroy(&D);
138: /* Test MatProductCreateWithMat() and reuse C and B for B = A*C */
139: flg = PETSC_FALSE;
140: PetscOptionsHasName(NULL, NULL, "-test_userAPI", &flg);
141: if (flg) {
142: /* user driver */
143: MatMatMult(A, C, MAT_REUSE_MATRIX, fill, &B);
144: } else {
145: /* clear internal data structures related with previous products to avoid circular references */
146: MatProductClear(A);
147: MatProductClear(B);
148: MatProductClear(C);
149: MatProductCreateWithMat(A, C, NULL, B);
150: MatProductSetType(B, MATPRODUCT_AB);
151: MatProductSetFromOptions(B);
152: MatProductSymbolic(B);
153: MatProductNumeric(B);
154: }
156: MatDestroy(&C);
157: MatDestroy(&B);
158: MatDestroy(&A);
159: PetscFinalize();
160: return 0;
161: }
163: /*TEST
165: test:
166: args: -M 10 -N 10
167: output_file: output/ex109.out
169: test:
170: suffix: 2
171: nsize: 3
172: output_file: output/ex109.out
174: test:
175: suffix: 3
176: nsize: 2
177: args: -matmattransmult_mpidense_mpidense_via cyclic
178: output_file: output/ex109.out
180: test:
181: suffix: 4
182: args: -test_userAPI
183: output_file: output/ex109.out
185: test:
186: suffix: 5
187: nsize: 3
188: args: -test_userAPI
189: output_file: output/ex109.out
191: TEST*/