1 static char help[] = "Tests the various sequential routines in MatSBAIJ format. Same as ex74.c except diagonal entries of the matrices are zeros.\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Vec x, y, b, s1, s2;
8 Mat A; /* linear system matrix */
9 Mat sA; /* symmetric part of the matrices */
10 PetscInt n, mbs = 16, bs = 1, nz = 3, prob = 2, i, j, col[3], row, Ii, J, n1;
11 const PetscInt *ip_ptr;
12 PetscScalar neg_one = -1.0, value[3], alpha = 0.1;
13 PetscMPIInt size;
14 IS ip, isrow, iscol;
15 PetscRandom rdm;
16 PetscBool reorder = PETSC_FALSE;
17 MatInfo minfo1, minfo2;
18 PetscReal norm1, norm2, tol = 10 * PETSC_SMALL;
19
20 PetscFunctionBeginUser;
21 PetscCall(PetscInitialize(&argc, &args, NULL, help));
22 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
23 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
24 PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
25 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL));
26
27 n = mbs * bs;
28 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_WORLD, bs, n, n, nz, NULL, &A));
29 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
30 PetscCall(MatCreateSeqSBAIJ(PETSC_COMM_WORLD, bs, n, n, nz, NULL, &sA));
31 PetscCall(MatSetOption(sA, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
32
33 /* Test MatGetOwnershipRange() */
34 PetscCall(MatGetOwnershipRange(A, &Ii, &J));
35 PetscCall(MatGetOwnershipRange(sA, &i, &j));
36 if (i - Ii || j - J) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetOwnershipRange() in MatSBAIJ format\n"));
37
38 /* Assemble matrix */
39 if (bs == 1) {
40 PetscCall(PetscOptionsGetInt(NULL, NULL, "-test_problem", &prob, NULL));
41 if (prob == 1) { /* tridiagonal matrix */
42 value[0] = -1.0;
43 value[1] = 2.0;
44 value[2] = -1.0;
45 for (i = 1; i < n - 1; i++) {
46 col[0] = i - 1;
47 col[1] = i;
48 col[2] = i + 1;
49 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
50 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES));
51 }
52 i = n - 1;
53 col[0] = 0;
54 col[1] = n - 2;
55 col[2] = n - 1;
56
57 value[0] = 0.1;
58 value[1] = -1;
59 value[2] = 2;
60 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
61 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES));
62
63 i = 0;
64 col[0] = 0;
65 col[1] = 1;
66 col[2] = n - 1;
67
68 value[0] = 2.0;
69 value[1] = -1.0;
70 value[2] = 0.1;
71 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
72 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES));
73 } else if (prob == 2) { /* matrix for the five point stencil */
74 n1 = (PetscInt)(PetscSqrtReal((PetscReal)n) + 0.001);
75 PetscCheck(n1 * n1 == n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "sqrt(n) must be a positive integer!");
76 for (i = 0; i < n1; i++) {
77 for (j = 0; j < n1; j++) {
78 Ii = j + n1 * i;
79 if (i > 0) {
80 J = Ii - n1;
81 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
82 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
83 }
84 if (i < n1 - 1) {
85 J = Ii + n1;
86 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
87 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
88 }
89 if (j > 0) {
90 J = Ii - 1;
91 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
92 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
93 }
94 if (j < n1 - 1) {
95 J = Ii + 1;
96 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
97 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
98 }
99 }
100 }
101 }
102 } else { /* bs > 1 */
103 #if defined(DIAGB)
104 for (block = 0; block < n / bs; block++) {
105 /* diagonal blocks */
106 value[0] = -1.0;
107 value[1] = 4.0;
108 value[2] = -1.0;
109 for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) {
110 col[0] = i - 1;
111 col[1] = i;
112 col[2] = i + 1;
113 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
114 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES));
115 }
116 i = bs - 1 + block * bs;
117 col[0] = bs - 2 + block * bs;
118 col[1] = bs - 1 + block * bs;
119
120 value[0] = -1.0;
121 value[1] = 4.0;
122 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
123 PetscCall(MatSetValues(sA, 1, &i, 2, col, value, INSERT_VALUES));
124
125 i = 0 + block * bs;
126 col[0] = 0 + block * bs;
127 col[1] = 1 + block * bs;
128
129 value[0] = 4.0;
130 value[1] = -1.0;
131 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
132 PetscCall(MatSetValues(sA, 1, &i, 2, col, value, INSERT_VALUES));
133 }
134 #endif
135 /* off-diagonal blocks */
136 value[0] = -1.0;
137 for (i = 0; i < (n / bs - 1) * bs; i++) {
138 col[0] = i + bs;
139 PetscCall(MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES));
140 PetscCall(MatSetValues(sA, 1, &i, 1, col, value, INSERT_VALUES));
141 col[0] = i;
142 row = i + bs;
143 PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES));
144 PetscCall(MatSetValues(sA, 1, &row, 1, col, value, INSERT_VALUES));
145 }
146 }
147 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
148 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
149
150 PetscCall(MatAssemblyBegin(sA, MAT_FINAL_ASSEMBLY));
151 PetscCall(MatAssemblyEnd(sA, MAT_FINAL_ASSEMBLY));
152
153 /* Test MatNorm() */
154 PetscCall(MatNorm(A, NORM_FROBENIUS, &norm1));
155 PetscCall(MatNorm(sA, NORM_FROBENIUS, &norm2));
156 norm1 -= norm2;
157 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm(), fnorm1-fnorm2=%16.14e\n", (double)norm1));
158 PetscCall(MatNorm(A, NORM_INFINITY, &norm1));
159 PetscCall(MatNorm(sA, NORM_INFINITY, &norm2));
160 norm1 -= norm2;
161 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm(), inf_norm1-inf_norm2=%16.14e\n", (double)norm1));
162
163 /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
164 PetscCall(MatGetInfo(A, MAT_LOCAL, &minfo1));
165 PetscCall(MatGetInfo(sA, MAT_LOCAL, &minfo2));
166 i = (int)(minfo1.nz_used - minfo2.nz_used);
167 j = (int)(minfo1.nz_allocated - minfo2.nz_allocated);
168 if (i < 0 || j < 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetInfo()\n"));
169
170 PetscCall(MatGetSize(A, &Ii, &J));
171 PetscCall(MatGetSize(sA, &i, &j));
172 if (i - Ii || j - J) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetSize()\n"));
173
174 PetscCall(MatGetBlockSize(A, &Ii));
175 PetscCall(MatGetBlockSize(sA, &i));
176 if (i - Ii) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetBlockSize()\n"));
177
178 /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
179 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
180 PetscCall(PetscRandomSetFromOptions(rdm));
181 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &x));
182 PetscCall(VecDuplicate(x, &s1));
183 PetscCall(VecDuplicate(x, &s2));
184 PetscCall(VecDuplicate(x, &y));
185 PetscCall(VecDuplicate(x, &b));
186
187 PetscCall(VecSetRandom(x, rdm));
188
189 PetscCall(MatDiagonalScale(A, x, x));
190 PetscCall(MatDiagonalScale(sA, x, x));
191
192 PetscCall(MatGetDiagonal(A, s1));
193 PetscCall(MatGetDiagonal(sA, s2));
194 PetscCall(VecNorm(s1, NORM_1, &norm1));
195 PetscCall(VecNorm(s2, NORM_1, &norm2));
196 norm1 -= norm2;
197 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatGetDiagonal() \n"));
198
199 PetscCall(MatScale(A, alpha));
200 PetscCall(MatScale(sA, alpha));
201
202 /* Test MatMult(), MatMultAdd() */
203 for (i = 0; i < 40; i++) {
204 PetscCall(VecSetRandom(x, rdm));
205 PetscCall(MatMult(A, x, s1));
206 PetscCall(MatMult(sA, x, s2));
207 PetscCall(VecNorm(s1, NORM_1, &norm1));
208 PetscCall(VecNorm(s2, NORM_1, &norm2));
209 norm1 -= norm2;
210 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMult(), MatDiagonalScale() or MatScale()\n"));
211 }
212
213 for (i = 0; i < 40; i++) {
214 PetscCall(VecSetRandom(x, rdm));
215 PetscCall(VecSetRandom(y, rdm));
216 PetscCall(MatMultAdd(A, x, y, s1));
217 PetscCall(MatMultAdd(sA, x, y, s2));
218 PetscCall(VecNorm(s1, NORM_1, &norm1));
219 PetscCall(VecNorm(s2, NORM_1, &norm2));
220 norm1 -= norm2;
221 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMultAdd(), MatDiagonalScale() or MatScale() \n"));
222 }
223
224 /* Test MatReordering() */
225 PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &isrow, &iscol));
226 ip = isrow;
227
228 if (reorder) {
229 IS nip;
230 PetscInt *nip_ptr;
231 PetscCall(PetscMalloc1(mbs, &nip_ptr));
232 PetscCall(ISGetIndices(ip, &ip_ptr));
233 PetscCall(PetscArraycpy(nip_ptr, ip_ptr, mbs));
234 i = nip_ptr[1];
235 nip_ptr[1] = nip_ptr[mbs - 2];
236 nip_ptr[mbs - 2] = i;
237 i = nip_ptr[0];
238 nip_ptr[0] = nip_ptr[mbs - 1];
239 nip_ptr[mbs - 1] = i;
240 PetscCall(ISRestoreIndices(ip, &ip_ptr));
241 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, mbs, nip_ptr, PETSC_COPY_VALUES, &nip));
242 PetscCall(PetscFree(nip_ptr));
243
244 PetscCall(MatReorderingSeqSBAIJ(sA, ip));
245 PetscCall(ISDestroy(&nip));
246 }
247
248 PetscCall(ISDestroy(&iscol));
249 PetscCall(ISDestroy(&isrow));
250 PetscCall(MatDestroy(&A));
251 PetscCall(MatDestroy(&sA));
252 PetscCall(VecDestroy(&x));
253 PetscCall(VecDestroy(&y));
254 PetscCall(VecDestroy(&s1));
255 PetscCall(VecDestroy(&s2));
256 PetscCall(VecDestroy(&b));
257 PetscCall(PetscRandomDestroy(&rdm));
258
259 PetscCall(PetscFinalize());
260 return 0;
261 }
262
263 /*TEST
264
265 test:
266 args: -bs {{1 2 3 4 5 6 7 8}}
267 output_file: output/empty.out
268
269 TEST*/
270