Actual source code: ex123.c

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

  3: #include <petscmat.h>

  5: int main(int argc, char **args)
  6: {
  7:   Mat                    A, At, AAt, T = NULL;
  8:   Vec                    x, y, z;
  9:   ISLocalToGlobalMapping rl2g, cl2g;
 10:   IS                     is;
 11:   PetscLayout            rmap, cmap;
 12:   PetscInt              *it, *jt;
 13:   PetscInt               n1 = 11, n2 = 9;
 14:   PetscInt               i1[] = {7, 6, 2, 0, 4, 1, 1, 0, 2, 2, 1, -1, -1};
 15:   PetscInt               j1[] = {1, 4, 3, 5, 3, 3, 4, 5, 0, 3, 1, -1, -1};
 16:   PetscInt               i2[] = {7, 6, 2, 0, 4, 1, 1, 2, 1, -1, -1};
 17:   PetscInt               j2[] = {1, 4, 3, 5, 3, 3, 4, 0, 1, -1, -1};
 18:   PetscScalar            v1[] = {-1., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., PETSC_MAX_REAL, PETSC_MAX_REAL};
 19:   PetscScalar            v2[] = {1., -1., -2., -3., -4., -5., -6., -7., -8., -9., -10., PETSC_MAX_REAL, PETSC_MAX_REAL};
 20:   PetscInt               N = 6, m = 8, M, rstart, cstart, i;
 21:   PetscMPIInt            size;
 22:   PetscBool              loc      = PETSC_FALSE;
 23:   PetscBool              locdiag  = PETSC_TRUE;
 24:   PetscBool              localapi = PETSC_FALSE;
 25:   PetscBool              neg      = PETSC_FALSE;
 26:   PetscBool              ismatis, ismpiaij, ishypre;

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

 53:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &ismatis));
 54:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre));

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

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

120:   /* INSERT_VALUES will overwrite matrix entries but
121:      still perform the sum of the repeated entries */
122:   PetscCall(MatSetValuesCOO(A, v2, INSERT_VALUES));
123:   PetscCall(MatView(A, NULL));

125:   /* test with unique entries */
126:   PetscCall(PetscArraycpy(it, i2, n2));
127:   PetscCall(PetscArraycpy(jt, j2, n2));
128:   if (!localapi) {
129:     PetscCall(MatSetPreallocationCOO(A, n2, it, jt));
130:   } else {
131:     PetscCall(MatSetPreallocationCOOLocal(A, n2, it, jt));
132:   }
133:   PetscCall(MatSetValuesCOO(A, v1, ADD_VALUES));
134:   PetscCall(MatMult(A, x, y));
135:   PetscCall(MatView(A, NULL));
136:   PetscCall(VecView(y, NULL));
137:   PetscCall(MatSetValuesCOO(A, v2, ADD_VALUES));
138:   PetscCall(MatMultAdd(A, x, y, z));
139:   PetscCall(MatView(A, NULL));
140:   PetscCall(VecView(z, NULL));
141:   PetscCall(PetscArraycpy(it, i2, n2));
142:   PetscCall(PetscArraycpy(jt, j2, n2));
143:   if (!localapi) {
144:     PetscCall(MatSetPreallocationCOO(A, n2, it, jt));
145:   } else {
146:     PetscCall(MatSetPreallocationCOOLocal(A, n2, it, jt));
147:   }
148:   PetscCall(MatSetValuesCOO(A, v1, INSERT_VALUES));
149:   PetscCall(MatMult(A, x, y));
150:   PetscCall(MatView(A, NULL));
151:   PetscCall(VecView(y, NULL));
152:   PetscCall(MatSetValuesCOO(A, v2, INSERT_VALUES));
153:   PetscCall(MatMultAdd(A, x, y, z));
154:   PetscCall(MatView(A, NULL));
155:   PetscCall(VecView(z, NULL));
156:   T = A;
157:   if (ishypre) PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &T));
158:   PetscCall(MatTranspose(T, MAT_INITIAL_MATRIX, &At));
159:   if (!ismatis) {
160:     PetscCall(MatMatMult(T, At, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt));
161:     PetscCall(MatView(AAt, NULL));
162:     PetscCall(MatDestroy(&AAt));
163:     PetscCall(MatMatMult(At, T, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt));
164:     PetscCall(MatView(AAt, NULL));
165:     PetscCall(MatDestroy(&AAt));
166:   }
167:   PetscCall(MatDestroy(&At));
168:   if (ishypre) PetscCall(MatDestroy(&T));

