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;
237: PetscFunctionBeginUser;
238: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
239: PetscCall(PetscOptionsGetInt(NULL, NULL, "-test_id", &testid, NULL));
240: switch (testid) {
241: case 0:
242: PetscCall(ex1_nonsquare_bs1());
243: break;
244: case 1:
245: PetscCall(ex2_square_bsvariable());
246: break;
247: default:
248: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Invalid value for -test_id. Must be {0,1}");
249: }
250: PetscCall(PetscFinalize());
251: return 0;
252: }
254: /*TEST
256: test:
257: suffix: t0_a_aij
258: nsize: 1
259: args: -test_id 0 -mat_type aij
261: test:
262: suffix: t0_b_aij
263: nsize: 6
264: args: -test_id 0 -mat_type aij
266: test:
267: suffix: t1_a_aij
268: nsize: 1
269: args: -test_id 1 -mat_type aij
271: test:
272: suffix: t2_a_baij
273: nsize: 1
274: args: -test_id 1 -mat_type baij
276: test:
277: suffix: t3_a_sbaij
278: nsize: 1
279: args: -test_id 1 -mat_type sbaij
281: test:
282: suffix: t4_a_aij_bs3
283: nsize: 1
284: args: -test_id 1 -mat_type aij -block_size 3
286: test:
287: suffix: t5_a_baij_bs3
288: nsize: 1
289: args: -test_id 1 -mat_type baij -block_size 3
291: test:
292: suffix: t6_a_sbaij_bs3
293: nsize: 1
294: args: -test_id 1 -mat_type sbaij -block_size 3
296: test:
297: suffix: t4_b_aij_bs3
298: nsize: 6
299: args: -test_id 1 -mat_type aij -block_size 3
301: test:
302: suffix: t5_b_baij_bs3
303: nsize: 6
304: args: -test_id 1 -mat_type baij -block_size 3
306: test:
307: suffix: t6_b_sbaij_bs3
308: nsize: 6
309: args: -test_id 1 -mat_type sbaij -block_size 3
311: TEST*/