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; jj = 3;
 45:     MatSetValues(preallocator,1,&ii,1,&jj,&vv,INSERT_VALUES);

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

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

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

 63:   /*
 64:      We no longer require the preallocator object so destroy it.
 65:   */
 66:   MatDestroy(&preallocator);

 68:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

 70:   /*
 71:      Insert non-zero values into A.
 72:   */
 73:   {
 74:     PetscInt    ii,jj;
 75:     PetscScalar vv;

 77:     ii = 3; jj = 3; vv = 0.3;
 78:     MatSetValue(A,ii,jj,vv,INSERT_VALUES);

 80:     ii = 7; jj = 4; vv = 3.3;
 81:     MatSetValues(A,1,&ii,1,&jj,&vv,INSERT_VALUES);

 83:     ii = 9; jj = 7; vv = 4.3;
 84:     MatSetValues(A,1,&ii,1,&jj,&vv,INSERT_VALUES);
 85:   }
 86:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 87:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 89:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

 91:   MatDestroy(&A);
 92:   return 0;
 93: }

 95: PetscErrorCode ex2_square_bsvariable(void)
 96: {
 97:   Mat            A,preallocator;
 98:   PetscInt       M,N,m,n,bs = 1;

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

102:   /*
103:      Create the Jacobian matrix.
104:   */
105:   M = 10 * bs;
106:   N = 10 * bs;
107:   MatCreate(PETSC_COMM_WORLD,&A);
108:   MatSetType(A,MATAIJ);
109:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
110:   MatSetBlockSize(A,bs);
111:   MatSetFromOptions(A);

113:   /*
114:      Get the sizes of the jacobian.
115:   */
116:   MatGetLocalSize(A,&m,&n);
117:   MatGetBlockSize(A,&bs);

119:   /*
120:      Create a preallocator matrix with dimensions matching the jacobian A.
121:   */
122:   MatCreate(PetscObjectComm((PetscObject)A),&preallocator);
123:   MatSetType(preallocator,MATPREALLOCATOR);
124:   MatSetSizes(preallocator,m,n,M,N);
125:   MatSetBlockSize(preallocator,bs);
126:   MatSetUp(preallocator);

128:   /*
129:      Insert non-zero pattern (e.g. perform a sweep over the grid).
130:      You can use MatSetValues(), MatSetValuesBlocked() or MatSetValue().
131:   */
132:   {
133:     PetscInt  ii,jj;
134:     PetscScalar *vv;

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

138:     ii = 0; jj = 9;
139:     MatSetValue(preallocator,ii,jj,vv[0],INSERT_VALUES);

141:     ii = 0; jj = 0;
142:     MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);

144:     ii = 2; jj = 4;
145:     MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);

147:     ii = 4; jj = 2;
148:     MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);

150:     ii = 4; jj = 4;
151:     MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);

153:     ii = 9; jj = 9;
154:     MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);

156:     PetscFree(vv);
157:   }
158:   MatAssemblyBegin(preallocator,MAT_FINAL_ASSEMBLY);
159:   MatAssemblyEnd(preallocator,MAT_FINAL_ASSEMBLY);

161:   /*
162:      Push non-zero pattern defined from preallocator into A.
163:      Internally this will call MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE).
164:      The arg fill = PETSC_TRUE will insert zeros in the matrix A automatically.
165:   */
166:   MatPreallocatorPreallocate(preallocator,PETSC_TRUE,A);

168:   /*
169:      We no longer require the preallocator object so destroy it.
170:   */
171:   MatDestroy(&preallocator);

173:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

175:   {
176:     PetscInt    ii,jj;
177:     PetscScalar *vv;

179:     PetscCalloc1(bs*bs,&vv);
180:     for (ii=0; ii<bs*bs; ii++) {
181:       vv[ii] = (PetscReal)(ii+1);
182:     }

184:     ii = 0; jj = 9;
185:     MatSetValue(A,ii,jj,33.3,INSERT_VALUES);

187:     ii = 0; jj = 0;
188:     MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);

190:     ii = 2; jj = 4;
191:     MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);

193:     ii = 4; jj = 2;
194:     MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);

196:     ii = 4; jj = 4;
197:     MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);

199:     ii = 9; jj = 9;
200:     MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);

202:     PetscFree(vv);
203:   }
204:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
205:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

207:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

209:   MatDestroy(&A);
210:   return 0;
211: }

213: int main(int argc, char **args)
214: {
215:   PetscInt testid = 0;
216:   PetscInitialize(&argc,&args,(char*)0,help);
217:   PetscOptionsGetInt(NULL,NULL,"-test_id",&testid,NULL);
218:   switch (testid) {
219:     case 0:
220:       ex1_nonsquare_bs1();
221:       break;
222:     case 1:
223:       ex2_square_bsvariable();
224:       break;
225:     default:
226:       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Invalid value for -test_id. Must be {0,1}");
227:   }
228:   PetscFinalize();
229:   return 0;
230: }

232: /*TEST

234:    test:
235:      suffix: t0_a_aij
236:      nsize: 1
237:      args: -test_id 0 -mat_type aij

239:    test:
240:      suffix: t0_b_aij
241:      nsize: 6
242:      args: -test_id 0 -mat_type aij

244:    test:
245:      suffix: t1_a_aij
246:      nsize: 1
247:      args: -test_id 1 -mat_type aij

249:    test:
250:      suffix: t2_a_baij
251:      nsize: 1
252:      args: -test_id 1 -mat_type baij

254:    test:
255:      suffix: t3_a_sbaij
256:      nsize: 1
257:      args: -test_id 1 -mat_type sbaij

259:    test:
260:      suffix: t4_a_aij_bs3
261:      nsize: 1
262:      args: -test_id 1 -mat_type aij -block_size 3

264:    test:
265:      suffix: t5_a_baij_bs3
266:      nsize: 1
267:      args: -test_id 1 -mat_type baij -block_size 3

269:    test:
270:      suffix: t6_a_sbaij_bs3
271:      nsize: 1
272:      args: -test_id 1 -mat_type sbaij -block_size 3

274:    test:
275:      suffix: t4_b_aij_bs3
276:      nsize: 6
277:      args: -test_id 1 -mat_type aij -block_size 3

279:    test:
280:      suffix: t5_b_baij_bs3
281:      nsize: 6
282:      args: -test_id 1 -mat_type baij -block_size 3
283:      filter: grep -v Mat_

285:    test:
286:      suffix: t6_b_sbaij_bs3
287:      nsize: 6
288:      args: -test_id 1 -mat_type sbaij -block_size 3
289:      filter: grep -v Mat_

291: TEST*/