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*/