1 static char help[] = "Tests Elemental interface.\n\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Mat Cdense, B, C, Ct;
8 Vec d;
9 PetscInt i, j, m = 5, n, nrows, ncols;
10 const PetscInt *rows, *cols;
11 IS isrows, iscols;
12 PetscScalar *v;
13 PetscMPIInt rank, size;
14 PetscReal Cnorm;
15 PetscBool flg, mats_view = PETSC_FALSE;
16
17 PetscFunctionBeginUser;
18 PetscCall(PetscInitialize(&argc, &args, NULL, help));
19 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
20 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
21 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
22 n = m;
23 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
24 PetscCall(PetscOptionsHasName(NULL, NULL, "-mats_view", &mats_view));
25
26 PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
27 PetscCall(MatSetSizes(C, m, n, PETSC_DECIDE, PETSC_DECIDE));
28 PetscCall(MatSetType(C, MATELEMENTAL));
29 PetscCall(MatSetFromOptions(C));
30 PetscCall(MatSetUp(C));
31 PetscCall(MatGetOwnershipIS(C, &isrows, &iscols));
32 PetscCall(ISGetLocalSize(isrows, &nrows));
33 PetscCall(ISGetIndices(isrows, &rows));
34 PetscCall(ISGetLocalSize(iscols, &ncols));
35 PetscCall(ISGetIndices(iscols, &cols));
36 PetscCall(PetscMalloc1(nrows * ncols, &v));
37 #if defined(PETSC_USE_COMPLEX)
38 PetscRandom rand;
39 PetscScalar rval;
40 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
41 PetscCall(PetscRandomSetFromOptions(rand));
42 for (i = 0; i < nrows; i++) {
43 for (j = 0; j < ncols; j++) {
44 PetscCall(PetscRandomGetValue(rand, &rval));
45 v[i * ncols + j] = rval;
46 }
47 }
48 PetscCall(PetscRandomDestroy(&rand));
49 #else
50 for (i = 0; i < nrows; i++) {
51 for (j = 0; j < ncols; j++) v[i * ncols + j] = (PetscReal)(10000 * rank + 100 * rows[i] + cols[j]);
52 }
53 #endif
54 PetscCall(MatSetValues(C, nrows, rows, ncols, cols, v, INSERT_VALUES));
55 PetscCall(ISRestoreIndices(isrows, &rows));
56 PetscCall(ISRestoreIndices(iscols, &cols));
57 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
58 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
59 PetscCall(ISDestroy(&isrows));
60 PetscCall(ISDestroy(&iscols));
61
62 /* Test MatView(), MatDuplicate() and out-of-place MatConvert() */
63 PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &B));
64 if (mats_view) {
65 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Duplicated C:\n"));
66 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
67 }
68 PetscCall(MatDestroy(&B));
69 PetscCall(MatConvert(C, MATMPIDENSE, MAT_INITIAL_MATRIX, &Cdense));
70 PetscCall(MatMultEqual(C, Cdense, 5, &flg));
71 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "Cdense != C. MatConvert() fails");
72
73 /* Test MatNorm() */
74 PetscCall(MatNorm(C, NORM_1, &Cnorm));
75
76 /* Test MatTranspose(), MatZeroEntries() and MatGetDiagonal() */
77 PetscCall(MatTranspose(C, MAT_INITIAL_MATRIX, &Ct));
78 PetscCall(MatConjugate(Ct));
79 if (mats_view) {
80 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C's Transpose Conjugate:\n"));
81 PetscCall(MatView(Ct, PETSC_VIEWER_STDOUT_WORLD));
82 }
83 PetscCall(MatZeroEntries(Ct));
84 PetscCall(VecCreate(PETSC_COMM_WORLD, &d));
85 PetscCall(VecSetSizes(d, m > n ? n : m, PETSC_DECIDE));
86 PetscCall(VecSetFromOptions(d));
87 PetscCall(MatGetDiagonal(C, d));
88 if (mats_view) {
89 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Diagonal of C:\n"));
90 PetscCall(VecView(d, PETSC_VIEWER_STDOUT_WORLD));
91 }
92 if (m > n) {
93 PetscCall(MatDiagonalScale(C, NULL, d));
94 } else {
95 PetscCall(MatDiagonalScale(C, d, NULL));
96 }
97 if (mats_view) {
98 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Diagonal Scaled C:\n"));
99 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
100 }
101
102 /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */
103 PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
104 PetscCall(MatSetSizes(B, m, n, PETSC_DECIDE, PETSC_DECIDE));
105 PetscCall(MatSetType(B, MATELEMENTAL));
106 PetscCall(MatSetFromOptions(B));
107 PetscCall(MatSetUp(B));
108 PetscCall(MatGetOwnershipIS(B, &isrows, &iscols));
109 PetscCall(ISGetLocalSize(isrows, &nrows));
110 PetscCall(ISGetIndices(isrows, &rows));
111 PetscCall(ISGetLocalSize(iscols, &ncols));
112 PetscCall(ISGetIndices(iscols, &cols));
113 for (i = 0; i < nrows; i++) {
114 for (j = 0; j < ncols; j++) v[i * ncols + j] = (PetscReal)(1000 * rows[i] + cols[j]);
115 }
116 PetscCall(MatSetValues(B, nrows, rows, ncols, cols, v, INSERT_VALUES));
117 PetscCall(PetscFree(v));
118 PetscCall(ISRestoreIndices(isrows, &rows));
119 PetscCall(ISRestoreIndices(iscols, &cols));
120 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
121 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
122 PetscCall(MatAXPY(B, 2.5, C, SAME_NONZERO_PATTERN));
123 PetscCall(MatAYPX(B, 3.75, C, SAME_NONZERO_PATTERN));
124 PetscCall(MatConvert(B, MATDENSE, MAT_INPLACE_MATRIX, &B));
125 if (mats_view) {
126 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B after MatAXPY and MatAYPX:\n"));
127 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
128 }
129 PetscCall(ISDestroy(&isrows));
130 PetscCall(ISDestroy(&iscols));
131 PetscCall(MatDestroy(&B));
132
133 /* Test MatMatTransposeMult(): B = C*C^T */
134 PetscCall(MatMatTransposeMult(C, C, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &B));
135 PetscCall(MatScale(C, 2.0));
136 PetscCall(MatMatTransposeMult(C, C, MAT_REUSE_MATRIX, PETSC_DETERMINE, &B));
137 PetscCall(MatMatTransposeMultEqual(C, C, B, 10, &flg));
138 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatMatTransposeMult: B != C*B^T");
139
140 if (mats_view) {
141 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C MatMatTransposeMult C:\n"));
142 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
143 }
144
145 PetscCall(MatDestroy(&Cdense));
146 PetscCall(PetscFree(v));
147 PetscCall(MatDestroy(&B));
148 PetscCall(MatDestroy(&C));
149 PetscCall(MatDestroy(&Ct));
150 PetscCall(VecDestroy(&d));
151 PetscCall(PetscFinalize());
152 return 0;
153 }
154
155 /*TEST
156
157 test:
158 nsize: 2
159 args: -m 3 -n 2
160 requires: elemental
161 output_file: output/empty.out
162
163 test:
164 suffix: 2
165 nsize: 6
166 args: -m 2 -n 3
167 requires: elemental
168 output_file: output/empty.out
169
170 test:
171 suffix: 3
172 nsize: 1
173 args: -m 2 -n 3
174 requires: elemental
175 output_file: output/empty.out
176
177 TEST*/
178