Actual source code: ex49.c
1: static char help[] = "Tests MatTranspose(), MatNorm(), and MatAXPY().\n\n";
3: #include <petscmat.h>
5: int main(int argc, char **argv)
6: {
7: Mat mat, tmat = 0;
8: PetscInt m = 4, n, i, j;
9: PetscMPIInt size, rank;
10: PetscInt rstart, rend, rect = 0;
11: PetscBool flg;
12: PetscScalar v;
13: PetscReal normf, normi, norm1;
14: MatInfo info;
16: PetscFunctionBeginUser;
17: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
18: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
19: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
20: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
21: n = m;
22: PetscCall(PetscOptionsHasName(NULL, NULL, "-rect1", &flg));
23: if (flg) {
24: n += 2;
25: rect = 1;
26: }
27: PetscCall(PetscOptionsHasName(NULL, NULL, "-rect2", &flg));
28: if (flg) {
29: n -= 2;
30: rect = 1;
31: }
33: /* Create and assemble matrix */
34: PetscCall(MatCreate(PETSC_COMM_WORLD, &mat));
35: PetscCall(MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, m, n));
36: PetscCall(MatSetFromOptions(mat));
37: PetscCall(MatSetUp(mat));
38: PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
39: for (i = rstart; i < rend; i++) {
40: for (j = 0; j < n; j++) {
41: v = 10 * i + j;
42: PetscCall(MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES));
43: }
44: }
45: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
46: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
48: /* Print info about original matrix */
49: PetscCall(MatGetInfo(mat, MAT_GLOBAL_SUM, &info));
50: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "original matrix nonzeros = %" PetscInt_FMT ", allocated nonzeros = %" PetscInt_FMT "\n", (PetscInt)info.nz_used, (PetscInt)info.nz_allocated));
51: PetscCall(MatNorm(mat, NORM_FROBENIUS, &normf));
52: PetscCall(MatNorm(mat, NORM_1, &norm1));
53: PetscCall(MatNorm(mat, NORM_INFINITY, &normi));
54: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "original: Frobenius norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi));
55: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
57: /* Form matrix transpose */
58: PetscCall(PetscOptionsHasName(NULL, NULL, "-in_place", &flg));
59: if (flg) {
60: PetscCall(MatTranspose(mat, MAT_INPLACE_MATRIX, &mat)); /* in-place transpose */
61: tmat = mat;
62: mat = 0;
63: } else { /* out-of-place transpose */
64: PetscCall(MatTranspose(mat, MAT_INITIAL_MATRIX, &tmat));
65: }
67: /* Print info about transpose matrix */
68: PetscCall(MatGetInfo(tmat, MAT_GLOBAL_SUM, &info));
69: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "transpose matrix nonzeros = %" PetscInt_FMT ", allocated nonzeros = %" PetscInt_FMT "\n", (PetscInt)info.nz_used, (PetscInt)info.nz_allocated));
70: PetscCall(MatNorm(tmat, NORM_FROBENIUS, &normf));
71: PetscCall(MatNorm(tmat, NORM_1, &norm1));
72: PetscCall(MatNorm(tmat, NORM_INFINITY, &normi));
73: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "transpose: Frobenius norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi));
74: PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
76: /* Test MatAXPY */
77: if (mat && !rect) {
78: PetscScalar alpha = 1.0;
79: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-alpha", &alpha, NULL));
80: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "matrix addition: B = B + alpha * A\n"));
81: PetscCall(MatAXPY(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN));
82: PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
83: }
85: /* Free data structures */
86: PetscCall(MatDestroy(&tmat));
87: if (mat) PetscCall(MatDestroy(&mat));
89: PetscCall(PetscFinalize());
90: return 0;
91: }
93: /*TEST
95: test:
97: testset:
98: args: -rect1
99: test:
100: suffix: r1
101: output_file: output/ex49_r1.out
102: test:
103: suffix: r1_inplace
104: args: -in_place
105: output_file: output/ex49_r1.out
106: test:
107: suffix: r1_par
108: nsize: 2
109: output_file: output/ex49_r1_par.out
110: test:
111: suffix: r1_par_inplace
112: args: -in_place
113: nsize: 2
114: output_file: output/ex49_r1_par.out
116: TEST*/