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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

218: /*TEST

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

379: TEST*/