Actual source code: ex123.c

  1: static char help[] = "Test MatSetPreallocationCOO and MatSetValuesCOO\n\n";

  3: #include <petscmat.h>
  4: #define MyMatView(a,b) (PetscPrintf(PetscObjectComm((PetscObject)(a)),"LINE %d\n",__LINE__),MatView(a,b));
  5: #define MyVecView(a,b) (PetscPrintf(PetscObjectComm((PetscObject)(a)),"LINE %d\n",__LINE__),VecView(a,b));
  6: int main(int argc,char **args)
  7: {
  8:   Mat                    A,At,AAt;
  9:   Vec                    x,y,z;
 10:   ISLocalToGlobalMapping rl2g,cl2g;
 11:   IS                     is;
 12:   PetscLayout            rmap,cmap;
 13:   PetscInt               *it,*jt;
 14:   PetscInt               n1 = 11, n2 = 9;
 15:   PetscInt               i1[] = {   7,  6,  2,  0,  4,  1,  1,  0,  2,  2,  1 , -1, -1};
 16:   PetscInt               j1[] = {   1,  4,  3,  5,  3,  3,  4,  5,  0,  3,  1 , -1, -1};
 17:   PetscInt               i2[] = {   7,  6,  2,  0,  4,  1,  1,  2,  1, -1, -1};
 18:   PetscInt               j2[] = {   1,  4,  3,  5,  3,  3,  4,  0,  1, -1, -1};
 19:   PetscScalar            v1[] = { -1., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., PETSC_MAX_REAL, PETSC_MAX_REAL};
 20:   PetscScalar            v2[] = {  1.,-1.,-2.,-3.,-4.,-5.,-6.,-7.,-8.,-9.,-10., PETSC_MAX_REAL, PETSC_MAX_REAL};
 21:   PetscInt               N = 6, m = 8, M, rstart, cstart, i;
 22:   PetscMPIInt            size;
 23:   PetscBool              loc = PETSC_FALSE;
 24:   PetscBool              locdiag = PETSC_TRUE;
 25:   PetscBool              localapi = PETSC_FALSE;
 26:   PetscBool              neg = PETSC_FALSE;
 27:   PetscBool              ismatis, ismpiaij;

 29:   PetscInitialize(&argc,&args,(char*)0,help);
 30:   PetscOptionsGetBool(NULL,NULL,"-neg",&neg,NULL);
 31:   PetscOptionsGetBool(NULL,NULL,"-loc",&loc,NULL);
 32:   PetscOptionsGetBool(NULL,NULL,"-locdiag",&locdiag,NULL);
 33:   PetscOptionsGetBool(NULL,NULL,"-localapi",&localapi,NULL);
 34:   MatCreate(PETSC_COMM_WORLD,&A);
 35:   if (loc) {
 36:     if (locdiag) {
 37:       MatSetSizes(A,m,N,PETSC_DECIDE,PETSC_DECIDE);
 38:     } else {
 39:       MatSetSizes(A,m,m+N,PETSC_DECIDE,PETSC_DECIDE);
 40:     }
 41:   } else {
 42:     MatSetSizes(A,m,PETSC_DECIDE,PETSC_DECIDE,N);
 43:   }
 44:   MatSetFromOptions(A);
 45:   MatGetLayouts(A,&rmap,&cmap);
 46:   PetscLayoutSetUp(rmap);
 47:   PetscLayoutSetUp(cmap);
 48:   PetscLayoutGetRange(rmap,&rstart,NULL);
 49:   PetscLayoutGetRange(cmap,&cstart,NULL);
 50:   PetscLayoutGetSize(rmap,&M);
 51:   PetscLayoutGetSize(cmap,&N);

 53:   PetscObjectTypeCompare((PetscObject)A,MATIS,&ismatis);

 55:   /* create fake l2g maps to test the local API */
 56:   ISCreateStride(PETSC_COMM_WORLD,M-rstart,rstart,1,&is);
 57:   ISLocalToGlobalMappingCreateIS(is,&rl2g);
 58:   ISDestroy(&is);
 59:   ISCreateStride(PETSC_COMM_WORLD,N,0,1,&is);
 60:   ISLocalToGlobalMappingCreateIS(is,&cl2g);
 61:   ISDestroy(&is);
 62:   MatSetLocalToGlobalMapping(A,rl2g,cl2g);
 63:   ISLocalToGlobalMappingDestroy(&rl2g);
 64:   ISLocalToGlobalMappingDestroy(&cl2g);

 66:   MatCreateVecs(A,&x,&y);
 67:   MatCreateVecs(A,NULL,&z);
 68:   VecSet(x,1.);
 69:   VecSet(z,2.);
 70:   if (!localapi) for (i = 0; i < n1; i++) i1[i] += rstart;
 71:   if (!localapi) for (i = 0; i < n2; i++) i2[i] += rstart;
 72:   if (loc) {
 73:     if (locdiag) {
 74:       for (i = 0; i < n1; i++) j1[i] += cstart;
 75:       for (i = 0; i < n2; i++) j2[i] += cstart;
 76:     } else {
 77:       for (i = 0; i < n1; i++) j1[i] += cstart + m;
 78:       for (i = 0; i < n2; i++) j2[i] += cstart + m;
 79:     }
 80:   }
 81:   if (neg) { n1 += 2; n2 += 2; }
 82:   /* MatSetPreallocationCOOLocal maps the indices! */
 83:   PetscMalloc2(PetscMax(n1,n2),&it,PetscMax(n1,n2),&jt);
 84:   /* test with repeated entries */
 85:   if (!localapi) {
 86:     MatSetPreallocationCOO(A,n1,i1,j1);
 87:   } else {
 88:     PetscArraycpy(it,i1,n1);
 89:     PetscArraycpy(jt,j1,n1);
 90:     MatSetPreallocationCOOLocal(A,n1,it,jt);
 91:   }
 92:   MatSetValuesCOO(A,v1,ADD_VALUES);
 93:   MatMult(A,x,y);
 94:   MyMatView(A,NULL);
 95:   MyVecView(y,NULL);
 96:   MatSetValuesCOO(A,v2,ADD_VALUES);
 97:   MatMultAdd(A,x,y,y);
 98:   MyMatView(A,NULL);
 99:   MyVecView(y,NULL);
100:   MatTranspose(A,MAT_INITIAL_MATRIX,&At);
101:   if (!ismatis) {
102:     MatMatMult(A,At,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);
103:     MyMatView(AAt,NULL);
104:     MatDestroy(&AAt);
105:     MatMatMult(At,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);
106:     MyMatView(AAt,NULL);
107:     MatDestroy(&AAt);
108:   }
109:   MatDestroy(&At);

111:   /* INSERT_VALUES will overwrite matrix entries but
112:      still perform the sum of the repeated entries */
113:   MatSetValuesCOO(A,v2,INSERT_VALUES);
114:   MyMatView(A,NULL);

116:   /* test with unique entries */
117:   if (!localapi) {
118:     MatSetPreallocationCOO(A,n2,i2,j2);
119:   } else {
120:     PetscArraycpy(it,i2,n2);
121:     PetscArraycpy(jt,j2,n2);
122:     MatSetPreallocationCOOLocal(A,n2,it,jt);
123:   }
124:   MatSetValuesCOO(A,v1,ADD_VALUES);
125:   MatMult(A,x,y);
126:   MyMatView(A,NULL);
127:   MyVecView(y,NULL);
128:   MatSetValuesCOO(A,v2,ADD_VALUES);
129:   MatMultAdd(A,x,y,z);
130:   MyMatView(A,NULL);
131:   MyVecView(z,NULL);
132:   if (!localapi) {
133:     MatSetPreallocationCOO(A,n2,i2,j2);
134:   } else {
135:     PetscArraycpy(it,i2,n2);
136:     PetscArraycpy(jt,j2,n2);
137:     MatSetPreallocationCOOLocal(A,n2,it,jt);
138:   }
139:   MatSetValuesCOO(A,v1,INSERT_VALUES);
140:   MatMult(A,x,y);
141:   MyMatView(A,NULL);
142:   MyVecView(y,NULL);
143:   MatSetValuesCOO(A,v2,INSERT_VALUES);
144:   MatMultAdd(A,x,y,z);
145:   MyMatView(A,NULL);
146:   MyVecView(z,NULL);
147:   MatTranspose(A,MAT_INITIAL_MATRIX,&At);
148:   if (!ismatis) {
149:     MatMatMult(A,At,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);
150:     MyMatView(AAt,NULL);
151:     MatDestroy(&AAt);
152:     MatMatMult(At,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);
153:     MyMatView(AAt,NULL);
154:     MatDestroy(&AAt);
155:   }
156:   MatDestroy(&At);

158:   /* test providing diagonal first, then offdiagonal */
159:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
160:   PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);
161:   if (ismpiaij && size > 1) {
162:     Mat               lA,lB;
163:     const PetscInt    *garray,*iA,*jA,*iB,*jB;
164:     const PetscScalar *vA,*vB;
165:     PetscScalar       *coo_v;
166:     PetscInt          *coo_i,*coo_j;
167:     PetscInt          i,j,nA,nB,nnz;
168:     PetscBool         flg;

170:     MatMPIAIJGetSeqAIJ(A,&lA,&lB,&garray);
171:     MatSeqAIJGetArrayRead(lA,&vA);
172:     MatSeqAIJGetArrayRead(lB,&vB);
173:     MatGetRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nA,&iA,&jA,&flg);
174:     MatGetRowIJ(lB,0,PETSC_FALSE,PETSC_FALSE,&nB,&iB,&jB,&flg);
175:     nnz  = iA[nA] + iB[nB];
176:     PetscMalloc3(nnz,&coo_i,nnz,&coo_j,nnz,&coo_v);
177:     nnz  = 0;
178:     for (i=0;i<nA;i++) {
179:       for (j=iA[i];j<iA[i+1];j++,nnz++) {
180:         coo_i[nnz] = i+rstart;
181:         coo_j[nnz] = jA[j]+cstart;
182:         coo_v[nnz] = vA[j];
183:       }
184:     }
185:     for (i=0;i<nB;i++) {
186:       for (j=iB[i];j<iB[i+1];j++,nnz++) {
187:         coo_i[nnz] = i+rstart;
188:         coo_j[nnz] = garray[jB[j]];
189:         coo_v[nnz] = vB[j];
190:       }
191:     }
192:     MatRestoreRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nA,&iA,&jA,&flg);
193:     MatRestoreRowIJ(lB,0,PETSC_FALSE,PETSC_FALSE,&nB,&iB,&jB,&flg);
194:     MatSeqAIJRestoreArrayRead(lA,&vA);
195:     MatSeqAIJRestoreArrayRead(lB,&vB);

