1 static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for parallel AIJ and BAIJ formats.\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)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;
16
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));
21
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));
27
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));
37
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));
46
47 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
48 PetscCall(PetscRandomSetFromOptions(rdm));
49
50 PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
51 PetscCall(MatGetSize(A, &M, &N));
52 Mbs = M / bs;
53
54 PetscCall(PetscMalloc1(bs, &rows));
55 PetscCall(PetscMalloc1(bs, &cols));
56 PetscCall(PetscMalloc1(bs * bs, &vals));
57 PetscCall(PetscMalloc1(M, &idx));
58
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 }
74
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));
86
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));
96
97 /* Test MatIncreaseOverlap() */
98 PetscCall(PetscMalloc1(nd, &is1));
99 PetscCall(PetscMalloc1(nd, &is2));
100
101 emptynd = PETSC_FALSE;
102 if (rank == 0 && test_nd0) emptynd = PETSC_TRUE; /* test case */
103
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));
117
118 for (i = 0; i < nd; ++i) {
119 PetscCall(ISEqual(is1[i], is2[i], &flg));
120
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 }
123
124 for (i = 0; i < nd; ++i) {
125 PetscCall(ISSort(is1[i]));
126 PetscCall(ISSort(is2[i]));
127 }
128
129 PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_INITIAL_MATRIX, &submatB));
130 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA));
131
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 }
151
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));
155
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 }
175
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));
183
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 }
196
197 /*TEST
198
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/empty.out
203
204 test:
205 suffix: 2
206 args: -nd 2 -test_nd0
207 output_file: output/empty.out
208
209 test:
210 suffix: 3
211 nsize: 3
212 args: -nd 2 -test_nd0
213 output_file: output/empty.out
214
215 TEST*/
216