1 static char help[] = "Tests MatTranspose(), MatNorm(), and MatAXPY().\n\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** argv)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;
15
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 }
32
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));
47
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));
56
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 }
66
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));
75
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 }
84
85 /* Free data structures */
86 PetscCall(MatDestroy(&tmat));
87 if (mat) PetscCall(MatDestroy(&mat));
88
89 PetscCall(PetscFinalize());
90 return 0;
91 }
92
93 /*TEST
94
95 test:
96
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
115
116 TEST*/
117