1 static char help[] = "Tests ScaLAPACK 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, Caij, B, C, Ct, Asub;
8 Vec d;
9 PetscInt i, j, M = 5, N, mb = 2, nb, nrows, ncols, mloc, nloc;
10 const PetscInt *rows, *cols;
11 IS isrows, iscols;
12 PetscScalar *v;
13 PetscMPIInt rank, color;
14 PetscReal Cnorm;
15 PetscBool flg, mats_view = PETSC_FALSE;
16 MPI_Comm subcomm;
17
18 PetscFunctionBeginUser;
19 PetscCall(PetscInitialize(&argc, &args, NULL, help));
20 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
21 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
22 N = M;
23 PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
24 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mb", &mb, NULL));
25 nb = mb;
26 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nb", &nb, NULL));
27 PetscCall(PetscOptionsHasName(NULL, NULL, "-mats_view", &mats_view));
28
29 PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
30 PetscCall(MatSetType(C, MATSCALAPACK));
31 mloc = PETSC_DECIDE;
32 PetscCall(PetscSplitOwnershipEqual(PETSC_COMM_WORLD, &mloc, &M));
33 nloc = PETSC_DECIDE;
34 PetscCall(PetscSplitOwnershipEqual(PETSC_COMM_WORLD, &nloc, &N));
35 PetscCall(MatSetSizes(C, mloc, nloc, M, N));
36 PetscCall(MatScaLAPACKSetBlockSizes(C, mb, nb));
37 PetscCall(MatSetFromOptions(C));
38 PetscCall(MatSetUp(C));
39 /*PetscCall(MatCreateScaLAPACK(PETSC_COMM_WORLD,mb,nb,M,N,0,0,&C)); */
40
41 PetscCall(MatGetOwnershipIS(C, &isrows, &iscols));
42 PetscCall(ISGetLocalSize(isrows, &nrows));
43 PetscCall(ISGetIndices(isrows, &rows));
44 PetscCall(ISGetLocalSize(iscols, &ncols));
45 PetscCall(ISGetIndices(iscols, &cols));
46 PetscCall(PetscMalloc1(nrows * ncols, &v));
47 for (i = 0; i < nrows; i++) {
48 for (j = 0; j < ncols; j++) v[i * ncols + j] = (PetscReal)(rows[i] + 1 + (cols[j] + 1) * 0.01);
49 }
50 PetscCall(MatSetValues(C, nrows, rows, ncols, cols, v, INSERT_VALUES));
51 PetscCall(PetscFree(v));
52 PetscCall(ISRestoreIndices(isrows, &rows));
53 PetscCall(ISRestoreIndices(iscols, &cols));
54 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
55 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
56 PetscCall(ISDestroy(&isrows));
57 PetscCall(ISDestroy(&iscols));
58
59 /* Test MatView(), MatDuplicate() and out-of-place MatConvert() */
60 PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &B));
61 if (mats_view) {
62 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Duplicated C:\n"));
63 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
64 }
65 PetscCall(MatDestroy(&B));
66 PetscCall(MatConvert(C, MATDENSE, MAT_INITIAL_MATRIX, &Cdense));
67 PetscCall(MatMultEqual(C, Cdense, 5, &flg));
68 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Check fails: Cdense != C");
69
70 /* Test MatNorm() */
71 PetscCall(MatNorm(C, NORM_1, &Cnorm));
72
73 /* Test MatTranspose(), MatZeroEntries() and MatGetDiagonal() */
74 PetscCall(MatTranspose(C, MAT_INITIAL_MATRIX, &Ct));
75 PetscCall(MatConjugate(Ct));
76 if (mats_view) {
77 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C's Transpose Conjugate:\n"));
78 PetscCall(MatView(Ct, PETSC_VIEWER_STDOUT_WORLD));
79 }
80 PetscCall(MatZeroEntries(Ct));
81 if (M > N) PetscCall(MatCreateVecs(C, &d, NULL));
82 else PetscCall(MatCreateVecs(C, NULL, &d));
83 PetscCall(MatGetDiagonal(C, d));
84 if (mats_view) {
85 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Diagonal of C:\n"));
86 PetscCall(VecView(d, PETSC_VIEWER_STDOUT_WORLD));
87 }
88 if (M > N) {
89 PetscCall(MatDiagonalScale(C, NULL, d));
90 } else {
91 PetscCall(MatDiagonalScale(C, d, NULL));
92 }
93 if (mats_view) {
94 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Diagonal Scaled C:\n"));
95 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
96 }
97
98 /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */
99 PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
100 PetscCall(MatSetType(B, MATSCALAPACK));
101 PetscCall(MatSetSizes(B, mloc, nloc, PETSC_DECIDE, PETSC_DECIDE));
102 PetscCall(MatScaLAPACKSetBlockSizes(B, mb, nb));
103 PetscCall(MatSetFromOptions(B));
104 PetscCall(MatSetUp(B));
105 /* PetscCall(MatCreateScaLAPACK(PETSC_COMM_WORLD,mb,nb,M,N,0,0,&B)); */
106 PetscCall(MatGetOwnershipIS(B, &isrows, &iscols));
107 PetscCall(ISGetLocalSize(isrows, &nrows));
108 PetscCall(ISGetIndices(isrows, &rows));
109 PetscCall(ISGetLocalSize(iscols, &ncols));
110 PetscCall(ISGetIndices(iscols, &cols));
111 PetscCall(PetscMalloc1(nrows * ncols, &v));
112 for (i = 0; i < nrows; i++) {
113 for (j = 0; j < ncols; j++) v[i * ncols + j] = (PetscReal)(1000 * rows[i] + cols[j]);
114 }
115 PetscCall(MatSetValues(B, nrows, rows, ncols, cols, v, INSERT_VALUES));
116 PetscCall(PetscFree(v));
117 PetscCall(ISRestoreIndices(isrows, &rows));
118 PetscCall(ISRestoreIndices(iscols, &cols));
119 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
120 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
121 if (mats_view) {
122 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B:\n"));
123 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
124 }
125 PetscCall(MatAXPY(B, 2.5, C, SAME_NONZERO_PATTERN));
126 PetscCall(MatAYPX(B, 3.75, C, SAME_NONZERO_PATTERN));
127 PetscCall(MatConvert(B, MATDENSE, MAT_INPLACE_MATRIX, &B));
128 if (mats_view) {
129 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B after MatAXPY and MatAYPX:\n"));
130 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
131 }
132 PetscCall(ISDestroy(&isrows));
133 PetscCall(ISDestroy(&iscols));
134 PetscCall(MatDestroy(&B));
135
136 /* Test MatMatTransposeMult(): B = C*C^T */
137 PetscCall(MatMatTransposeMult(C, C, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &B));
138 PetscCall(MatScale(C, 2.0));
139 PetscCall(MatMatTransposeMult(C, C, MAT_REUSE_MATRIX, PETSC_DETERMINE, &B));
140 PetscCall(MatMatTransposeMultEqual(C, C, B, 10, &flg));
141 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Check fails: B != C*C^T");
142
143 if (mats_view) {
144 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C MatMatTransposeMult C:\n"));
145 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
146 }
147 PetscCall(MatDestroy(&B));
148
149 /* Test MatTransposeMatMult(): B = C^T*C */
150 PetscCall(MatTransposeMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &B));
151 PetscCall(MatScale(C, 2.0));
152 PetscCall(MatTransposeMatMult(C, C, MAT_REUSE_MATRIX, PETSC_DETERMINE, &B));
153 PetscCall(MatTransposeMatMultEqual(C, C, B, 10, &flg));
154 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Check fails: B != C^T*C");
155
156 if (mats_view) {
157 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C MatTransposeMatMult C:\n"));
158 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
159 }
160
161 /* Test MatMult() */
162 PetscCall(MatComputeOperator(C, MATAIJ, &Caij));
163 PetscCall(MatMultEqual(C, Caij, 5, &flg));
164 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultEqual() fails");
165 PetscCall(MatMultTransposeEqual(C, Caij, 5, &flg));
166 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultTransposeEqual() fails");
167
168 /* Test MatMultAdd() and MatMultTransposeAddEqual() */
169 PetscCall(MatMultAddEqual(C, Caij, 5, &flg));
170 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultAddEqual() fails");
171 PetscCall(MatMultTransposeAddEqual(C, Caij, 5, &flg));
172 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultTransposeAddEqual() fails");
173
174 /* Test MatMatMult() */
175 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_matmatmult", &flg));
176 if (flg) {
177 Mat CC, CCaij;
178 PetscCall(MatMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &CC));
179 PetscCall(MatMatMult(Caij, Caij, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &CCaij));
180 PetscCall(MatMultEqual(CC, CCaij, 5, &flg));
181 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "CC != CCaij. MatMatMult() fails");
182 PetscCall(MatDestroy(&CCaij));
183 PetscCall(MatDestroy(&CC));
184 }
185
186 /* Test MatCreate() on subcomm */
187 color = rank % 2;
188 PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, color, 0, &subcomm));
189 if (color == 0) {
190 PetscCall(MatCreate(subcomm, &Asub));
191 PetscCall(MatSetType(Asub, MATSCALAPACK));
192 mloc = PETSC_DECIDE;
193 PetscCall(PetscSplitOwnershipEqual(subcomm, &mloc, &M));
194 nloc = PETSC_DECIDE;
195 PetscCall(PetscSplitOwnershipEqual(subcomm, &nloc, &N));
196 PetscCall(MatSetSizes(Asub, mloc, nloc, M, N));
197 PetscCall(MatScaLAPACKSetBlockSizes(Asub, mb, nb));
198 PetscCall(MatSetFromOptions(Asub));
199 PetscCall(MatSetUp(Asub));
200 PetscCall(MatDestroy(&Asub));
201 }
202
203 PetscCall(MatDestroy(&Cdense));
204 PetscCall(MatDestroy(&Caij));
205 PetscCall(MatDestroy(&B));
206 PetscCall(MatDestroy(&C));
207 PetscCall(MatDestroy(&Ct));
208 PetscCall(VecDestroy(&d));
209 PetscCallMPI(MPI_Comm_free(&subcomm));
210 PetscCall(PetscFinalize());
211 return 0;
212 }
213
214 /*TEST
215
216 build:
217 requires: scalapack double
218
219 testset:
220 output_file: output/empty.out
221
222 test:
223 nsize: 2
224 args: -mb 5 -nb 5 -M 12 -N 10
225
226 test:
227 suffix: 2
228 nsize: 6
229 args: -mb 8 -nb 6 -M 20 -N 50
230
231 test:
232 suffix: 3
233 nsize: 3
234 args: -mb 2 -nb 2 -M 20 -N 20 -test_matmatmult
235
236 TEST*/
237