Actual source code: ex2.c


  2: static char help[] = "Tests MatTranspose(), MatNorm(), MatAXPY() and MatAYPX().\n\n";

  4: #include <petscmat.h>

  6: static PetscErrorCode TransposeAXPY(Mat C, PetscScalar alpha, Mat mat, PetscErrorCode (*f)(Mat, Mat *))
  7: {
  8:   Mat     D, E, F, G;
  9:   MatType mtype;

 11:   MatGetType(mat, &mtype);
 12:   if (f == MatCreateTranspose) {
 13:     PetscPrintf(PETSC_COMM_WORLD, "\nMatAXPY:  (C^T)^T = (C^T)^T + alpha * A, C=A, SAME_NONZERO_PATTERN\n");
 14:   } else {
 15:     PetscPrintf(PETSC_COMM_WORLD, "\nMatAXPY:  (C^H)^H = (C^H)^H + alpha * A, C=A, SAME_NONZERO_PATTERN\n");
 16:   }
 17:   MatDuplicate(mat, MAT_COPY_VALUES, &C);
 18:   f(C, &D);
 19:   f(D, &E);
 20:   MatAXPY(E, alpha, mat, SAME_NONZERO_PATTERN);
 21:   MatConvert(E, mtype, MAT_INPLACE_MATRIX, &E);
 22:   MatView(E, PETSC_VIEWER_STDOUT_WORLD);
 23:   MatDestroy(&E);
 24:   MatDestroy(&D);
 25:   MatDestroy(&C);
 26:   if (f == MatCreateTranspose) {
 27:     PetscPrintf(PETSC_COMM_WORLD, "MatAXPY:  C = C + alpha * (A^T)^T, C=A, SAME_NONZERO_PATTERN\n");
 28:   } else {
 29:     PetscPrintf(PETSC_COMM_WORLD, "MatAXPY:  C = C + alpha * (A^H)^H, C=A, SAME_NONZERO_PATTERN\n");
 30:   }
 31:   if (f == MatCreateTranspose) {
 32:     MatTranspose(mat, MAT_INITIAL_MATRIX, &D);
 33:   } else {
 34:     MatHermitianTranspose(mat, MAT_INITIAL_MATRIX, &D);
 35:   }
 36:   f(D, &E);
 37:   MatDuplicate(mat, MAT_COPY_VALUES, &C);
 38:   MatAXPY(C, alpha, E, SAME_NONZERO_PATTERN);
 39:   MatView(C, PETSC_VIEWER_STDOUT_WORLD);
 40:   MatDestroy(&E);
 41:   MatDestroy(&D);
 42:   MatDestroy(&C);
 43:   if (f == MatCreateTranspose) {
 44:     PetscPrintf(PETSC_COMM_WORLD, "MatAXPY:  (C^T)^T = (C^T)^T + alpha * (A^T)^T, C=A, SAME_NONZERO_PATTERN\n");
 45:   } else {
 46:     PetscPrintf(PETSC_COMM_WORLD, "MatAXPY:  (C^H)^H = (C^H)^H + alpha * (A^H)^H, C=A, SAME_NONZERO_PATTERN\n");
 47:   }
 48:   MatDuplicate(mat, MAT_COPY_VALUES, &C);
 49:   f(C, &D);
 50:   f(D, &E);
 51:   f(mat, &F);
 52:   f(F, &G);
 53:   MatAXPY(E, alpha, G, SAME_NONZERO_PATTERN);
 54:   MatConvert(E, mtype, MAT_INPLACE_MATRIX, &E);
 55:   MatView(E, PETSC_VIEWER_STDOUT_WORLD);
 56:   MatDestroy(&G);
 57:   MatDestroy(&F);
 58:   MatDestroy(&E);
 59:   MatDestroy(&D);
 60:   MatDestroy(&C);

 62:   /* Call f on a matrix that does not implement the transposition */
 63:   PetscPrintf(PETSC_COMM_WORLD, "MatAXPY:  Now without the transposition operation\n");
 64:   MatConvert(mat, MATSHELL, MAT_INITIAL_MATRIX, &C);
 65:   f(C, &D);
 66:   f(D, &E);
 67:   /* XXX cannot use MAT_INPLACE_MATRIX, it leaks mat */
 68:   MatConvert(E, mtype, MAT_INITIAL_MATRIX, &F);
 69:   MatAXPY(F, alpha, mat, SAME_NONZERO_PATTERN);
 70:   MatView(F, PETSC_VIEWER_STDOUT_WORLD);
 71:   MatDestroy(&F);
 72:   MatDestroy(&E);
 73:   MatDestroy(&D);
 74:   MatDestroy(&C);
 75:   return 0;
 76: }

 78: int main(int argc, char **argv)
 79: {
 80:   Mat         mat, tmat = 0;
 81:   PetscInt    m = 7, n, i, j, rstart, rend, rect = 0;
 82:   PetscMPIInt size, rank;
 83:   PetscBool   flg;
 84:   PetscScalar v, alpha;
 85:   PetscReal   normf, normi, norm1;

 88:   PetscInitialize(&argc, &argv, (char *)0, help);
 89:   PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON);
 90:   PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
 91:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 92:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 93:   n = m;
 94:   PetscOptionsHasName(NULL, NULL, "-rectA", &flg);
 95:   if (flg) {
 96:     n += 2;
 97:     rect = 1;
 98:   }
 99:   PetscOptionsHasName(NULL, NULL, "-rectB", &flg);