197:     MatSetPreallocationCOO(A,nnz,coo_i,coo_j);
198:     MatSetValuesCOO(A,coo_v,ADD_VALUES);
199:     MatMult(A,x,y);
200:     MyMatView(A,NULL);
201:     MyVecView(y,NULL);
202:     MatSetValuesCOO(A,coo_v,INSERT_VALUES);
203:     MatMult(A,x,y);
204:     MyMatView(A,NULL);
205:     MyVecView(y,NULL);
206:     MatTranspose(A,MAT_INITIAL_MATRIX,&At);
207:     MatMatMult(A,At,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);
208:     MyMatView(AAt,NULL);
209:     MatDestroy(&AAt);
210:     MatMatMult(At,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AAt);
211:     MyMatView(AAt,NULL);
212:     MatDestroy(&AAt);
213:     MatDestroy(&At);

215:     PetscFree3(coo_i,coo_j,coo_v);
216:   }
217:   PetscFree2(it,jt);
218:   VecDestroy(&z);
219:   VecDestroy(&x);
220:   VecDestroy(&y);
221:   MatDestroy(&A);
222:   PetscFinalize();
223:   return 0;
224: }

226: /*TEST

228:    test:
229:      suffix: 1
230:      filter: grep -v type
231:      diff_args: -j
232:      args: -mat_type {{seqaij mpiaij}} -localapi {{0 1}} -neg {{0 1}}

234:    test:
235:      requires: cuda
236:      suffix: 1_cuda
237:      filter: grep -v type
238:      diff_args: -j
239:      args: -mat_type {{seqaijcusparse mpiaijcusparse}} -localapi {{0 1}} -neg {{0 1}}
240:      output_file: output/ex123_1.out

242:    test:
243:      requires: kokkos_kernels !sycl
244:      suffix: 1_kokkos
245:      filter: grep -v type
246:      diff_args: -j
247:      args: -mat_type {{seqaijkokkos mpiaijkokkos}} -localapi {{0 1}} -neg {{0 1}}
248:      output_file: output/ex123_1.out

250:    test:
251:      suffix: 2
252:      nsize: 7
253:      filter: grep -v type
254:      diff_args: -j
255:      args: -mat_type mpiaij -localapi {{0 1}} -neg {{0 1}}

257:    test:
258:      requires: cuda
259:      suffix: 2_cuda
260:      nsize: 7
261:      filter: grep -v type
262:      diff_args: -j
263:      args: -mat_type mpiaijcusparse -localapi {{0 1}} -neg {{0 1}}
264:      output_file: output/ex123_2.out

266:    test:
267:      requires: kokkos_kernels !sycl
268:      suffix: 2_kokkos
269:      nsize: 7
270:      filter: grep -v type
271:      diff_args: -j
272:      args: -mat_type mpiaijkokkos -localapi {{0 1}} -neg {{0 1}}
273:      output_file: output/ex123_2.out

275:    test:
276:      suffix: 3
277:      nsize: 3
278:      filter: grep -v type
279:      diff_args: -j
280:      args: -mat_type mpiaij -loc -localapi {{0 1}} -neg {{0 1}}

282:    test:
283:      requires: cuda
284:      suffix: 3_cuda
285:      nsize: 3
286:      filter: grep -v type
287:      diff_args: -j
288:      args: -mat_type mpiaijcusparse -loc -localapi {{0 1}} -neg {{0 1}}
289:      output_file: output/ex123_3.out

291:    test:
292:      requires: !sycl kokkos_kernels
293:      suffix: 3_kokkos
294:      nsize: 3
295:      filter: grep -v type
296:      diff_args: -j
297:      args: -mat_type aijkokkos -loc -localapi {{0 1}} -neg {{0 1}}
298:      output_file: output/ex123_3.out

300:    test:
301:      suffix: 4
302:      nsize: 4
303:      filter: grep -v type
304:      diff_args: -j
305:      args: -mat_type mpiaij -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}}

307:    test:
308:      requires: cuda
309:      suffix: 4_cuda
310:      nsize: 4
311:      filter: grep -v type
312:      diff_args: -j
313:      args: -mat_type mpiaijcusparse -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}}
314:      output_file: output/ex123_4.out

316:    test:
317:      requires: !sycl kokkos_kernels
318:      suffix: 4_kokkos
319:      nsize: 4
320:      filter: grep -v type
321:      diff_args: -j
322:      args: -mat_type aijkokkos -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}}
323:      output_file: output/ex123_4.out

325:    test:
326:      suffix: matis
327:      nsize: 3
328:      filter: grep -v type
329:      diff_args: -j
330:      args: -mat_type is -localapi {{0 1}} -neg {{0 1}}

332: TEST*/