Actual source code: ex54.c
2: static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for parallel AIJ and BAIJ formats.\n";
4: #include <petscmat.h>
6: int main(int argc, char **args)
7: {
8: Mat A, B, *submatA, *submatB;
9: PetscInt bs = 1, m = 11, ov = 1, i, j, k, *rows, *cols, nd = 5, *idx, rstart, rend, sz, mm, nn, M, N, Mbs;
10: PetscMPIInt size, rank;
11: PetscScalar *vals, rval;
12: IS *is1, *is2;
13: PetscRandom rdm;
14: Vec xx, s1, s2;
15: PetscReal s1norm, s2norm, rnorm, tol = 100 * PETSC_SMALL;
16: PetscBool flg, test_nd0 = PETSC_FALSE, emptynd;
19: PetscInitialize(&argc, &args, (char *)0, help);
20: MPI_Comm_size(PETSC_COMM_WORLD, &size);
21: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
23: PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL);
24: PetscOptionsGetInt(NULL, NULL, "-mat_size", &m, NULL);
25: PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL);
26: PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL);
27: PetscOptionsGetBool(NULL, NULL, "-test_nd0", &test_nd0, NULL);
29: /* Create a AIJ matrix A */
30: MatCreate(PETSC_COMM_WORLD, &A);
31: MatSetSizes(A, m * bs, m * bs, PETSC_DECIDE, PETSC_DECIDE);
32: MatSetType(A, MATAIJ);
33: MatSeqAIJSetPreallocation(A, PETSC_DEFAULT, NULL);
34: MatMPIAIJSetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL);
35: MatSetFromOptions(A);
36: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
38: /* Create a BAIJ matrix B */
39: MatCreate(PETSC_COMM_WORLD, &B);
40: MatSetSizes(B, m * bs, m * bs, PETSC_DECIDE, PETSC_DECIDE);
41: MatSetType(B, MATBAIJ);
42: MatSeqBAIJSetPreallocation(B, bs, PETSC_DEFAULT, NULL);
43: MatMPIBAIJSetPreallocation(B, bs, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL);
44: MatSetFromOptions(B);
45: MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
47: PetscRandomCreate(PETSC_COMM_WORLD, &rdm);
48: PetscRandomSetFromOptions(rdm);
50: MatGetOwnershipRange(A, &rstart, &rend);
51: MatGetSize(A, &M, &N);
52: Mbs = M / bs;
54: PetscMalloc1(bs, &rows);
55: PetscMalloc1(bs, &cols);
56: PetscMalloc1(bs * bs, &vals);
57: PetscMalloc1(M, &idx);
59: /* Now set blocks of values */
60: for (i = 0; i < 40 * bs; i++) {
61: PetscRandomGetValue(rdm, &rval);
62: cols[0] = bs * (int)(PetscRealPart(rval) * Mbs);
63: PetscRandomGetValue(rdm, &rval);
64: rows[0] = rstart + bs * (int)(PetscRealPart(rval) * m);
65: for (j = 1; j < bs; j++) {
66: rows[j] = rows[j - 1] + 1;
67: cols[j] = cols[j - 1] + 1;
68: }
70: for (j = 0; j < bs * bs; j++) {
71: PetscRandomGetValue(rdm, &rval);
72: vals[j] = rval;
73: }
74: MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES);
75: MatSetValues(B, bs, rows, bs, cols, vals, ADD_VALUES);
76: }
77: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
78: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
79: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
80: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
82: /* Test MatIncreaseOverlap() */
83: PetscMalloc1(nd, &is1);
84: PetscMalloc1(nd, &is2);
86: emptynd = PETSC_FALSE;
87: if (rank == 0 && test_nd0) emptynd = PETSC_TRUE; /* test case */
89: for (i = 0; i < nd; i++) {
90: PetscRandomGetValue(rdm, &rval);
91: sz = (int)(PetscRealPart(rval) * m);
92: for (j = 0; j < sz; j++) {
93: PetscRandomGetValue(rdm, &rval);
94: idx[j * bs] = bs * (int)(PetscRealPart(rval) * Mbs);
95: for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k;
96: }
97: ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is1 + i);
98: ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is2 + i);
99: }
100: MatIncreaseOverlap(A, nd, is1, ov);
101: MatIncreaseOverlap(B, nd, is2, ov);
103: for (i = 0; i < nd; ++i) {
104: ISEqual(is1[i], is2[i], &flg);
106: if (!flg) PetscPrintf(PETSC_COMM_SELF, "i=%" PetscInt_FMT ", flg=%d :bs=%" PetscInt_FMT " m=%" PetscInt_FMT " ov=%" PetscInt_FMT " nd=%" PetscInt_FMT " np=%d\n", i, flg, bs, m, ov, nd, size);
107: }
109: for (i = 0; i < nd; ++i) {
110: ISSort(is1[i]);
111: ISSort(is2[i]);
112: }
114: MatCreateSubMatrices(B, nd, is2, is2, MAT_INITIAL_MATRIX, &submatB);
115: MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA);
117: /* Test MatMult() */
118: for (i = 0; i < nd; i++) {
119: MatGetSize(submatA[i], &mm, &nn);
120: VecCreateSeq(PETSC_COMM_SELF, mm, &xx);
121: VecDuplicate(xx, &s1);
122: VecDuplicate(xx, &s2);
123: for (j = 0; j < 3; j++) {
124: VecSetRandom(xx, rdm);
125: MatMult(submatA[i], xx, s1);
126: MatMult(submatB[i], xx, s2);
127: VecNorm(s1, NORM_2, &s1norm);
128: VecNorm(s2, NORM_2, &s2norm);
129: rnorm = s2norm - s1norm;
130: if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", rank, (double)s1norm, (double)s2norm);
131: }
132: VecDestroy(&xx);
133: VecDestroy(&s1);
134: VecDestroy(&s2);
135: }
137: /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
138: MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA);
139: MatCreateSubMatrices(B, nd, is2, is2, MAT_REUSE_MATRIX, &submatB);
141: /* Test MatMult() */
142: for (i = 0; i < nd; i++) {
143: MatGetSize(submatA[i], &mm, &nn);
144: VecCreateSeq(PETSC_COMM_SELF, mm, &xx);
145: VecDuplicate(xx, &s1);
146: VecDuplicate(xx, &s2);
147: for (j = 0; j < 3; j++) {
148: VecSetRandom(xx, rdm);
149: MatMult(submatA[i], xx, s1);
150: MatMult(submatB[i], xx, s2);
151: VecNorm(s1, NORM_2, &s1norm);
152: VecNorm(s2, NORM_2, &s2norm);
153: rnorm = s2norm - s1norm;
154: if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", rank, (double)s1norm, (double)s2norm);
155: }
156: VecDestroy(&xx);
157: VecDestroy(&s1);
158: VecDestroy(&s2);
159: }
161: /* Free allocated memory */
162: for (i = 0; i < nd; ++i) {
163: ISDestroy(&is1[i]);
164: ISDestroy(&is2[i]);
165: }
166: MatDestroySubMatrices(nd, &submatA);
167: MatDestroySubMatrices(nd, &submatB);
169: PetscFree(is1);
170: PetscFree(is2);
171: PetscFree(idx);
172: PetscFree(rows);
173: PetscFree(cols);
174: PetscFree(vals);
175: MatDestroy(&A);
176: MatDestroy(&B);
177: PetscRandomDestroy(&rdm);
178: PetscFinalize();
179: return 0;
180: }
182: /*TEST
184: test:
185: nsize: {{1 3}}
186: args: -mat_block_size {{1 3 4 6 8}} -ov {{1 3}} -mat_size {{11 13}} -nd 7
187: output_file: output/ex54.out
189: test:
190: suffix: 2
191: args: -nd 2 -test_nd0
192: output_file: output/ex54.out
194: test:
195: suffix: 3
196: nsize: 3
197: args: -nd 2 -test_nd0
198: output_file: output/ex54.out
200: TEST*/