Actual source code: ex54.c
1: static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for parallel AIJ and BAIJ formats.\n";
3: #include <petscmat.h>
5: int main(int argc, char **args)
6: {
7: Mat E, A, B, Bt, *submatA, *submatB;
8: PetscInt bs = 1, m = 11, ov = 1, i, j, k, *rows, *cols, nd = 5, *idx, rstart, rend, sz, mm, nn, M, N, Mbs;
9: PetscMPIInt size, rank;
10: PetscScalar *vals, rval;
11: IS *is1, *is2;
12: PetscRandom rdm;
13: Vec xx, s1, s2;
14: PetscReal s1norm, s2norm, rnorm, tol = 100 * PETSC_SMALL;
15: PetscBool flg, test_nd0 = PETSC_FALSE, emptynd;
17: PetscFunctionBeginUser;
18: PetscCall(PetscInitialize(&argc, &args, NULL, help));
19: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
20: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
22: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL));
23: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_size", &m, NULL));
24: PetscCall(PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL));
25: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL));
26: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nd0", &test_nd0, NULL));
28: /* Create a AIJ matrix A */
29: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
30: PetscCall(MatSetSizes(A, m * bs, m * bs, PETSC_DECIDE, PETSC_DECIDE));
31: PetscCall(MatSetType(A, MATAIJ));
32: PetscCall(MatSetBlockSize(A, bs));
33: PetscCall(MatSeqAIJSetPreallocation(A, PETSC_DEFAULT, NULL));
34: PetscCall(MatMPIAIJSetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));
35: PetscCall(MatSetFromOptions(A));
36: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
38: /* Create a BAIJ matrix B */
39: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
40: PetscCall(MatSetSizes(B, m * bs, m * bs, PETSC_DECIDE, PETSC_DECIDE));
41: PetscCall(MatSetType(B, MATBAIJ));
42: PetscCall(MatSeqBAIJSetPreallocation(B, bs, PETSC_DEFAULT, NULL));
43: PetscCall(MatMPIBAIJSetPreallocation(B, bs, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));
44: PetscCall(MatSetFromOptions(B));
45: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
47: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
48: PetscCall(PetscRandomSetFromOptions(rdm));
50: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
51: PetscCall(MatGetSize(A, &M, &N));
52: Mbs = M / bs;
54: PetscCall(PetscMalloc1(bs, &rows));
55: PetscCall(PetscMalloc1(bs, &cols));
56: PetscCall(PetscMalloc1(bs * bs, &vals));
57: PetscCall(PetscMalloc1(M, &idx));
59: /* Now set blocks of values */
60: for (i = 0; i < 40 * bs; i++) {
61: PetscInt nr = 1, nc = 1;
62: PetscCall(PetscRandomGetValue(rdm, &rval));
63: cols[0] = bs * (int)(PetscRealPart(rval) * Mbs);
64: PetscCall(PetscRandomGetValue(rdm, &rval));
65: rows[0] = rstart + bs * (int)(PetscRealPart(rval) * m);
66: for (j = 1; j < bs; j++) {
67: PetscCall(PetscRandomGetValue(rdm, &rval));
68: if (PetscRealPart(rval) > .5) rows[nr++] = rows[0] + j - 1;
69: }
70: for (j = 1; j < bs; j++) {
71: PetscCall(PetscRandomGetValue(rdm, &rval));
72: if (PetscRealPart(rval) > .5) cols[nc++] = cols[0] + j - 1;
73: }
75: for (j = 0; j < nr * nc; j++) {
76: PetscCall(PetscRandomGetValue(rdm, &rval));
77: vals[j] = rval;
78: }
79: PetscCall(MatSetValues(A, nr, rows, nc, cols, vals, ADD_VALUES));
80: PetscCall(MatSetValues(B, nr, rows, nc, cols, vals, ADD_VALUES));
81: }
82: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
83: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
84: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
85: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
87: /* Test MatConvert_MPIAIJ_MPI(S)BAIJ handles incompletely filled blocks */
88: PetscCall(MatConvert(A, MATBAIJ, MAT_INITIAL_MATRIX, &E));
89: PetscCall(MatDestroy(&E));
90: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Bt));
91: PetscCall(MatAXPY(Bt, 1.0, B, DIFFERENT_NONZERO_PATTERN));
92: PetscCall(MatSetOption(Bt, MAT_SYMMETRIC, PETSC_TRUE));
93: PetscCall(MatConvert(Bt, MATSBAIJ, MAT_INITIAL_MATRIX, &E));
94: PetscCall(MatDestroy(&E));
95: PetscCall(MatDestroy(&Bt));
97: /* Test MatIncreaseOverlap() */
98: PetscCall(PetscMalloc1(nd, &is1));
99: PetscCall(PetscMalloc1(nd, &is2));
101: emptynd = PETSC_FALSE;
102: if (rank == 0 && test_nd0) emptynd = PETSC_TRUE; /* test case */
104: for (i = 0; i < nd; i++) {
105: PetscCall(PetscRandomGetValue(rdm, &rval));
106: sz = (int)(PetscRealPart(rval) * m);
107: for (j = 0; j < sz; j++) {
108: PetscCall(PetscRandomGetValue(rdm, &rval));
109: idx[j * bs] = bs * (int)(PetscRealPart(rval) * Mbs);
110: for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k;
111: }
112: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is1 + i));
113: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is2 + i));
114: }
115: PetscCall(MatIncreaseOverlap(A, nd, is1, ov));
116: PetscCall(MatIncreaseOverlap(B, nd, is2, ov));
118: for (i = 0; i < nd; ++i) {
119: PetscCall(ISEqual(is1[i], is2[i], &flg));
121: 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));
122: }
124: for (i = 0; i < nd; ++i) {
125: PetscCall(ISSort(is1[i]));
126: PetscCall(ISSort(is2[i]));
127: }
129: PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_INITIAL_MATRIX, &submatB));
130: PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA));
132: /* Test MatMult() */
133: for (i = 0; i < nd; i++) {
134: PetscCall(MatGetSize(submatA[i], &mm, &nn));
135: PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx));
136: PetscCall(VecDuplicate(xx, &s1));
137: PetscCall(VecDuplicate(xx, &s2));
138: for (j = 0; j < 3; j++) {
139: PetscCall(VecSetRandom(xx, rdm));
140: PetscCall(MatMult(submatA[i], xx, s1));
141: PetscCall(MatMult(submatB[i], xx, s2));
142: PetscCall(VecNorm(s1, NORM_2, &s1norm));
143: PetscCall(VecNorm(s2, NORM_2, &s2norm));
144: rnorm = s2norm - s1norm;
145: 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));
146: }
147: PetscCall(VecDestroy(&xx));
148: PetscCall(VecDestroy(&s1));
149: PetscCall(VecDestroy(&s2));
150: }
152: /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
153: PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA));
154: PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_REUSE_MATRIX, &submatB));
156: /* Test MatMult() */
157: for (i = 0; i < nd; i++) {
158: PetscCall(MatGetSize(submatA[i], &mm, &nn));
159: PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx));
160: PetscCall(VecDuplicate(xx, &s1));
161: PetscCall(VecDuplicate(xx, &s2));
162: for (j = 0; j < 3; j++) {
163: PetscCall(VecSetRandom(xx, rdm));
164: PetscCall(MatMult(submatA[i], xx, s1));
165: PetscCall(MatMult(submatB[i], xx, s2));
166: PetscCall(VecNorm(s1, NORM_2, &s1norm));
167: PetscCall(VecNorm(s2, NORM_2, &s2norm));
168: rnorm = s2norm - s1norm;
169: 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));
170: }
171: PetscCall(VecDestroy(&xx));
172: PetscCall(VecDestroy(&s1));
173: PetscCall(VecDestroy(&s2));
174: }
176: /* Free allocated memory */
177: for (i = 0; i < nd; ++i) {
178: PetscCall(ISDestroy(&is1[i]));
179: PetscCall(ISDestroy(&is2[i]));
180: }
181: PetscCall(MatDestroySubMatrices(nd, &submatA));
182: PetscCall(MatDestroySubMatrices(nd, &submatB));
184: PetscCall(PetscFree(is1));
185: PetscCall(PetscFree(is2));
186: PetscCall(PetscFree(idx));
187: PetscCall(PetscFree(rows));
188: PetscCall(PetscFree(cols));
189: PetscCall(PetscFree(vals));
190: PetscCall(MatDestroy(&A));
191: PetscCall(MatDestroy(&B));
192: PetscCall(PetscRandomDestroy(&rdm));
193: PetscCall(PetscFinalize());
194: return 0;
195: }
197: /*TEST
199: test:
200: nsize: {{1 3}}
201: args: -mat_block_size {{1 3 4 6 8}} -ov {{1 3}} -mat_size {{11 13}} -nd 7
202: output_file: output/ex54.out
204: test:
205: suffix: 2
206: args: -nd 2 -test_nd0
207: output_file: output/ex54.out
209: test:
210: suffix: 3
211: nsize: 3
212: args: -nd 2 -test_nd0
213: output_file: output/ex54.out
215: TEST*/