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