Actual source code: ex230.c
1: static char help[] = "Example of using MatPreallocator\n\n";
3: #include <petscmat.h>
5: PetscErrorCode ex1_nonsquare_bs1(void)
6: {
7: Mat A, preallocator;
8: PetscInt M, N, m, n, bs;
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));
22: /*
23: Get the sizes of the jacobian.
24: */
25: PetscCall(MatGetLocalSize(A, &m, &n));
26: PetscCall(MatGetBlockSize(A, &bs));
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));
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;
45: ii = 3;
46: jj = 3;
47: PetscCall(MatSetValues(preallocator, 1, &ii, 1, &jj, &vv, INSERT_VALUES));
49: ii = 7;
50: jj = 4;
51: PetscCall(MatSetValues(preallocator, 1, &ii, 1, &jj, &vv, INSERT_VALUES));
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));
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));
67: /*
68: We no longer require the preallocator object so destroy it.
69: */
70: PetscCall(MatDestroy(&preallocator));
72: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
74: /*
75: Insert non-zero values into A.
76: */
77: {
78: PetscInt ii, jj;
79: PetscScalar vv;
81: ii = 3;
82: jj = 3;
83: vv = 0.3;
84: PetscCall(MatSetValue(A, ii, jj, vv, INSERT_VALUES));
86: ii = 7;
87: jj = 4;
88: vv = 3.3;
89: PetscCall(MatSetValues(A, 1, &ii, 1, &jj, &vv, INSERT_VALUES));
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));
99: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
101: PetscCall(MatDestroy(&A));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: PetscErrorCode ex2_square_bsvariable(void)
106: {
107: Mat A, preallocator;
108: PetscInt M, N, m, n, bs = 1;
110: PetscFunctionBegin;
111: PetscCall(PetscOptionsGetInt(NULL, NULL, "-block_size", &bs, NULL));
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));
124: /*
125: Get the sizes of the jacobian.
126: */
127: PetscCall(MatGetLocalSize(A, &m, &n));
128: PetscCall(MatGetBlockSize(A, &bs));
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));
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;
147: PetscCall(PetscCalloc1(bs * bs, &vv));
149: ii = 0;
150: jj = 9;
151: PetscCall(MatSetValue(preallocator, ii, jj, vv[0], INSERT_VALUES));
153: ii = 0;
154: jj = 0;
155: PetscCall(MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES));
157: ii = 2;
158: jj = 4;
159: PetscCall(MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES));
161: ii = 4;
162: jj = 2;
163: PetscCall(MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES));
165: ii = 4;
166: jj = 4;
167: PetscCall(MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES));
169: ii = 9;
170: jj = 9;
171: PetscCall(MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES));
173: PetscCall(PetscFree(vv));
174: }
175: PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
176: PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
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));
185: /*
186: We no longer require the preallocator object so destroy it.
187: */
188: PetscCall(MatDestroy(&preallocator));
190: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
192: {
193: PetscInt ii, jj;
194: PetscScalar *vv;
196: PetscCall(PetscCalloc1(bs * bs, &vv));
197: for (ii = 0; ii < bs * bs; ii++) vv[ii] = (PetscReal)(ii + 1);
199: ii = 0;
200: jj = 9;
201: PetscCall(MatSetValue(A, ii, jj, 33.3, INSERT_VALUES));
203: ii = 0;
204: jj = 0;
205: PetscCall(MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES));
207: ii = 2;
208: jj = 4;
209: PetscCall(MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES));
211: ii = 4;
212: jj = 2;
213: PetscCall(MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES));
215: ii = 4;
216: jj = 4;
217: PetscCall(MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES));
219: ii = 9;
220: jj = 9;
221: PetscCall(MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES));
223: PetscCall(PetscFree(vv));
224: }
225: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
226: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
228: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
230: PetscCall(MatDestroy(&A));
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: int main(int argc, char **args)
235: {
236: PetscInt testid = 0;
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: }
255: /*TEST
257: test:
258: suffix: t0_a_aij
259: nsize: 1
260: args: -test_id 0 -mat_type aij
262: test:
263: suffix: t0_b_aij
264: nsize: 6
265: args: -test_id 0 -mat_type aij
267: test:
268: suffix: t1_a_aij
269: nsize: 1
270: args: -test_id 1 -mat_type aij
272: test:
273: suffix: t2_a_baij
274: nsize: 1
275: args: -test_id 1 -mat_type baij
277: test:
278: suffix: t3_a_sbaij
279: nsize: 1
280: args: -test_id 1 -mat_type sbaij
282: test:
283: suffix: t4_a_aij_bs3
284: nsize: 1
285: args: -test_id 1 -mat_type aij -block_size 3
287: test:
288: suffix: t5_a_baij_bs3
289: nsize: 1
290: args: -test_id 1 -mat_type baij -block_size 3
292: test:
293: suffix: t6_a_sbaij_bs3
294: nsize: 1
295: args: -test_id 1 -mat_type sbaij -block_size 3
297: test:
298: suffix: t4_b_aij_bs3
299: nsize: 6
300: args: -test_id 1 -mat_type aij -block_size 3
302: test:
303: suffix: t5_b_baij_bs3
304: nsize: 6
305: args: -test_id 1 -mat_type baij -block_size 3
307: test:
308: suffix: t6_b_sbaij_bs3
309: nsize: 6
310: args: -test_id 1 -mat_type sbaij -block_size 3
312: TEST*/