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:   M = 10;
 14:   N = 8;
 15:   MatCreate(PETSC_COMM_WORLD, &A);
 16:   MatSetType(A, MATAIJ);
 17:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N);
 18:   MatSetBlockSize(A, 1);
 19:   MatSetFromOptions(A);

 21:   /*
 22:      Get the sizes of the jacobian.
 23:   */
 24:   MatGetLocalSize(A, &m, &n);
 25:   MatGetBlockSize(A, &bs);

 27:   /*
 28:      Create a preallocator matrix with sizes (local and global) matching the jacobian A.
 29:   */
 30:   MatCreate(PetscObjectComm((PetscObject)A), &preallocator);
 31:   MatSetType(preallocator, MATPREALLOCATOR);
 32:   MatSetSizes(preallocator, m, n, M, N);
 33:   MatSetBlockSize(preallocator, bs);
 34:   MatSetUp(preallocator);

 36:   /*
 37:      Insert non-zero pattern (e.g. perform a sweep over the grid).
 38:      You can use MatSetValues(), MatSetValuesBlocked() or MatSetValue().
 39:   */
 40:   {
 41:     PetscInt    ii, jj;
 42:     PetscScalar vv = 0.0;

 44:     ii = 3;
 45:     jj = 3;
 46:     MatSetValues(preallocator, 1, &ii, 1, &jj, &vv, INSERT_VALUES);

 48:     ii = 7;
 49:     jj = 4;
 50:     MatSetValues(preallocator, 1, &ii, 1, &jj, &vv, INSERT_VALUES);

 52:     ii = 9;
 53:     jj = 7;
 54:     MatSetValue(preallocator, ii, jj, vv, INSERT_VALUES);
 55:   }
 56:   MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
 57:   MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);

 59:   /*
 60:      Push the non-zero pattern defined within preallocator into A.
 61:      Internally this will call MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE).
 62:      The arg fill = PETSC_TRUE will insert zeros in the matrix A automatically.
 63:   */
 64:   MatPreallocatorPreallocate(preallocator, PETSC_TRUE, A);

 66:   /*
 67:      We no longer require the preallocator object so destroy it.
 68:   */
 69:   MatDestroy(&preallocator);

 71:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);

 73:   /*
 74:      Insert non-zero values into A.
 75:   */
 76:   {
 77:     PetscInt    ii, jj;
 78:     PetscScalar vv;

 80:     ii = 3;
 81:     jj = 3;
 82:     vv = 0.3;
 83:     MatSetValue(A, ii, jj, vv, INSERT_VALUES);

 85:     ii = 7;
 86:     jj = 4;
 87:     vv = 3.3;
 88:     MatSetValues(A, 1, &ii, 1, &jj, &vv, INSERT_VALUES);

 90:     ii = 9;
 91:     jj = 7;
 92:     vv = 4.3;
 93:     MatSetValues(A, 1, &ii, 1, &jj, &vv, INSERT_VALUES);
 94:   }
 95:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 96:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

 98:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);

100:   MatDestroy(&A);
101:   return 0;
102: }

104: PetscErrorCode ex2_square_bsvariable(void)
105: {
106:   Mat      A, preallocator;
107:   PetscInt M, N, m, n, bs = 1;

109:   PetscOptionsGetInt(NULL, NULL, "-block_size", &bs, NULL);

111:   /*
112:      Create the Jacobian matrix.
113:   */
114:   M = 10 * bs;
115:   N = 10 * bs;
116:   MatCreate(PETSC_COMM_WORLD, &A);
117:   MatSetType(A, MATAIJ);
118:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N);
119:   MatSetBlockSize(A, bs);
120:   MatSetFromOptions(A);

122:   /*
123:      Get the sizes of the jacobian.
124:   */
125:   MatGetLocalSize(A, &m, &n);
126:   MatGetBlockSize(A, &bs);

128:   /*
129:      Create a preallocator matrix with dimensions matching the jacobian A.
130:   */
131:   MatCreate(PetscObjectComm((PetscObject)A), &preallocator);
132:   MatSetType(preallocator, MATPREALLOCATOR);
133:   MatSetSizes(preallocator, m, n, M, N);
134:   MatSetBlockSize(preallocator, bs);
135:   MatSetUp(preallocator);

