1 static char help[] = "Tests various routines in MatMPISBAIJ format.\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Vec x, y, u, s1, s2;
8 Mat A, sA, sB;
9 PetscRandom rctx;
10 PetscReal r1, r2, rnorm, tol = PETSC_SQRT_MACHINE_EPSILON;
11 PetscScalar one = 1.0, neg_one = -1.0, value[3], four = 4.0, alpha = 0.1;
12 PetscInt n, col[3], n1, block, row, i, j, i2, j2, Ii, J, rstart, rend, bs = 1, mbs = 16, d_nz = 3, o_nz = 3, prob = 1;
13 PetscMPIInt size, rank;
14 PetscBool flg;
15
16 PetscFunctionBeginUser;
17 PetscCall(PetscInitialize(&argc, &args, NULL, help));
18 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL));
19 PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
20 PetscCall(PetscOptionsGetInt(NULL, NULL, "-prob", &prob, NULL));
21
22 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
23 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
24
25 /* Create a BAIJ matrix A */
26 n = mbs * bs;
27 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
28 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
29 PetscCall(MatSetType(A, MATBAIJ));
30 PetscCall(MatSetFromOptions(A));
31 PetscCall(MatMPIBAIJSetPreallocation(A, bs, d_nz, NULL, o_nz, NULL));
32 PetscCall(MatSeqBAIJSetPreallocation(A, bs, d_nz, NULL));
33 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
34
35 if (bs == 1) {
36 if (prob == 1) { /* tridiagonal matrix */
37 value[0] = -1.0;
38 value[1] = 0.0;
39 value[2] = -1.0;
40 for (i = 1; i < n - 1; i++) {
41 col[0] = i - 1;
42 col[1] = i;
43 col[2] = i + 1;
44 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
45 }
46 i = n - 1;
47 col[0] = 0;
48 col[1] = n - 2;
49 col[2] = n - 1;
50 value[0] = 0.1;
51 value[1] = -1.0;
52 value[2] = 0.0;
53 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
54
55 i = 0;
56 col[0] = 0;
57 col[1] = 1;
58 col[2] = n - 1;
59 value[0] = 0.0;
60 value[1] = -1.0;
61 value[2] = 0.1;
62 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
63 } else if (prob == 2) { /* matrix for the five point stencil */
64 n1 = (int)PetscSqrtReal((PetscReal)n);
65 for (i = 0; i < n1; i++) {
66 for (j = 0; j < n1; j++) {
67 Ii = j + n1 * i;
68 if (i > 0) {
69 J = Ii - n1;
70 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
71 }
72 if (i < n1 - 1) {
73 J = Ii + n1;
74 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
75 }
76 if (j > 0) {
77 J = Ii - 1;
78 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
79 }
80 if (j < n1 - 1) {
81 J = Ii + 1;
82 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
83 }
84 PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &four, INSERT_VALUES));
85 }
86 }
87 }
88 /* end of if (bs == 1) */
89 } else { /* bs > 1 */
90 for (block = 0; block < n / bs; block++) {
91 /* diagonal blocks */
92 value[0] = -1.0;
93 value[1] = 4.0;
94 value[2] = -1.0;
95 for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) {
96 col[0] = i - 1;
97 col[1] = i;
98 col[2] = i + 1;
99 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
100 }
101 i = bs - 1 + block * bs;
102 col[0] = bs - 2 + block * bs;
103 col[1] = bs - 1 + block * bs;
104 value[0] = -1.0;
105 value[1] = 4.0;
106 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
107
108 i = 0 + block * bs;
109 col[0] = 0 + block * bs;
110 col[1] = 1 + block * bs;
111 value[0] = 4.0;
112 value[1] = -1.0;
113 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
114 }
115 /* off-diagonal blocks */
116 value[0] = -1.0;
117 for (i = 0; i < (n / bs - 1) * bs; i++) {
118 col[0] = i + bs;
119 PetscCall(MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES));
120 col[0] = i;
121 row = i + bs;
122 PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES));
123 }
124 }
125 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
126 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
127 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
128
129 /* Get SBAIJ matrix sA from A */
130 PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &sA));
131
132 /* Test MatGetSize(), MatGetLocalSize() */
133 PetscCall(MatGetSize(sA, &i, &j));
134 PetscCall(MatGetSize(A, &i2, &j2));
135 i -= i2;
136 j -= j2;
137 if (i || j) {
138 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatGetSize()\n", rank));
139 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
140 }
141
142 PetscCall(MatGetLocalSize(sA, &i, &j));
143 PetscCall(MatGetLocalSize(A, &i2, &j2));
144 i2 -= i;
145 j2 -= j;
146 if (i2 || j2) {
147 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatGetLocalSize()\n", rank));
148 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
149 }
150
151 /* vectors */
152 /*--------------------*/
153 /* i is obtained from MatGetLocalSize() */
154 PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
155 PetscCall(VecSetSizes(x, i, PETSC_DECIDE));
156 PetscCall(VecSetFromOptions(x));
157 PetscCall(VecDuplicate(x, &y));
158 PetscCall(VecDuplicate(x, &u));
159 PetscCall(VecDuplicate(x, &s1));
160 PetscCall(VecDuplicate(x, &s2));
161
162 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
163 PetscCall(PetscRandomSetFromOptions(rctx));
164 PetscCall(VecSetRandom(x, rctx));
165 PetscCall(VecSet(u, one));
166
167 /* Test MatNorm() */
168 PetscCall(MatNorm(A, NORM_FROBENIUS, &r1));
169 PetscCall(MatNorm(sA, NORM_FROBENIUS, &r2));
170 rnorm = PetscAbsReal(r1 - r2) / r2;
171 if (rnorm > tol && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs));
172 PetscCall(MatNorm(A, NORM_INFINITY, &r1));
173 PetscCall(MatNorm(sA, NORM_INFINITY, &r2));
174 rnorm = PetscAbsReal(r1 - r2) / r2;
175 if (rnorm > tol && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs));
176 PetscCall(MatNorm(A, NORM_1, &r1));
177 PetscCall(MatNorm(sA, NORM_1, &r2));
178 rnorm = PetscAbsReal(r1 - r2) / r2;
179 if (rnorm > tol && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: MatNorm_1(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs));
180
181 /* Test MatGetOwnershipRange() */
182 PetscCall(MatGetOwnershipRange(sA, &rstart, &rend));
183 PetscCall(MatGetOwnershipRange(A, &i2, &j2));
184 i2 -= rstart;
185 j2 -= rend;
186 if (i2 || j2) {
187 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MaGetOwnershipRange()\n", rank));
188 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
189 }
190
191 /* Test MatDiagonalScale() */
192 PetscCall(MatDiagonalScale(A, x, x));
193 PetscCall(MatDiagonalScale(sA, x, x));
194 PetscCall(MatMultEqual(A, sA, 10, &flg));
195 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "Error in MatDiagonalScale");
196
197 /* Test MatGetDiagonal(), MatScale() */
198 PetscCall(MatGetDiagonal(A, s1));
199 PetscCall(MatGetDiagonal(sA, s2));
200 PetscCall(VecNorm(s1, NORM_1, &r1));
201 PetscCall(VecNorm(s2, NORM_1, &r2));
202 r1 -= r2;
203 if (r1 < -tol || r1 > tol) {
204 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n", rank, (double)r1));
205 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
206 }
207
208 PetscCall(MatScale(A, alpha));
209 PetscCall(MatScale(sA, alpha));
210
211 /* Test MatGetRowMaxAbs() */
212 PetscCall(MatGetRowMaxAbs(A, s1, NULL));
213 PetscCall(MatGetRowMaxAbs(sA, s2, NULL));
214
215 PetscCall(VecNorm(s1, NORM_1, &r1));
216 PetscCall(VecNorm(s2, NORM_1, &r2));
217 r1 -= r2;
218 if (r1 < -tol || r1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetRowMaxAbs() \n"));
219
220 /* Test MatMult(), MatMultAdd() */
221 PetscCall(MatMultEqual(A, sA, 10, &flg));
222 if (!flg) {
223 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMult() or MatScale()\n", rank));
224 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
225 }
226
227 PetscCall(MatMultAddEqual(A, sA, 10, &flg));
228 if (!flg) {
229 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMultAdd()\n", rank));
230 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
231 }
232
233 /* Test MatMultTranspose(), MatMultTransposeAdd() */
234 for (i = 0; i < 10; i++) {
235 PetscCall(VecSetRandom(x, rctx));
236 PetscCall(MatMultTranspose(A, x, s1));
237 PetscCall(MatMultTranspose(sA, x, s2));
238 PetscCall(VecNorm(s1, NORM_1, &r1));
239 PetscCall(VecNorm(s2, NORM_1, &r2));
240 r1 -= r2;
241 if (r1 < -tol || r1 > tol) {
242 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMult() or MatScale(), err=%g\n", rank, (double)r1));
243 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
244 }
245 }
246 for (i = 0; i < 10; i++) {
247 PetscCall(VecSetRandom(x, rctx));
248 PetscCall(VecSetRandom(y, rctx));
249 PetscCall(MatMultTransposeAdd(A, x, y, s1));
250 PetscCall(MatMultTransposeAdd(sA, x, y, s2));
251 PetscCall(VecNorm(s1, NORM_1, &r1));
252 PetscCall(VecNorm(s2, NORM_1, &r2));
253 r1 -= r2;
254 if (r1 < -tol || r1 > tol) {
255 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMultAdd(), err=%g \n", rank, (double)r1));
256 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
257 }
258 }
259
260 /* Test MatDuplicate() */
261 PetscCall(MatDuplicate(sA, MAT_COPY_VALUES, &sB));
262 PetscCall(MatEqual(sA, sB, &flg));
263 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error in MatDuplicate(), sA != sB \n"));
264 PetscCall(MatMultEqual(sA, sB, 5, &flg));
265 if (!flg) {
266 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatDuplicate() or MatMult()\n", rank));
267 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
268 }
269 PetscCall(MatMultAddEqual(sA, sB, 5, &flg));
270 if (!flg) {
271 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatDuplicate() or MatMultAdd(()\n", rank));
272 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
273 }
274 PetscCall(MatDestroy(&sB));
275 PetscCall(VecDestroy(&u));
276 PetscCall(VecDestroy(&x));
277 PetscCall(VecDestroy(&y));
278 PetscCall(VecDestroy(&s1));
279 PetscCall(VecDestroy(&s2));
280 PetscCall(MatDestroy(&sA));
281 PetscCall(MatDestroy(&A));
282 PetscCall(PetscRandomDestroy(&rctx));
283 PetscCall(PetscFinalize());
284 return 0;
285 }
286
287 /*TEST
288
289 test:
290 nsize: {{1 3}}
291 args: -bs {{1 2 3 5 7 8}} -mat_ignore_lower_triangular -prob {{1 2}}
292 output_file: output/empty.out
293
294 TEST*/
295