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 E, A, B, Bt, *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;
18: PetscFunctionBeginUser;
19: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
20: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
21: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
23: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL));
24: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_size", &m, NULL));
25: PetscCall(PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL));
26: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL));
27: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nd0", &test_nd0, NULL));
29: /* Create a AIJ matrix A */
30: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
31: PetscCall(MatSetSizes(A, m * bs, m * bs, PETSC_DECIDE, PETSC_DECIDE));
32: PetscCall(MatSetType(A, MATAIJ));
33: PetscCall(MatSetBlockSize(A, bs));
34: PetscCall(MatSeqAIJSetPreallocation(A, PETSC_DEFAULT, NULL));
35: PetscCall(MatMPIAIJSetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));
36: PetscCall(MatSetFromOptions(A));
37: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
39: /* Create a BAIJ matrix B */
40: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
41: PetscCall(MatSetSizes(B, m * bs, m * bs, PETSC_DECIDE, PETSC_DECIDE));
42: PetscCall(MatSetType(B, MATBAIJ));
43: PetscCall(MatSeqBAIJSetPreallocation(B, bs, PETSC_DEFAULT, NULL));
44: PetscCall(MatMPIBAIJSetPreallocation(B, bs, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));
45: PetscCall(MatSetFromOptions(B));
46: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
48: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
49: PetscCall(PetscRandomSetFromOptions(rdm));
51: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
52: PetscCall(MatGetSize(A, &M, &N));
53: Mbs = M / bs;
55: PetscCall(PetscMalloc1(bs, &rows));
56: PetscCall(PetscMalloc1(bs, &cols));
57: PetscCall(PetscMalloc1(bs * bs, &vals));
58: PetscCall(PetscMalloc1(M, &idx));
60: /* Now set blocks of values */
61: for (i = 0; i < 40 * bs; i++) {
62: PetscInt nr = 1, nc = 1;
63: PetscCall(PetscRandomGetValue(rdm, &rval));
64: cols[0] = bs * (int)(PetscRealPart(rval) * Mbs);
65: PetscCall(PetscRandomGetValue(rdm, &rval));
66: rows[0] = rstart + bs * (int)(PetscRealPart(rval) * m);
67: for (j = 1; j < bs; j++) {
68: PetscCall(PetscRandomGetValue(rdm, &rval));
69: if (PetscRealPart(rval) > .5) rows[nr++] = rows[0] + j - 1;
70: }
71: for (j = 1; j < bs; j++) {
72: PetscCall(PetscRandomGetValue(rdm, &rval));
73: if (PetscRealPart(rval) > .5) cols[nc++] = cols[0] + j - 1;
74: }
76: for (j = 0; j < nr * nc; j++) {
77: PetscCall(PetscRandomGetValue(rdm, &rval));
78: vals[j] = rval;
79: }
80: PetscCall(MatSetValues(A, nr, rows, nc, cols, vals, ADD_VALUES));
81: PetscCall(MatSetValues(B, nr, rows, nc, cols, vals, ADD_VALUES));
82: }
83: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
84: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
85: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
86: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
88: /* Test MatConvert_MPIAIJ_MPI(S)BAIJ handles incompletely filled blocks */
89: PetscCall(MatConvert(A, MATBAIJ, MAT_INITIAL_MATRIX, &E));
90: PetscCall(MatDestroy(&E));
91: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Bt));
92: PetscCall(MatAXPY(Bt, 1.0, B, DIFFERENT_NONZERO_PATTERN));
93: PetscCall(MatSetOption(Bt, MAT_SYMMETRIC, PETSC_TRUE));
94: PetscCall(MatConvert(Bt, MATSBAIJ, MAT_INITIAL_MATRIX, &E));
95: PetscCall(MatDestroy(&E));
96: PetscCall(MatDestroy(&Bt));
98: /* Test MatIncreaseOverlap() */
99: PetscCall(PetscMalloc1(nd, &is1));
100: PetscCall(PetscMalloc1(nd, &is2));
102: emptynd = PETSC_FALSE;
103: if (rank == 0 && test_nd0) emptynd = PETSC_TRUE; /* test case */
105: for (i = 0; i < nd; i++) {
106: PetscCall(PetscRandomGetValue(rdm, &rval));
107: sz = (int)(PetscRealPart(rval) * m);
108: for (j = 0; j < sz; j++) {
109: PetscCall(PetscRandomGetValue(rdm, &rval));
110: idx[j * bs] = bs * (int)(PetscRealPart(rval) * Mbs);
111: for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k;
112: }
113: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is1 + i));
114: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is2 + i));
115: }
116: PetscCall(MatIncreaseOverlap(A, nd, is1, ov));
117: PetscCall(MatIncreaseOverlap(B, nd, is2, ov));
119: for (i = 0; i < nd; ++i) {
120: PetscCall(ISEqual(is1[i], is2[i], &flg));
122: if (!flg) PetscCall(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));
123: }
125: for (i = 0; i < nd; ++i) {
126: PetscCall(ISSort(is1[i]));
127: PetscCall(ISSort(is2[i]));
128: }
130: PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_INITIAL_MATRIX, &submatB));
131: PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA));
133: /* Test MatMult() */
134: for (i = 0; i < nd; i++) {
135: PetscCall(MatGetSize(submatA[i], &mm, &nn));
136: PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx));
137: PetscCall(VecDuplicate(xx, &s1));
138: PetscCall(VecDuplicate(xx, &s2));
139: for (j = 0; j < 3; j++) {
140: PetscCall(VecSetRandom(xx, rdm));
141: PetscCall(MatMult(submatA[i], xx, s1));
142: PetscCall(MatMult(submatB[i], xx, s2));
143: PetscCall(VecNorm(s1, NORM_2, &s1norm));
144: PetscCall(VecNorm(s2, NORM_2, &s2norm));
145: rnorm = s2norm - s1norm;
146: if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", rank, (double)s1norm, (double)s2norm));
147: }
148: PetscCall(VecDestroy(&xx));
149: PetscCall(VecDestroy(&s1));
150: PetscCall(VecDestroy(&s2));
151: }
153: /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
154: PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA));
155: PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_REUSE_MATRIX, &submatB));
157: /* Test MatMult() */
158: for (i = 0; i < nd; i++) {
159: PetscCall(MatGetSize(submatA[i], &mm, &nn));
160: PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx));
161: PetscCall(VecDuplicate(xx, &s1));
162: PetscCall(VecDuplicate(xx, &s2));
163: for (j = 0; j < 3; j++) {
164: PetscCall(VecSetRandom(xx, rdm));
165: PetscCall(MatMult(submatA[i], xx, s1));
166: PetscCall(MatMult(submatB[i], xx, s2));
167: PetscCall(VecNorm(s1, NORM_2, &s1norm));
168: PetscCall(VecNorm(s2, NORM_2, &s2norm));
169: rnorm = s2norm - s1norm;
170: if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", rank, (double)s1norm, (double)s2norm));
171: }
172: PetscCall(VecDestroy(&xx));
173: PetscCall(VecDestroy(&s1));
174: PetscCall(VecDestroy(&s2));
175: }
177: /* Free allocated memory */
178: for (i = 0; i < nd; ++i) {
179: PetscCall(ISDestroy(&is1[i]));
180: PetscCall(ISDestroy(&is2[i]));
181: }
182: PetscCall(MatDestroySubMatrices(nd, &submatA));
183: PetscCall(MatDestroySubMatrices(nd, &submatB));
185: PetscCall(PetscFree(is1));
186: PetscCall(PetscFree(is2));
187: PetscCall(PetscFree(idx));
188: PetscCall(PetscFree(rows));
189: PetscCall(PetscFree(cols));
190: PetscCall(PetscFree(vals));
191: PetscCall(MatDestroy(&A));
192: PetscCall(MatDestroy(&B));
193: PetscCall(PetscRandomDestroy(&rdm));
194: PetscCall(PetscFinalize());
195: return 0;
196: }
198: /*TEST
200: test:
201: nsize: {{1 3}}
202: args: -mat_block_size {{1 3 4 6 8}} -ov {{1 3}} -mat_size {{11 13}} -nd 7
203: output_file: output/ex54.out
205: test:
206: suffix: 2
207: args: -nd 2 -test_nd0
208: output_file: output/ex54.out
210: test:
211: suffix: 3
212: nsize: 3
213: args: -nd 2 -test_nd0
214: output_file: output/ex54.out
216: TEST*/