1 static char help[] = "Tests LU, Cholesky, and QR factorization and MatMatSolve() for a sequential dense matrix. \n\
2 For MATSEQDENSE matrix, the factorization is just a thin wrapper to LAPACK. \n\
3 For MATSEQDENSECUDA, it uses cusolverDn routines \n\n";
4
5 #include <petscmat.h>
6
createMatsAndVecs(PetscInt m,PetscInt n,PetscInt nrhs,PetscBool full,Mat * _mat,Mat * _RHS,Mat * _SOLU,Vec * _x,Vec * _y,Vec * _b)7 static PetscErrorCode createMatsAndVecs(PetscInt m, PetscInt n, PetscInt nrhs, PetscBool full, Mat *_mat, Mat *_RHS, Mat *_SOLU, Vec *_x, Vec *_y, Vec *_b)
8 {
9 PetscRandom rand;
10 Mat mat, RHS, SOLU;
11 PetscInt rstart, rend;
12 PetscInt cstart, cend;
13 PetscScalar value = 1.0;
14 Vec x, y, b;
15
16 PetscFunctionBegin;
17 /* create multiple vectors RHS and SOLU */
18 PetscCall(MatCreate(PETSC_COMM_WORLD, &RHS));
19 PetscCall(MatSetSizes(RHS, PETSC_DECIDE, PETSC_DECIDE, m, nrhs));
20 PetscCall(MatSetType(RHS, MATDENSE));
21 PetscCall(MatSetOptionsPrefix(RHS, "rhs_"));
22 PetscCall(MatSetFromOptions(RHS));
23 PetscCall(MatSeqDenseSetPreallocation(RHS, NULL));
24
25 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
26 PetscCall(PetscRandomSetFromOptions(rand));
27 PetscCall(MatSetRandom(RHS, rand));
28
29 if (m == n) {
30 PetscCall(MatDuplicate(RHS, MAT_DO_NOT_COPY_VALUES, &SOLU));
31 } else {
32 PetscCall(MatCreate(PETSC_COMM_WORLD, &SOLU));
33 PetscCall(MatSetSizes(SOLU, PETSC_DECIDE, PETSC_DECIDE, n, nrhs));
34 PetscCall(MatSetType(SOLU, MATDENSE));
35 PetscCall(MatSeqDenseSetPreallocation(SOLU, NULL));
36 }
37 PetscCall(MatSetRandom(SOLU, rand));
38
39 /* create matrix */
40 PetscCall(MatCreate(PETSC_COMM_WORLD, &mat));
41 PetscCall(MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, m, n));
42 PetscCall(MatSetType(mat, MATDENSE));
43 PetscCall(MatSetFromOptions(mat));
44 PetscCall(MatSetUp(mat));
45 PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
46 PetscCall(MatGetOwnershipRangeColumn(mat, &cstart, &cend));
47 if (!full) {
48 for (PetscInt i = rstart; i < rend; i++) {
49 if (m == n) {
50 value = (PetscReal)i + 1;
51 PetscCall(MatSetValues(mat, 1, &i, 1, &i, &value, INSERT_VALUES));
52 } else {
53 for (PetscInt j = cstart; j < cend; j++) {
54 value = ((PetscScalar)i + 1.) / (PetscSqr(i - j) + 1.);
55 PetscCall(MatSetValues(mat, 1, &i, 1, &j, &value, INSERT_VALUES));
56 }
57 }
58 }
59 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
60 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
61 } else {
62 PetscCall(MatSetRandom(mat, rand));
63 if (m == n) {
64 Mat T;
65
66 PetscCall(MatMatTransposeMult(mat, mat, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &T));
67 PetscCall(MatDestroy(&mat));
68 mat = T;
69 }
70 }
71
72 /* create single vectors */
73 PetscCall(MatCreateVecs(mat, &x, &b));
74 PetscCall(VecDuplicate(x, &y));
75 PetscCall(VecSet(x, value));
76 PetscCall(PetscRandomDestroy(&rand));
77 *_mat = mat;
78 *_RHS = RHS;
79 *_SOLU = SOLU;
80 *_x = x;
81 *_y = y;
82 *_b = b;
83 PetscFunctionReturn(PETSC_SUCCESS);
84 }
85
main(int argc,char ** argv)86 int main(int argc, char **argv)
87 {
88 Mat mat, F, RHS, SOLU;
89 MatInfo info;
90 PetscInt m = 15, n = 10, i, j, nrhs = 2;
91 Vec x, y, b, ytmp;
92 IS perm;
93 PetscReal norm, tol = PETSC_SMALL;
94 PetscMPIInt size;
95 char solver[64];
96 PetscBool inplace, full = PETSC_FALSE, ldl = PETSC_TRUE, qr = PETSC_TRUE;
97
98 PetscFunctionBeginUser;
99 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
100 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
101 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
102 PetscCall(PetscStrncpy(solver, MATSOLVERPETSC, sizeof(solver)));
103 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
104 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
105 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
106 PetscCall(PetscOptionsGetBool(NULL, NULL, "-ldl", &ldl, NULL));
107 PetscCall(PetscOptionsGetBool(NULL, NULL, "-qr", &qr, NULL));
108 PetscCall(PetscOptionsGetBool(NULL, NULL, "-full", &full, NULL));
109 PetscCall(PetscOptionsGetString(NULL, NULL, "-solver_type", solver, sizeof(solver), NULL));
110
111 PetscCall(createMatsAndVecs(n, n, nrhs, full, &mat, &RHS, &SOLU, &x, &y, &b));
112 PetscCall(VecDuplicate(y, &ytmp));
113
114 /* Only SeqDense* support in-place factorizations and NULL permutations */
115 PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQDENSE, &inplace));
116 PetscCall(MatGetLocalSize(mat, &i, NULL));
117 PetscCall(MatGetOwnershipRange(mat, &j, NULL));
118 PetscCall(ISCreateStride(PETSC_COMM_WORLD, i, j, 1, &perm));
119
120 PetscCall(MatGetInfo(mat, MAT_LOCAL, &info));
121 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "matrix nonzeros = %" PetscInt_FMT ", allocated nonzeros = %" PetscInt_FMT "\n", (PetscInt)info.nz_used, (PetscInt)info.nz_allocated));
122 PetscCall(MatMult(mat, x, b));
123
124 /* Cholesky factorization - perm and factinfo are ignored by LAPACK */
125 /* in-place Cholesky */
126 if (inplace) {
127 Mat RHS2;
128
129 if (!ldl) PetscCall(MatSetOption(mat, MAT_SPD, PETSC_TRUE));
130 PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &F));
131 PetscCall(MatCholeskyFactor(F, perm, 0));
132 PetscCall(MatSolve(F, b, y));
133 PetscCall(VecAXPY(y, -1.0, x));
134 PetscCall(VecNorm(y, NORM_2, &norm));
135 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for in-place Cholesky %g\n", (double)norm));
136
137 PetscCall(MatMatSolve(F, RHS, SOLU));
138 PetscCall(MatMatMult(mat, SOLU, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &RHS2));
139 PetscCall(MatAXPY(RHS, -1.0, RHS2, SAME_NONZERO_PATTERN));
140 PetscCall(MatNorm(RHS, NORM_FROBENIUS, &norm));
141 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: Norm of residual for in-place Cholesky (MatMatSolve) %g\n", (double)norm));
142 PetscCall(MatDestroy(&F));
143 PetscCall(MatDestroy(&RHS2));
144 }
145
146 /* out-of-place Cholesky */
147 if (!ldl) PetscCall(MatSetOption(mat, MAT_SPD, PETSC_TRUE));
148 PetscCall(MatGetFactor(mat, solver, MAT_FACTOR_CHOLESKY, &F));
149 PetscCall(MatCholeskyFactorSymbolic(F, mat, perm, 0));
150 PetscCall(MatCholeskyFactorNumeric(F, mat, 0));
151 PetscCall(MatSolve(F, b, y));
152 PetscCall(VecAXPY(y, -1.0, x));
153 PetscCall(VecNorm(y, NORM_2, &norm));
154 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for out-of-place Cholesky %g\n", (double)norm));
155 PetscCall(MatDestroy(&F));
156
157 /* LU factorization - perms and factinfo are ignored by LAPACK */
158 i = n - 1;
159 PetscCall(MatZeroRows(mat, 1, &i, -1.0, NULL, NULL));
160 PetscCall(MatMult(mat, x, b));
161
162 /* in-place LU */
163 if (inplace) {
164 Mat RHS2;
165
166 PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &F));
167 PetscCall(MatLUFactor(F, perm, perm, 0));
168 PetscCall(MatSolve(F, b, y));
169 PetscCall(VecAXPY(y, -1.0, x));
170 PetscCall(VecNorm(y, NORM_2, &norm));
171 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for in-place LU %g\n", (double)norm));
172 PetscCall(MatMatSolve(F, RHS, SOLU));
173 PetscCall(MatMatMult(mat, SOLU, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &RHS2));
174 PetscCall(MatAXPY(RHS, -1.0, RHS2, SAME_NONZERO_PATTERN));
175 PetscCall(MatNorm(RHS, NORM_FROBENIUS, &norm));
176 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: Norm of residual for in-place LU (MatMatSolve) %g\n", (double)norm));
177 PetscCall(MatDestroy(&F));
178 PetscCall(MatDestroy(&RHS2));
179 }
180
181 /* out-of-place LU */
182 PetscCall(MatGetFactor(mat, solver, MAT_FACTOR_LU, &F));
183 PetscCall(MatLUFactorSymbolic(F, mat, perm, perm, 0));
184 PetscCall(MatLUFactorNumeric(F, mat, 0));
185 PetscCall(MatSolve(F, b, y));
186 PetscCall(VecAXPY(y, -1.0, x));
187 PetscCall(VecNorm(y, NORM_2, &norm));
188 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for out-of-place LU %g\n", (double)norm));
189
190 /* free space */
191 PetscCall(ISDestroy(&perm));
192 PetscCall(MatDestroy(&F));
193 PetscCall(MatDestroy(&mat));
194 PetscCall(MatDestroy(&RHS));
195 PetscCall(MatDestroy(&SOLU));
196 PetscCall(VecDestroy(&x));
197 PetscCall(VecDestroy(&b));
198 PetscCall(VecDestroy(&y));
199 PetscCall(VecDestroy(&ytmp));
200
201 if (qr) {
202 /* setup rectangular */
203 PetscCall(createMatsAndVecs(m, n, nrhs, full, &mat, &RHS, &SOLU, &x, &y, &b));
204 PetscCall(VecDuplicate(y, &ytmp));
205
206 /* QR factorization - perms and factinfo are ignored by LAPACK */
207 PetscCall(MatMult(mat, x, b));
208
209 /* in-place QR */
210 if (inplace) {
211 Mat SOLU2;
212
213 PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &F));
214 PetscCall(MatQRFactor(F, NULL, 0));
215 PetscCall(MatSolve(F, b, y));
216 PetscCall(VecAXPY(y, -1.0, x));
217 PetscCall(VecNorm(y, NORM_2, &norm));
218 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for in-place QR %g\n", (double)norm));
219 PetscCall(MatMatMult(mat, SOLU, MAT_REUSE_MATRIX, PETSC_DETERMINE, &RHS));
220 PetscCall(MatDuplicate(SOLU, MAT_DO_NOT_COPY_VALUES, &SOLU2));
221 PetscCall(MatMatSolve(F, RHS, SOLU2));
222 PetscCall(MatAXPY(SOLU2, -1.0, SOLU, SAME_NONZERO_PATTERN));
223 PetscCall(MatNorm(SOLU2, NORM_FROBENIUS, &norm));
224 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: Norm of error for in-place QR (MatMatSolve) %g\n", (double)norm));
225 PetscCall(MatDestroy(&F));
226 PetscCall(MatDestroy(&SOLU2));
227 }
228
229 /* out-of-place QR */
230 PetscCall(MatGetFactor(mat, solver, MAT_FACTOR_QR, &F));
231 PetscCall(MatQRFactorSymbolic(F, mat, NULL, NULL));
232 PetscCall(MatQRFactorNumeric(F, mat, NULL));
233 PetscCall(MatSolve(F, b, y));
234 PetscCall(VecAXPY(y, -1.0, x));
235 PetscCall(VecNorm(y, NORM_2, &norm));
236 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for out-of-place QR %g\n", (double)norm));
237
238 if (m == n) {
239 /* out-of-place MatSolveTranspose */
240 PetscCall(MatMultTranspose(mat, x, b));
241 PetscCall(MatSolveTranspose(F, b, y));
242 PetscCall(VecAXPY(y, -1.0, x));
243 PetscCall(VecNorm(y, NORM_2, &norm));
244 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for out-of-place QR %g\n", (double)norm));
245 }
246
247 /* free space */
248 PetscCall(MatDestroy(&F));
249 PetscCall(MatDestroy(&mat));
250 PetscCall(MatDestroy(&RHS));
251 PetscCall(MatDestroy(&SOLU));
252 PetscCall(VecDestroy(&x));
253 PetscCall(VecDestroy(&b));
254 PetscCall(VecDestroy(&y));
255 PetscCall(VecDestroy(&ytmp));
256 }
257 PetscCall(PetscFinalize());
258 return 0;
259 }
260
261 /*TEST
262
263 test:
264
265 test:
266 requires: cuda
267 suffix: seqdensecuda
268 args: -mat_type seqdensecuda -rhs_mat_type seqdensecuda -ldl 0 -solver_type {{petsc cuda}}
269 output_file: output/ex1_1.out
270
271 test:
272 requires: cuda
273 suffix: seqdensecuda_2
274 args: -ldl 0 -solver_type cuda
275 output_file: output/ex1_1.out
276
277 test:
278 requires: cuda
279 suffix: seqdensecuda_seqaijcusparse
280 args: -mat_type seqaijcusparse -rhs_mat_type seqdensecuda -qr 0
281 output_file: output/ex1_2.out
282
283 test:
284 requires: cuda viennacl
285 suffix: seqdensecuda_seqaijviennacl
286 args: -mat_type seqaijviennacl -rhs_mat_type seqdensecuda -qr 0
287 output_file: output/ex1_2.out
288
289 test:
290 suffix: 4
291 args: -m 10 -n 10
292 output_file: output/ex1_1.out
293
294 TEST*/
295