1 static char help[] = "Example of using MatPreallocator\n\n";
2
3 #include <petscmat.h>
4
ex1_nonsquare_bs1(void)5 PetscErrorCode ex1_nonsquare_bs1(void)
6 {
7 Mat A, preallocator;
8 PetscInt M, N, m, n, bs;
9
10 /*
11 Create the Jacobian matrix
12 */
13 PetscFunctionBegin;
14 M = 10;
15 N = 8;
16 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
17 PetscCall(MatSetType(A, MATAIJ));
18 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N));
19 PetscCall(MatSetBlockSize(A, 1));
20 PetscCall(MatSetFromOptions(A));
21
22 /*
23 Get the sizes of the jacobian.
24 */
25 PetscCall(MatGetLocalSize(A, &m, &n));
26 PetscCall(MatGetBlockSize(A, &bs));
27
28 /*
29 Create a preallocator matrix with sizes (local and global) matching the jacobian A.
30 */
31 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &preallocator));
32 PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
33 PetscCall(MatSetSizes(preallocator, m, n, M, N));
34 PetscCall(MatSetBlockSize(preallocator, bs));
35 PetscCall(MatSetUp(preallocator));
36
37 /*
38 Insert non-zero pattern (e.g. perform a sweep over the grid).
39 You can use MatSetValues(), MatSetValuesBlocked() or MatSetValue().
40 */
41 {
42 PetscInt ii, jj;
43 PetscScalar vv = 0.0;
44
45 ii = 3;
46 jj = 3;
47 PetscCall(MatSetValues(preallocator, 1, &ii, 1, &jj, &vv, INSERT_VALUES));
48
49 ii = 7;
50 jj = 4;
51 PetscCall(MatSetValues(preallocator, 1, &ii, 1, &jj, &vv, INSERT_VALUES));
52
53 ii = 9;
54 jj = 7;
55 PetscCall(MatSetValue(preallocator, ii, jj, vv, INSERT_VALUES));
56 }
57 PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
58 PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
59
60 /*
61 Push the non-zero pattern defined within preallocator into A.
62 Internally this will call MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE).
63 The arg fill = PETSC_TRUE will insert zeros in the matrix A automatically.
64 */
65 PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_TRUE, A));
66
67 /*
68 We no longer require the preallocator object so destroy it.
69 */
70 PetscCall(MatDestroy(&preallocator));
71
72 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
73
74 /*
75 Insert non-zero values into A.
76 */
77 {
78 PetscInt ii, jj;
79 PetscScalar vv;
80
81 ii = 3;
82 jj = 3;
83 vv = 0.3;
84 PetscCall(MatSetValue(A, ii, jj, vv, INSERT_VALUES));
85
86 ii = 7;
87 jj = 4;
88 vv = 3.3;
89 PetscCall(MatSetValues(A, 1, &ii, 1, &jj, &vv, INSERT_VALUES));
90
91 ii = 9;
92 jj = 7;
93 vv = 4.3;
94 PetscCall(MatSetValues(A, 1, &ii, 1, &jj, &vv, INSERT_VALUES));
95 }
96 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
97 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
98
99 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
100
101 PetscCall(MatDestroy(&A));
102 PetscFunctionReturn(PETSC_SUCCESS);
103 }
104
ex2_square_bsvariable(void)105 PetscErrorCode ex2_square_bsvariable(void)
106 {
107 Mat A, preallocator;
108 PetscInt M, N, m, n, bs = 1;
109
110 PetscFunctionBegin;
111 PetscCall(PetscOptionsGetInt(NULL, NULL, "-block_size", &bs, NULL));
112
113 /*
114 Create the Jacobian matrix.
115 */
116 M = 10 * bs;
117 N = 10 * bs;
118 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
119 PetscCall(MatSetType(A, MATAIJ));
120 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N));
121 PetscCall(MatSetBlockSize(A, bs));
122 PetscCall(MatSetFromOptions(A));
123
124 /*
125 Get the sizes of the jacobian.
126 */
127 PetscCall(MatGetLocalSize(A, &m, &n));
128 PetscCall(MatGetBlockSize(A, &bs));
129
130 /*
131 Create a preallocator matrix with dimensions matching the jacobian A.
132 */
133 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &preallocator));
134 PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
135 PetscCall(MatSetSizes(preallocator, m, n, M, N));
136 PetscCall(MatSetBlockSize(preallocator, bs));
137 PetscCall(MatSetUp(preallocator));
138
139 /*
140 Insert non-zero pattern (e.g. perform a sweep over the grid).
141 You can use MatSetValues(), MatSetValuesBlocked() or MatSetValue().
142 */
143 {
144 PetscInt ii, jj;
145 PetscScalar *vv;
146
147 PetscCall(PetscCalloc1(bs * bs, &vv));
148
149 ii = 0;
150 jj = 9;
151 PetscCall(MatSetValue(preallocator, ii, jj, vv[0], INSERT_VALUES));
152
153 ii = 0;
154 jj = 0;
155 PetscCall(MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES));
156
157 ii = 2;
158 jj = 4;
159 PetscCall(MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES));
160
161 ii = 4;
162 jj = 2;
163 PetscCall(MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES));
164
165 ii = 4;
166 jj = 4;
167 PetscCall(MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES));
168
169 ii = 9;
170 jj = 9;
171 PetscCall(MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES));
172
173 PetscCall(PetscFree(vv));
174 }
175 PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
176 PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
177
178 /*
179 Push non-zero pattern defined from preallocator into A.
180 Internally this will call MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE).
181 The arg fill = PETSC_TRUE will insert zeros in the matrix A automatically.
182 */
183 PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_TRUE, A));
184
185 /*
186 We no longer require the preallocator object so destroy it.
187 */
188 PetscCall(MatDestroy(&preallocator));
189
190 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
191
192 {
193 PetscInt ii, jj;
194 PetscScalar *vv;
195
196 PetscCall(PetscCalloc1(bs * bs, &vv));
197 for (ii = 0; ii < bs * bs; ii++) vv[ii] = (PetscReal)(ii + 1);
198
199 ii = 0;
200 jj = 9;
201 PetscCall(MatSetValue(A, ii, jj, 33.3, INSERT_VALUES));
202
203 ii = 0;
204 jj = 0;
205 PetscCall(MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES));
206
207 ii = 2;
208 jj = 4;
209 PetscCall(MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES));
210
211 ii = 4;
212 jj = 2;
213 PetscCall(MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES));
214
215 ii = 4;
216 jj = 4;
217 PetscCall(MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES));
218
219 ii = 9;
220 jj = 9;
221 PetscCall(MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES));
222
223 PetscCall(PetscFree(vv));
224 }
225 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
226 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
227
228 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
229
230 PetscCall(MatDestroy(&A));
231 PetscFunctionReturn(PETSC_SUCCESS);
232 }
233
main(int argc,char ** args)234 int main(int argc, char **args)
235 {
236 PetscInt testid = 0;
237
238 PetscFunctionBeginUser;
239 PetscCall(PetscInitialize(&argc, &args, NULL, help));
240 PetscCall(PetscOptionsGetInt(NULL, NULL, "-test_id", &testid, NULL));
241 switch (testid) {
242 case 0:
243 PetscCall(ex1_nonsquare_bs1());
244 break;
245 case 1:
246 PetscCall(ex2_square_bsvariable());
247 break;
248 default:
249 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Invalid value for -test_id. Must be {0,1}");
250 }
251 PetscCall(PetscFinalize());
252 return 0;
253 }
254
255 /*TEST
256
257 test:
258 suffix: t0_a_aij
259 nsize: 1
260 args: -test_id 0 -mat_type aij
261
262 test:
263 suffix: t0_b_aij
264 nsize: 6
265 args: -test_id 0 -mat_type aij
266
267 test:
268 suffix: t1_a_aij
269 nsize: 1
270 args: -test_id 1 -mat_type aij
271
272 test:
273 suffix: t2_a_baij
274 nsize: 1
275 args: -test_id 1 -mat_type baij
276
277 test:
278 suffix: t3_a_sbaij
279 nsize: 1
280 args: -test_id 1 -mat_type sbaij
281
282 test:
283 suffix: t4_a_aij_bs3
284 nsize: 1
285 args: -test_id 1 -mat_type aij -block_size 3
286
287 test:
288 suffix: t5_a_baij_bs3
289 nsize: 1
290 args: -test_id 1 -mat_type baij -block_size 3
291
292 test:
293 suffix: t6_a_sbaij_bs3
294 nsize: 1
295 args: -test_id 1 -mat_type sbaij -block_size 3
296
297 test:
298 suffix: t4_b_aij_bs3
299 nsize: 6
300 args: -test_id 1 -mat_type aij -block_size 3
301
302 test:
303 suffix: t5_b_baij_bs3
304 nsize: 6
305 args: -test_id 1 -mat_type baij -block_size 3
306
307 test:
308 suffix: t6_b_sbaij_bs3
309 nsize: 6
310 args: -test_id 1 -mat_type sbaij -block_size 3
311
312 TEST*/
313