170:   /* test providing diagonal first, then off-diagonal */
171:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
172:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
173:   if ((ismpiaij || ishypre) && size > 1) {
174:     Mat                lA, lB;
175:     const PetscInt    *garray, *iA, *jA, *iB, *jB;
176:     const PetscScalar *vA, *vB;
177:     PetscScalar       *coo_v;
178:     PetscInt          *coo_i, *coo_j;
179:     PetscInt           i, j, nA, nB, nnz;
180:     PetscBool          flg;

182:     T = A;
183:     if (ishypre) PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &T));
184:     PetscCall(MatMPIAIJGetSeqAIJ(T, &lA, &lB, &garray));
185:     PetscCall(MatSeqAIJGetArrayRead(lA, &vA));
186:     PetscCall(MatSeqAIJGetArrayRead(lB, &vB));
187:     PetscCall(MatGetRowIJ(lA, 0, PETSC_FALSE, PETSC_FALSE, &nA, &iA, &jA, &flg));
188:     PetscCall(MatGetRowIJ(lB, 0, PETSC_FALSE, PETSC_FALSE, &nB, &iB, &jB, &flg));
189:     nnz = iA[nA] + iB[nB];
190:     PetscCall(PetscMalloc3(nnz, &coo_i, nnz, &coo_j, nnz, &coo_v));
191:     nnz = 0;
192:     for (i = 0; i < nA; i++) {
193:       for (j = iA[i]; j < iA[i + 1]; j++, nnz++) {
194:         coo_i[nnz] = i + rstart;
195:         coo_j[nnz] = jA[j] + cstart;
196:         coo_v[nnz] = vA[j];
197:       }
198:     }
199:     for (i = 0; i < nB; i++) {
200:       for (j = iB[i]; j < iB[i + 1]; j++, nnz++) {
201:         coo_i[nnz] = i + rstart;
202:         coo_j[nnz] = garray[jB[j]];
203:         coo_v[nnz] = vB[j];
204:       }
205:     }
206:     PetscCall(MatRestoreRowIJ(lA, 0, PETSC_FALSE, PETSC_FALSE, &nA, &iA, &jA, &flg));
207:     PetscCall(MatRestoreRowIJ(lB, 0, PETSC_FALSE, PETSC_FALSE, &nB, &iB, &jB, &flg));
208:     PetscCall(MatSeqAIJRestoreArrayRead(lA, &vA));
209:     PetscCall(MatSeqAIJRestoreArrayRead(lB, &vB));
210:     if (ishypre) PetscCall(MatDestroy(&T));

212:     PetscCall(MatSetPreallocationCOO(A, nnz, coo_i, coo_j));
213:     PetscCall(MatSetValuesCOO(A, coo_v, ADD_VALUES));
214:     PetscCall(MatMult(A, x, y));
215:     PetscCall(MatView(A, NULL));
216:     PetscCall(VecView(y, NULL));
217:     PetscCall(MatSetValuesCOO(A, coo_v, INSERT_VALUES));
218:     PetscCall(MatMult(A, x, y));
219:     PetscCall(MatView(A, NULL));
220:     PetscCall(VecView(y, NULL));

222:     T = A;
223:     if (ishypre) PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &T));
224:     PetscCall(MatTranspose(T, MAT_INITIAL_MATRIX, &At));
225:     PetscCall(MatMatMult(T, At, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt));
226:     PetscCall(MatView(AAt, NULL));
227:     PetscCall(MatDestroy(&AAt));
228:     PetscCall(MatMatMult(At, T, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt));
229:     PetscCall(MatView(AAt, NULL));
230:     PetscCall(MatDestroy(&AAt));
231:     PetscCall(MatDestroy(&At));
232:     if (ishypre) PetscCall(MatDestroy(&T));

234:     PetscCall(PetscFree3(coo_i, coo_j, coo_v));
235:   }
236:   PetscCall(PetscFree2(it, jt));
237:   PetscCall(VecDestroy(&z));
238:   PetscCall(VecDestroy(&x));
239:   PetscCall(VecDestroy(&y));
240:   PetscCall(MatDestroy(&A));
241:   PetscCall(PetscFinalize());
242:   return 0;
243: }

