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;
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: MatSetBlockSize(A, bs);
34: MatSeqAIJSetPreallocation(A, PETSC_DEFAULT, NULL);
35: MatMPIAIJSetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL);
36: MatSetFromOptions(A);
37: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
39: /* Create a BAIJ matrix B */
40: MatCreate(PETSC_COMM_WORLD, &B);
41: MatSetSizes(B, m * bs, m * bs, PETSC_DECIDE, PETSC_DECIDE);
42: MatSetType(B, MATBAIJ);
43: MatSeqBAIJSetPreallocation(B, bs, PETSC_DEFAULT, NULL);
44: MatMPIBAIJSetPreallocation(B, bs, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL);
45: MatSetFromOptions(B);
46: MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
48: PetscRandomCreate(PETSC_COMM_WORLD, &rdm);
49: PetscRandomSetFromOptions(rdm);
51: MatGetOwnershipRange(A, &rstart, &rend);
52: MatGetSize(A, &M, &N);
53: Mbs = M / bs;
55: PetscMalloc1(bs, &rows);
56: PetscMalloc1(bs, &cols);
57: PetscMalloc1(bs * bs, &vals);
58: PetscMalloc1(M, &idx);
60: /* Now set blocks of values */
61: for (i = 0; i < 40 * bs; i++) {
62: PetscInt nr = 1, nc = 1;
63: PetscRandomGetValue(rdm, &rval);
64: cols[0] = bs * (int)(PetscRealPart(rval) * Mbs);
65: PetscRandomGetValue(rdm, &rval);
66: rows[0] = rstart + bs * (int)(PetscRealPart(rval) * m);
67: for (j = 1; j < bs; j++) {
68: PetscRandomGetValue(rdm, &rval);
69: if (PetscRealPart(rval) > .5) rows[nr++] = rows[0] + j - 1;
70: }
71: for (j = 1; j < bs; j++) {
72: PetscRandomGetValue(rdm, &rval);
73: if (PetscRealPart(rval) > .5) cols[nc++] = cols[0] + j - 1;
74: }
76: for (j = 0; j < nr * nc; j++) {
77: PetscRandomGetValue(rdm, &rval);
78: vals[j] = rval;
79: }
80: MatSetValues(A, nr, rows, nc, cols, vals, ADD_VALUES);
81: MatSetValues(B, nr, rows, nc, cols, vals, ADD_VALUES);
82: }
83: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
84: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
85: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
86: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
88: /* Test MatConvert_MPIAIJ_MPI(S)BAIJ handles incompletely filled blocks */
89: MatConvert(A, MATBAIJ, MAT_INITIAL_MATRIX, &E);
90: MatDestroy(&E);
91: MatTranspose(A, MAT_INITIAL_MATRIX, &Bt);
92: MatAXPY(Bt, 1.0, B, DIFFERENT_NONZERO_PATTERN);
93: MatSetOption(Bt, MAT_SYMMETRIC, PETSC_TRUE);
94: MatConvert(Bt, MATSBAIJ, MAT_INITIAL_MATRIX, &E);
95: MatDestroy(&E);
96: MatDestroy(&Bt);
98: /* Test MatIncreaseOverlap() */
99: PetscMalloc1(nd, &is1);
100: 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: PetscRandomGetValue(rdm, &rval);
107: sz = (int)(PetscRealPart(rval) * m);
108: for (j = 0; j < sz; j++) {
109: 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: ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is1 + i);
114: ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is2 + i);
115: }
116: MatIncreaseOverlap(A, nd, is1, ov);
117: MatIncreaseOverlap(B, nd, is2, ov);
119: for (i = 0; i < nd; ++i) {
120: ISEqual(is1[i], is2[i], &flg);
122: 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);
123: }
125: for (i = 0; i < nd; ++i) {
126: ISSort(is1[i]);
127: ISSort(is2[i]);
128: }
130: MatCreateSubMatrices(B, nd, is2, is2, MAT_INITIAL_MATRIX, &submatB);
131: MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA);
133: /* Test MatMult() */
134: for (i = 0; i < nd; i++) {
135: MatGetSize(submatA[i], &mm, &nn);
136: VecCreateSeq(PETSC_COMM_SELF, mm, &xx);
137: VecDuplicate(xx, &s1);
138: VecDuplicate(xx, &s2);
139: for (j = 0; j < 3; j++) {
140: VecSetRandom(xx, rdm);
141: MatMult(submatA[i], xx, s1);
142: MatMult(submatB[i], xx, s2);
143: VecNorm(s1, NORM_2, &s1norm);
144: VecNorm(s2, NORM_2, &s2norm);
145: rnorm = s2norm - s1norm;
146: if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", rank, (double)s1norm, (double)s2norm);
147: }
148: VecDestroy(&xx);
149: VecDestroy(&s1);
150: VecDestroy(&s2);
151: }
153: /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
154: MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA);
155: MatCreateSubMatrices(B, nd, is2, is2, MAT_REUSE_MATRIX, &submatB);
157: /* Test MatMult() */
158: for (i = 0; i < nd; i++) {
159: MatGetSize(submatA[i], &mm, &nn);
160: VecCreateSeq(PETSC_COMM_SELF, mm, &xx);
161: VecDuplicate(xx, &s1);
162: VecDuplicate(xx, &s2);
163: for (j = 0; j < 3; j++) {
164: VecSetRandom(xx, rdm);
165: MatMult(submatA[i], xx, s1);
166: MatMult(submatB[i], xx, s2);
167: VecNorm(s1, NORM_2, &s1norm);
168: VecNorm(s2, NORM_2, &s2norm);
169: rnorm = s2norm - s1norm;
170: if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", rank, (double)s1norm, (double)s2norm);
171: }
172: VecDestroy(&xx);
173: VecDestroy(&s1);
174: VecDestroy(&s2);
175: }
177: /* Free allocated memory */
178: for (i = 0; i < nd; ++i) {
179: ISDestroy(&is1[i]);
180: ISDestroy(&is2[i]);
181: }
182: MatDestroySubMatrices(nd, &submatA);
183: MatDestroySubMatrices(nd, &submatB);
185: PetscFree(is1);
186: PetscFree(is2);
187: PetscFree(idx);
188: PetscFree(rows);
189: PetscFree(cols);
190: PetscFree(vals);
191: MatDestroy(&A);
192: MatDestroy(&B);
193: PetscRandomDestroy(&rdm);
194: 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*/