137:   /*
138:      Insert non-zero pattern (e.g. perform a sweep over the grid).
139:      You can use MatSetValues(), MatSetValuesBlocked() or MatSetValue().
140:   */
141:   {
142:     PetscInt     ii, jj;
143:     PetscScalar *vv;

145:     PetscCalloc1(bs * bs, &vv);

147:     ii = 0;
148:     jj = 9;
149:     MatSetValue(preallocator, ii, jj, vv[0], INSERT_VALUES);

151:     ii = 0;
152:     jj = 0;
153:     MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES);

155:     ii = 2;
156:     jj = 4;
157:     MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES);

159:     ii = 4;
160:     jj = 2;
161:     MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES);

163:     ii = 4;
164:     jj = 4;
165:     MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES);

167:     ii = 9;
168:     jj = 9;
169:     MatSetValuesBlocked(preallocator, 1, &ii, 1, &jj, vv, INSERT_VALUES);

171:     PetscFree(vv);
172:   }
173:   MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
174:   MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);

176:   /*
177:      Push non-zero pattern defined from preallocator into A.
178:      Internally this will call MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE).
179:      The arg fill = PETSC_TRUE will insert zeros in the matrix A automatically.
180:   */
181:   MatPreallocatorPreallocate(preallocator, PETSC_TRUE, A);

183:   /*
184:      We no longer require the preallocator object so destroy it.
185:   */
186:   MatDestroy(&preallocator);

188:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);

190:   {
191:     PetscInt     ii, jj;
192:     PetscScalar *vv;

194:     PetscCalloc1(bs * bs, &vv);
195:     for (ii = 0; ii < bs * bs; ii++) vv[ii] = (PetscReal)(ii + 1);

197:     ii = 0;
198:     jj = 9;
199:     MatSetValue(A, ii, jj, 33.3, INSERT_VALUES);

201:     ii = 0;
202:     jj = 0;
203:     MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES);

205:     ii = 2;
206:     jj = 4;
207:     MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES);

209:     ii = 4;
210:     jj = 2;
211:     MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES);

213:     ii = 4;
214:     jj = 4;
215:     MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES);

217:     ii = 9;
218:     jj = 9;
219:     MatSetValuesBlocked(A, 1, &ii, 1, &jj, vv, INSERT_VALUES);

221:     PetscFree(vv);
222:   }
223:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
224:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

226:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);

228:   MatDestroy(&A);
229:   return 0;
230: }

232: int main(int argc, char **args)
233: {
234:   PetscInt testid = 0;
236:   PetscInitialize(&argc, &args, (char *)0, help);
237:   PetscOptionsGetInt(NULL, NULL, "-test_id", &testid, NULL);
238:   switch (testid) {
239:   case 0:
240:     ex1_nonsquare_bs1();
241:     break;
242:   case 1:
243:     ex2_square_bsvariable();
244:     break;
245:   default:
246:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Invalid value for -test_id. Must be {0,1}");
247:   }
248:   PetscFinalize();
249:   return 0;
250: }

252: /*TEST

254:    test:
255:      suffix: t0_a_aij
256:      nsize: 1
257:      args: -test_id 0 -mat_type aij

259:    test:
260:      suffix: t0_b_aij
261:      nsize: 6
262:      args: -test_id 0 -mat_type aij

264:    test:
265:      suffix: t1_a_aij
266:      nsize: 1
267:      args: -test_id 1 -mat_type aij

269:    test:
270:      suffix: t2_a_baij
271:      nsize: 1
272:      args: -test_id 1 -mat_type baij

274:    test:
275:      suffix: t3_a_sbaij
276:      nsize: 1
277:      args: -test_id 1 -mat_type sbaij

279:    test:
280:      suffix: t4_a_aij_bs3
281:      nsize: 1
282:      args: -test_id 1 -mat_type aij -block_size 3

284:    test:
285:      suffix: t5_a_baij_bs3
286:      nsize: 1
287:      args: -test_id 1 -mat_type baij -block_size 3

289:    test:
290:      suffix: t6_a_sbaij_bs3
291:      nsize: 1
292:      args: -test_id 1 -mat_type sbaij -block_size 3

294:    test:
295:      suffix: t4_b_aij_bs3
296:      nsize: 6
297:      args: -test_id 1 -mat_type aij -block_size 3

299:    test:
300:      suffix: t5_b_baij_bs3
301:      nsize: 6
302:      args: -test_id 1 -mat_type baij -block_size 3
303:      filter: grep -v Mat_

305:    test:
306:      suffix: t6_b_sbaij_bs3
307:      nsize: 6
308:      args: -test_id 1 -mat_type sbaij -block_size 3
309:      filter: grep -v Mat_

311: TEST*/