Actual source code: ex2.c

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

  3: #include <petscmat.h>

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

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

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

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

105:   /* ------- Assemble matrix --------- */
106:   PetscCall(MatCreate(PETSC_COMM_WORLD, &mat));
107:   PetscCall(MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, m, n));
108:   PetscCall(MatSetFromOptions(mat));
109:   PetscCall(MatSetUp(mat));
110:   PetscCall(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:       PetscCall(MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES));
115:     }
116:   }
117:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
118:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));

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

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

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

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

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

158:   {
159:     Mat C;
160:     alpha = 1.0;
161:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY:  C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n"));
162:     PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
163:     PetscCall(MatAXPY(C, alpha, mat, SAME_NONZERO_PATTERN));
164:     PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
165:     PetscCall(MatDestroy(&C));
166:     PetscCall(TransposeAXPY(C, alpha, mat, MatCreateTranspose));
167:     PetscCall(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:     PetscCall(MatCreate(PETSC_COMM_WORLD, &matB));
174:     PetscCall(MatSetSizes(matB, PETSC_DECIDE, PETSC_DECIDE, m, n));
175:     PetscCall(MatSetFromOptions(matB));
176:     PetscCall(MatSetUp(matB));
177:     PetscCall(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:         PetscCall(MatSetValues(matB, 1, &i, 1, &j, &v, INSERT_VALUES));
183:       }
184:     }
185:     PetscCall(MatAssemblyBegin(matB, MAT_FINAL_ASSEMBLY));
186:     PetscCall(MatAssemblyEnd(matB, MAT_FINAL_ASSEMBLY));
187:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " A: original matrix:\n"));
188:     PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
189:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " B(a subset of A):\n"));
190:     PetscCall(MatView(matB, PETSC_VIEWER_STDOUT_WORLD));
191:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY:  B = B + alpha * A, SUBSET_NONZERO_PATTERN\n"));
192:     PetscCall(MatAXPY(mat, alpha, matB, SUBSET_NONZERO_PATTERN));
193:     PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
194:     PetscCall(MatDestroy(&matB));
195:   }

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

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

209:   PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
210:   /* Free data structures */
211:   PetscCall(MatDestroy(&mat));
212:   PetscCall(MatDestroy(&tmat));
213:   PetscCall(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*/