245: /*TEST

247:    test:
248:      suffix: 1
249:      filter: grep -v type | grep -v "Mat Object"
250:      diff_args: -j
251:      args: -mat_type {{seqaij mpiaij}} -localapi {{0 1}} -neg {{0 1}}

253:    test:
254:      requires: hypre
255:      suffix: 1_hypre
256:      filter: grep -v type | grep -v "Mat Object"
257:      diff_args: -j
258:      args: -mat_type hypre -localapi {{0 1}} -neg {{0 1}}
259:      output_file: output/ex123_1.out

261:    test:
262:      requires: cuda
263:      suffix: 1_cuda
264:      filter: grep -v type | grep -v "Mat Object"
265:      diff_args: -j
266:      args: -mat_type {{seqaijcusparse mpiaijcusparse}} -localapi {{0 1}} -neg {{0 1}}
267:      output_file: output/ex123_1.out

269:    test:
270:      requires: kokkos_kernels
271:      suffix: 1_kokkos
272:      filter: grep -v type | grep -v "Mat Object"
273:      diff_args: -j
274:      args: -mat_type {{seqaijkokkos mpiaijkokkos}} -localapi {{0 1}} -neg {{0 1}}
275:      output_file: output/ex123_1.out

277:    test:
278:      suffix: 2
279:      nsize: 7
280:      filter: grep -v type | grep -v "Mat Object"
281:      diff_args: -j
282:      args: -mat_type mpiaij -localapi {{0 1}} -neg {{0 1}}

284:    test:
285:      requires: hypre
286:      suffix: 2_hypre
287:      nsize: 7
288:      filter: grep -v type | grep -v "Mat Object"
289:      diff_args: -j
290:      args: -mat_type hypre -localapi {{0 1}} -neg {{0 1}}
291:      output_file: output/ex123_2.out

293:    test:
294:      requires: cuda
295:      suffix: 2_cuda
296:      nsize: 7
297:      filter: grep -v type | grep -v "Mat Object"
298:      diff_args: -j
299:      args: -mat_type mpiaijcusparse -localapi {{0 1}} -neg {{0 1}}
300:      output_file: output/ex123_2.out

302:    test:
303:      requires: kokkos_kernels
304:      suffix: 2_kokkos
305:      nsize: 7
306:      filter: grep -v type | grep -v "Mat Object"
307:      diff_args: -j
308:      args: -mat_type mpiaijkokkos -localapi {{0 1}} -neg {{0 1}}
309:      output_file: output/ex123_2.out

311:    test:
312:      suffix: 3
313:      nsize: 3
314:      filter: grep -v type | grep -v "Mat Object"
315:      diff_args: -j
316:      args: -mat_type mpiaij -loc -localapi {{0 1}} -neg {{0 1}}

318:    test:
319:      requires: hypre
320:      suffix: 3_hypre
321:      nsize: 3
322:      filter: grep -v type | grep -v "Mat Object"
323:      diff_args: -j
324:      args: -mat_type hypre -loc -localapi {{0 1}} -neg {{0 1}}
325:      output_file: output/ex123_3.out

327:    test:
328:      requires: cuda
329:      suffix: 3_cuda
330:      nsize: 3
331:      filter: grep -v type | grep -v "Mat Object"
332:      diff_args: -j
333:      args: -mat_type mpiaijcusparse -loc -localapi {{0 1}} -neg {{0 1}}
334:      output_file: output/ex123_3.out

336:    test:
337:      requires: kokkos_kernels
338:      suffix: 3_kokkos
339:      nsize: 3
340:      filter: grep -v type | grep -v "Mat Object"
341:      diff_args: -j
342:      args: -mat_type aijkokkos -loc -localapi {{0 1}} -neg {{0 1}}
343:      output_file: output/ex123_3.out

345:    test:
346:      suffix: 4
347:      nsize: 4
348:      filter: grep -v type | grep -v "Mat Object"
349:      diff_args: -j
350:      args: -mat_type mpiaij -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}}

352:    test:
353:      requires: hypre
354:      suffix: 4_hypre
355:      nsize: 4
356:      filter: grep -v type | grep -v "Mat Object"
357:      diff_args: -j
358:      args: -mat_type hypre -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}}
359:      output_file: output/ex123_4.out

361:    test:
362:      requires: cuda
363:      suffix: 4_cuda
364:      nsize: 4
365:      filter: grep -v type | grep -v "Mat Object"
366:      diff_args: -j
367:      args: -mat_type mpiaijcusparse -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}}
368:      output_file: output/ex123_4.out

370:    test:
371:      requires: kokkos_kernels
372:      suffix: 4_kokkos
373:      nsize: 4
374:      filter: grep -v type | grep -v "Mat Object"
375:      diff_args: -j
376:      args: -mat_type aijkokkos -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}}
377:      output_file: output/ex123_4.out

379:    test:
380:      suffix: matis
381:      nsize: 3
382:      filter: grep -v type | grep -v "Mat Object"
383:      diff_args: -j
384:      args: -mat_type is -localapi {{0 1}} -neg {{0 1}}

386: TEST*/