Actual source code: ex230.c

  1: static char help[] = "Example of using MatPreallocator\n\n";

  3: /*T
  4:    Concepts: Mat^matrix preallocation
  5:    Processors: n
  6: T*/

  8: #include <petscmat.h>

 10: PetscErrorCode ex1_nonsquare_bs1(void)
 11: {
 12:   Mat            A,preallocator;
 13:   PetscInt       M,N,m,n,bs;

 15:   /*
 16:      Create the Jacobian matrix
 17:   */
 18:   M = 10;
 19:   N = 8;
 20:   MatCreate(PETSC_COMM_WORLD,&A);
 21:   MatSetType(A,MATAIJ);
 22:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
 23:   MatSetBlockSize(A,1);
 24:   MatSetFromOptions(A);

 26:   /*
 27:      Get the sizes of the jacobian.
 28:   */
 29:   MatGetLocalSize(A,&m,&n);
 30:   MatGetBlockSize(A,&bs);

 32:   /*
 33:      Create a preallocator matrix with sizes (local and global) matching the jacobian A.
 34:   */
 35:   MatCreate(PetscObjectComm((PetscObject)A),&preallocator);
 36:   MatSetType(preallocator,MATPREALLOCATOR);
 37:   MatSetSizes(preallocator,m,n,M,N);
 38:   MatSetBlockSize(preallocator,bs);
 39:   MatSetUp(preallocator);

 41:   /*
 42:      Insert non-zero pattern (e.g. perform a sweep over the grid).
 43:      You can use MatSetValues(), MatSetValuesBlocked() or MatSetValue().
 44:   */
 45:   {
 46:     PetscInt    ii,jj;
 47:     PetscScalar vv = 0.0;

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

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

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

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

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

 73:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

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

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

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

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

 94:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

 96:   MatDestroy(&A);
 97:   return 0;
 98: }

100: PetscErrorCode ex2_square_bsvariable(void)
101: {
102:   Mat            A,preallocator;
103:   PetscInt       M,N,m,n,bs = 1;

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

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

118:   /*
119:      Get the sizes of the jacobian.
120:   */
121:   MatGetLocalSize(A,&m,&n);
122:   MatGetBlockSize(A,&bs);

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

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

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

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

146:     ii = 0; jj = 0;
147:     MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);

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

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

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

158:     ii = 9; jj = 9;
159:     MatSetValuesBlocked(preallocator,1,&ii,1,&jj,vv,INSERT_VALUES);

161:     PetscFree(vv);
162:   }
163:   MatAssemblyBegin(preallocator,MAT_FINAL_ASSEMBLY);
164:   MatAssemblyEnd(preallocator,MAT_FINAL_ASSEMBLY);

166:   /*
167:      Push non-zero pattern defined from preallocator into A.
168:      Internally this will call MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE).
169:      The arg fill = PETSC_TRUE will insert zeros in the matrix A automatically.
170:   */
171:   MatPreallocatorPreallocate(preallocator,PETSC_TRUE,A);

173:   /*
174:      We no longer require the preallocator object so destroy it.
175:   */
176:   MatDestroy(&preallocator);

178:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

180:   {
181:     PetscInt    ii,jj;
182:     PetscScalar *vv;

184:     PetscCalloc1(bs*bs,&vv);
185:     for (ii=0; ii<bs*bs; ii++) {
186:       vv[ii] = (PetscReal)(ii+1);
187:     }

189:     ii = 0; jj = 9;
190:     MatSetValue(A,ii,jj,33.3,INSERT_VALUES);

192:     ii = 0; jj = 0;
193:     MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);

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

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

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

204:     ii = 9; jj = 9;
205:     MatSetValuesBlocked(A,1,&ii,1,&jj,vv,INSERT_VALUES);

207:     PetscFree(vv);
208:   }
209:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
210:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

212:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

214:   MatDestroy(&A);
215:   return 0;
216: }

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

237: /*TEST

239:    test:
240:      suffix: t0_a_aij
241:      nsize: 1
242:      args: -test_id 0 -mat_type aij

244:    test:
245:      suffix: t0_b_aij
246:      nsize: 6
247:      args: -test_id 0 -mat_type aij

249:    test:
250:      suffix: t1_a_aij
251:      nsize: 1
252:      args: -test_id 1 -mat_type aij

254:    test:
255:      suffix: t2_a_baij
256:      nsize: 1
257:      args: -test_id 1 -mat_type baij

259:    test:
260:      suffix: t3_a_sbaij
261:      nsize: 1
262:      args: -test_id 1 -mat_type sbaij

264:    test:
265:      suffix: t4_a_aij_bs3
266:      nsize: 1
267:      args: -test_id 1 -mat_type aij -block_size 3

269:    test:
270:      suffix: t5_a_baij_bs3
271:      nsize: 1
272:      args: -test_id 1 -mat_type baij -block_size 3

274:    test:
275:      suffix: t6_a_sbaij_bs3
276:      nsize: 1
277:      args: -test_id 1 -mat_type sbaij -block_size 3

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

284:    test:
285:      suffix: t5_b_baij_bs3
286:      nsize: 6
287:      args: -test_id 1 -mat_type baij -block_size 3
288:      filter: grep -v Mat_

290:    test:
291:      suffix: t6_b_sbaij_bs3
292:      nsize: 6
293:      args: -test_id 1 -mat_type sbaij -block_size 3
294:      filter: grep -v Mat_

296: TEST*/