100:   if (flg) {
101:     n -= 2;
102:     rect = 1;
103:   }

105:   /* ------- Assemble matrix --------- */
106:   MatCreate(PETSC_COMM_WORLD, &mat);
107:   MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, m, n);
108:   MatSetFromOptions(mat);
109:   MatSetUp(mat);
110:   MatGetOwnershipRange(mat, &rstart, &rend);
111:   for (i = rstart; i < rend; i++) {
112:     for (j = 0; j < n; j++) {
113:       v = 10.0 * i + j + 1.0;
114:       MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES);
115:     }
116:   }
117:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
118:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);

120:   /* ----------------- Test MatNorm()  ----------------- */
121:   MatNorm(mat, NORM_FROBENIUS, &normf);
122:   MatNorm(mat, NORM_1, &norm1);
123:   MatNorm(mat, NORM_INFINITY, &normi);
124:   PetscPrintf(PETSC_COMM_WORLD, "original A: Frobenious norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi);
125:   MatView(mat, PETSC_VIEWER_STDOUT_WORLD);

127:   /* --------------- Test MatTranspose()  -------------- */
128:   PetscOptionsHasName(NULL, NULL, "-in_place", &flg);
129:   if (!rect && flg) {
130:     MatTranspose(mat, MAT_REUSE_MATRIX, &mat); /* in-place transpose */
131:     tmat = mat;
132:     mat  = NULL;
133:   } else { /* out-of-place transpose */
134:     MatTranspose(mat, MAT_INITIAL_MATRIX, &tmat);
135:   }

137:   /* ----------------- Test MatNorm()  ----------------- */
138:   /* Print info about transpose matrix */
139:   MatNorm(tmat, NORM_FROBENIUS, &normf);
140:   MatNorm(tmat, NORM_1, &norm1);
141:   MatNorm(tmat, NORM_INFINITY, &normi);
142:   PetscPrintf(PETSC_COMM_WORLD, "B = A^T: Frobenious norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi);
143:   MatView(tmat, PETSC_VIEWER_STDOUT_WORLD);

145:   /* ----------------- Test MatAXPY(), MatAYPX()  ----------------- */
146:   if (mat && !rect) {
147:     alpha = 1.0;
148:     PetscOptionsGetScalar(NULL, NULL, "-alpha", &alpha, NULL);
149:     PetscPrintf(PETSC_COMM_WORLD, "MatAXPY:  B = B + alpha * A\n");
150:     MatAXPY(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN);
151:     MatView(tmat, PETSC_VIEWER_STDOUT_WORLD);

153:     PetscPrintf(PETSC_COMM_WORLD, "MatAYPX:  B = alpha*B + A\n");
154:     MatAYPX(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN);
155:     MatView(tmat, PETSC_VIEWER_STDOUT_WORLD);
156:   }

158:   {
159:     Mat C;
160:     alpha = 1.0;
161:     PetscPrintf(PETSC_COMM_WORLD, "MatAXPY:  C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n");
162:     MatDuplicate(mat, MAT_COPY_VALUES, &C);
163:     MatAXPY(C, alpha, mat, SAME_NONZERO_PATTERN);
164:     MatView(C, PETSC_VIEWER_STDOUT_WORLD);
165:     MatDestroy(&C);
166:     TransposeAXPY(C, alpha, mat, MatCreateTranspose);
167:     TransposeAXPY(C, alpha, mat, MatCreateHermitianTranspose);
168:   }

170:   {
171:     Mat matB;
172:     /* get matB that has nonzeros of mat in all even numbers of row and col */
173:     MatCreate(PETSC_COMM_WORLD, &matB);
174:     MatSetSizes(matB, PETSC_DECIDE, PETSC_DECIDE, m, n);
175:     MatSetFromOptions(matB);
176:     MatSetUp(matB);
177:     MatGetOwnershipRange(matB, &rstart, &rend);
178:     if (rstart % 2 != 0) rstart++;
179:     for (i = rstart; i < rend; i += 2) {
180:       for (j = 0; j < n; j += 2) {
181:         v = 10.0 * i + j + 1.0;
182:         MatSetValues(matB, 1, &i, 1, &j, &v, INSERT_VALUES);
183:       }
184:     }
185:     MatAssemblyBegin(matB, MAT_FINAL_ASSEMBLY);
186:     MatAssemblyEnd(matB, MAT_FINAL_ASSEMBLY);
187:     PetscPrintf(PETSC_COMM_WORLD, " A: original matrix:\n");
188:     MatView(mat, PETSC_VIEWER_STDOUT_WORLD);
189:     PetscPrintf(PETSC_COMM_WORLD, " B(a subset of A):\n");
190:     MatView(matB, PETSC_VIEWER_STDOUT_WORLD);
191:     PetscPrintf(PETSC_COMM_WORLD, "MatAXPY:  B = B + alpha * A, SUBSET_NONZERO_PATTERN\n");
192:     MatAXPY(mat, alpha, matB, SUBSET_NONZERO_PATTERN);
193:     MatView(mat, PETSC_VIEWER_STDOUT_WORLD);
194:     MatDestroy(&matB);
195:   }

197:   /* Test MatZeroRows */
198:   j = rstart - 1;
199:   if (j < 0) j = m - 1;
200:   PetscPrintf(PETSC_COMM_WORLD, "MatZeroRows:\n");
201:   MatZeroRows(mat, 1, &j, 0.0, NULL, NULL);
202:   MatView(mat, PETSC_VIEWER_STDOUT_WORLD);

204:   /* Test MatShift */
205:   PetscPrintf(PETSC_COMM_WORLD, "MatShift: B = B - 2*I\n");
206:   MatShift(mat, -2.0);
207:   MatView(mat, PETSC_VIEWER_STDOUT_WORLD);

209:   PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
210:   /* Free data structures */
211:   MatDestroy(&mat);
212:   MatDestroy(&tmat);
213:   PetscFinalize();
214:   return 0;
215: }

217: /*TEST

219:    test:
220:       suffix: 11_A
221:       args: -mat_type seqaij -rectA
222:       filter: grep -v "Mat Object"

224:    test:
225:       suffix: 12_A
226:       args: -mat_type seqdense -rectA
227:       filter: grep -v type | grep -v "Mat Object"

229:    test:
230:       requires: cuda
231:       suffix: 12_A_cuda
232:       args: -mat_type seqdensecuda -rectA
233:       output_file: output/ex2_12_A.out
234:       filter: grep -v type | grep -v "Mat Object"

236:    test:
237:       requires: kokkos_kernels
238:       suffix: 12_A_kokkos
239:       args: -mat_type aijkokkos -rectA
240:       output_file: output/ex2_12_A.out
241:       filter: grep -v type | grep -v "Mat Object"

243:    test:
244:       suffix: 11_B
245:       args: -mat_type seqaij -rectB
246:       filter: grep -v "Mat Object"

248:    test:
249:       suffix: 12_B
250:       args: -mat_type seqdense -rectB
251:       filter: grep -v type | grep -v "Mat Object"

253:    testset:
254:       args: -rectB
255:       output_file: output/ex2_12_B.out
256:       filter: grep -v type | grep -v "Mat Object"

258:       test:
259:          requires: cuda
260:          suffix: 12_B_cuda
261:          args: -mat_type {{seqdensecuda seqaijcusparse}}

263:       test:
264:          requires: kokkos_kernels
265:          suffix: 12_B_kokkos
266:          args: -mat_type aijkokkos

268:       test:
269:          suffix: 12_B_aij
270:          args: -mat_type aij
271:    test:
272:       suffix: 21
273:       args: -mat_type mpiaij
274:       filter: grep -v type | grep -v " MPI process"

276:    test:
277:       suffix: 22
278:       args: -mat_type mpidense
279:       filter: grep -v type | grep -v "Mat Object"

281:    test:
282:       requires: cuda
283:       suffix: 22_cuda
284:       output_file: output/ex2_22.out
285:       args: -mat_type mpidensecuda
286:       filter: grep -v type | grep -v "Mat Object"

288:    test:
289:       requires: kokkos_kernels
290:       suffix: 22_kokkos
291:       output_file: output/ex2_22.out
292:       args: -mat_type aijkokkos
293:       filter: grep -v type | grep -v "Mat Object"

295:    test:
296:       suffix: 23
297:       nsize: 3
298:       args: -mat_type mpiaij
299:       filter: grep -v type | grep -v " MPI process"

301:    test:
302:       suffix: 24
303:       nsize: 3
304:       args: -mat_type mpidense
305:       filter: grep -v type | grep -v "Mat Object"

307:    test:
308:       requires: cuda
309:       suffix: 24_cuda
310:       nsize: 3
311:       output_file: output/ex2_24.out
312:       args: -mat_type mpidensecuda
313:       filter: grep -v type | grep -v "Mat Object"

315:    test:
316:       suffix: 2_aijcusparse_1
317:       args: -mat_type mpiaijcusparse
318:       output_file: output/ex2_21.out
319:       requires: cuda
320:       filter: grep -v type | grep -v " MPI process"

322:    test:
323:       suffix: 2_aijkokkos_1
324:       args: -mat_type aijkokkos
325:       output_file: output/ex2_21.out
326:       requires: kokkos_kernels
327:       filter: grep -v type | grep -v " MPI process"

329:    test:
330:       suffix: 2_aijcusparse_2
331:       nsize: 3
332:       args: -mat_type mpiaijcusparse
333:       output_file: output/ex2_23.out
334:       requires: cuda
335:       filter: grep -v type | grep -v " MPI process"

337:    test:
338:       suffix: 2_aijkokkos_2
339:       nsize: 3
340:       args: -mat_type aijkokkos
341:       output_file: output/ex2_23.out
342:       # Turn off hip due to intermittent CI failures on hip.txcorp.com. Should re-enable this test when the machine is upgraded.
343:       requires: !hip kokkos_kernels
344:       filter: grep -v type | grep -v "MPI processes"

346:    test:
347:       suffix: 3
348:       nsize: 2
349:       args: -mat_type mpiaij -rectA

351:    test:
352:       suffix: 3_aijcusparse
353:       nsize: 2
354:       args: -mat_type mpiaijcusparse -rectA
355:       requires: cuda

357:    test:
358:       suffix: 4
359:       nsize: 2
360:       args: -mat_type mpidense -rectA
361:       filter: grep -v type | grep -v " MPI process"

363:    test:
364:       requires: cuda
365:       suffix: 4_cuda
366:       nsize: 2
367:       output_file: output/ex2_4.out
368:       args: -mat_type mpidensecuda -rectA
369:       filter: grep -v type | grep -v " MPI process"

371:    test:
372:       suffix: aijcusparse_1
373:       args: -mat_type seqaijcusparse -rectA
374:       filter: grep -v "Mat Object"
375:       output_file: output/ex2_11_A_aijcusparse.out
376:       requires: cuda

378: TEST*/