1 static char help[] = "Tests LU, Cholesky factorization and MatMatSolve() for an Elemental dense matrix.\n\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** argv)5 int main(int argc, char **argv)
6 {
7 Mat A, F, B, X, C, Aher, G;
8 Vec b, x, c, d, e;
9 PetscInt m = 5, n, p, i, j, nrows, ncols;
10 PetscScalar *v, *barray, rval;
11 PetscReal norm, tol = 1.e-11;
12 PetscMPIInt size, rank;
13 PetscRandom rand;
14 const PetscInt *rows, *cols;
15 IS isrows, iscols;
16 PetscBool mats_view = PETSC_FALSE;
17 MatFactorInfo finfo;
18
19 PetscFunctionBeginUser;
20 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
21 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
22 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
23
24 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
25 PetscCall(PetscRandomSetFromOptions(rand));
26
27 /* Get local dimensions of matrices */
28 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
29 n = m;
30 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
31 p = m / 2;
32 PetscCall(PetscOptionsGetInt(NULL, NULL, "-p", &p, NULL));
33 PetscCall(PetscOptionsHasName(NULL, NULL, "-mats_view", &mats_view));
34
35 /* Create matrix A */
36 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Create Elemental matrix A\n"));
37 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
38 PetscCall(MatSetSizes(A, m, n, PETSC_DECIDE, PETSC_DECIDE));
39 PetscCall(MatSetType(A, MATELEMENTAL));
40 PetscCall(MatSetFromOptions(A));
41 PetscCall(MatSetUp(A));
42 /* Set local matrix entries */
43 PetscCall(MatGetOwnershipIS(A, &isrows, &iscols));
44 PetscCall(ISGetLocalSize(isrows, &nrows));
45 PetscCall(ISGetIndices(isrows, &rows));
46 PetscCall(ISGetLocalSize(iscols, &ncols));
47 PetscCall(ISGetIndices(iscols, &cols));
48 PetscCall(PetscMalloc1(nrows * ncols, &v));
49 for (i = 0; i < nrows; i++) {
50 for (j = 0; j < ncols; j++) {
51 PetscCall(PetscRandomGetValue(rand, &rval));
52 v[i * ncols + j] = rval;
53 }
54 }
55 PetscCall(MatSetValues(A, nrows, rows, ncols, cols, v, INSERT_VALUES));
56 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
57 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
58 PetscCall(ISRestoreIndices(isrows, &rows));
59 PetscCall(ISRestoreIndices(iscols, &cols));
60 PetscCall(ISDestroy(&isrows));
61 PetscCall(ISDestroy(&iscols));
62 PetscCall(PetscFree(v));
63 if (mats_view) {
64 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "A: nrows %" PetscInt_FMT ", m %" PetscInt_FMT "; ncols %" PetscInt_FMT ", n %" PetscInt_FMT "\n", nrows, m, ncols, n));
65 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
66 }
67
68 /* Create rhs matrix B */
69 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Create rhs matrix B\n"));
70 PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
71 PetscCall(MatSetSizes(B, m, p, PETSC_DECIDE, PETSC_DECIDE));
72 PetscCall(MatSetType(B, MATELEMENTAL));
73 PetscCall(MatSetFromOptions(B));
74 PetscCall(MatSetUp(B));
75 PetscCall(MatGetOwnershipIS(B, &isrows, &iscols));
76 PetscCall(ISGetLocalSize(isrows, &nrows));
77 PetscCall(ISGetIndices(isrows, &rows));
78 PetscCall(ISGetLocalSize(iscols, &ncols));
79 PetscCall(ISGetIndices(iscols, &cols));
80 PetscCall(PetscMalloc1(nrows * ncols, &v));
81 for (i = 0; i < nrows; i++) {
82 for (j = 0; j < ncols; j++) {
83 PetscCall(PetscRandomGetValue(rand, &rval));
84 v[i * ncols + j] = rval;
85 }
86 }
87 PetscCall(MatSetValues(B, nrows, rows, ncols, cols, v, INSERT_VALUES));
88 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
89 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
90 PetscCall(ISRestoreIndices(isrows, &rows));
91 PetscCall(ISRestoreIndices(iscols, &cols));
92 PetscCall(ISDestroy(&isrows));
93 PetscCall(ISDestroy(&iscols));
94 PetscCall(PetscFree(v));
95 if (mats_view) {
96 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B: nrows %" PetscInt_FMT ", m %" PetscInt_FMT "; ncols %" PetscInt_FMT ", p %" PetscInt_FMT "\n", nrows, m, ncols, p));
97 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
98 }
99
100 /* Create rhs vector b and solution x (same size as b) */
101 PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
102 PetscCall(VecSetSizes(b, m, PETSC_DECIDE));
103 PetscCall(VecSetFromOptions(b));
104 PetscCall(VecGetArray(b, &barray));
105 for (j = 0; j < m; j++) {
106 PetscCall(PetscRandomGetValue(rand, &rval));
107 barray[j] = rval;
108 }
109 PetscCall(VecRestoreArray(b, &barray));
110 PetscCall(VecAssemblyBegin(b));
111 PetscCall(VecAssemblyEnd(b));
112 if (mats_view) {
113 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] b: m %" PetscInt_FMT "\n", rank, m));
114 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
115 PetscCall(VecView(b, PETSC_VIEWER_STDOUT_WORLD));
116 }
117 PetscCall(VecDuplicate(b, &x));
118
119 /* Create matrix X - same size as B */
120 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Create solution matrix X\n"));
121 PetscCall(MatCreate(PETSC_COMM_WORLD, &X));
122 PetscCall(MatSetSizes(X, m, p, PETSC_DECIDE, PETSC_DECIDE));
123 PetscCall(MatSetType(X, MATELEMENTAL));
124 PetscCall(MatSetFromOptions(X));
125 PetscCall(MatSetUp(X));
126 PetscCall(MatAssemblyBegin(X, MAT_FINAL_ASSEMBLY));
127 PetscCall(MatAssemblyEnd(X, MAT_FINAL_ASSEMBLY));
128
129 /* Cholesky factorization */
130 /*------------------------*/
131 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Create Elemental matrix Aher\n"));
132 PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &Aher));
133 PetscCall(MatAXPY(Aher, 1.0, A, SAME_NONZERO_PATTERN)); /* Aher = A + A^T */
134 if (rank == 0) { /* add 100.0 to diagonals of Aher to make it spd */
135
136 /* TODO: Replace this with a call to El::ShiftDiagonal( A, 100.),
137 or at least pre-allocate the right amount of space */
138 PetscInt M, N;
139 PetscCall(MatGetSize(Aher, &M, &N));
140 for (i = 0; i < M; i++) {
141 rval = 100.0;
142 PetscCall(MatSetValues(Aher, 1, &i, 1, &i, &rval, ADD_VALUES));
143 }
144 }
145 PetscCall(MatAssemblyBegin(Aher, MAT_FINAL_ASSEMBLY));
146 PetscCall(MatAssemblyEnd(Aher, MAT_FINAL_ASSEMBLY));
147 if (mats_view) {
148 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Aher:\n"));
149 PetscCall(MatView(Aher, PETSC_VIEWER_STDOUT_WORLD));
150 }
151
152 /* Cholesky factorization */
153 /*------------------------*/
154 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Test Cholesky Solver \n"));
155 /* In-place Cholesky */
156 /* Create matrix factor G, then copy Aher to G */
157 PetscCall(MatCreate(PETSC_COMM_WORLD, &G));
158 PetscCall(MatSetSizes(G, m, n, PETSC_DECIDE, PETSC_DECIDE));
159 PetscCall(MatSetType(G, MATELEMENTAL));
160 PetscCall(MatSetFromOptions(G));
161 PetscCall(MatSetUp(G));
162 PetscCall(MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY));
163 PetscCall(MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY));
164 PetscCall(MatCopy(Aher, G, SAME_NONZERO_PATTERN));
165
166 /* Only G = U^T * U is implemented for now */
167 PetscCall(MatCholeskyFactor(G, 0, 0));
168 if (mats_view) {
169 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Cholesky Factor G:\n"));
170 PetscCall(MatView(G, PETSC_VIEWER_STDOUT_WORLD));
171 }
172
173 /* Solve U^T * U x = b and U^T * U X = B */
174 PetscCall(MatSolve(G, b, x));
175 PetscCall(MatMatSolve(G, B, X));
176 PetscCall(MatDestroy(&G));
177
178 /* Out-place Cholesky */
179 PetscCall(MatGetFactor(Aher, MATSOLVERELEMENTAL, MAT_FACTOR_CHOLESKY, &G));
180 PetscCall(MatCholeskyFactorSymbolic(G, Aher, 0, &finfo));
181 PetscCall(MatCholeskyFactorNumeric(G, Aher, &finfo));
182 if (mats_view) PetscCall(MatView(G, PETSC_VIEWER_STDOUT_WORLD));
183 PetscCall(MatSolve(G, b, x));
184 PetscCall(MatMatSolve(G, B, X));
185 PetscCall(MatDestroy(&G));
186
187 /* Check norm(Aher*x - b) */
188 PetscCall(VecCreate(PETSC_COMM_WORLD, &c));
189 PetscCall(VecSetSizes(c, m, PETSC_DECIDE));
190 PetscCall(VecSetFromOptions(c));
191 PetscCall(MatMult(Aher, x, c));
192 PetscCall(VecAXPY(c, -1.0, b));
193 PetscCall(VecNorm(c, NORM_1, &norm));
194 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: |Aher*x - b| for Cholesky %g\n", (double)norm));
195
196 /* Check norm(Aher*X - B) */
197 PetscCall(MatMatMult(Aher, X, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
198 PetscCall(MatAXPY(C, -1.0, B, SAME_NONZERO_PATTERN));
199 PetscCall(MatNorm(C, NORM_1, &norm));
200 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: |Aher*X - B| for Cholesky %g\n", (double)norm));
201
202 /* LU factorization */
203 /*------------------*/
204 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Test LU Solver \n"));
205 /* In-place LU */
206 /* Create matrix factor F, then copy A to F */
207 PetscCall(MatCreate(PETSC_COMM_WORLD, &F));
208 PetscCall(MatSetSizes(F, m, n, PETSC_DECIDE, PETSC_DECIDE));
209 PetscCall(MatSetType(F, MATELEMENTAL));
210 PetscCall(MatSetFromOptions(F));
211 PetscCall(MatSetUp(F));
212 PetscCall(MatAssemblyBegin(F, MAT_FINAL_ASSEMBLY));
213 PetscCall(MatAssemblyEnd(F, MAT_FINAL_ASSEMBLY));
214 PetscCall(MatCopy(A, F, SAME_NONZERO_PATTERN));
215 /* Create vector d to test MatSolveAdd() */
216 PetscCall(VecDuplicate(x, &d));
217 PetscCall(VecCopy(x, d));
218
219 /* PF=LU or F=LU factorization - perms is ignored by Elemental;
220 set finfo.dtcol !0 or 0 to enable/disable partial pivoting */
221 finfo.dtcol = 0.1;
222 PetscCall(MatLUFactor(F, 0, 0, &finfo));
223
224 /* Solve LUX = PB or LUX = B */
225 PetscCall(MatSolveAdd(F, b, d, x));
226 PetscCall(MatMatSolve(F, B, X));
227 PetscCall(MatDestroy(&F));
228
229 /* Check norm(A*X - B) */
230 PetscCall(VecCreate(PETSC_COMM_WORLD, &e));
231 PetscCall(VecSetSizes(e, m, PETSC_DECIDE));
232 PetscCall(VecSetFromOptions(e));
233 PetscCall(MatMult(A, x, c));
234 PetscCall(MatMult(A, d, e));
235 PetscCall(VecAXPY(c, -1.0, e));
236 PetscCall(VecAXPY(c, -1.0, b));
237 PetscCall(VecNorm(c, NORM_1, &norm));
238 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: |A*x - b| for LU %g\n", (double)norm));
239 /* Reuse product C; replace Aher with A */
240 PetscCall(MatProductReplaceMats(A, NULL, NULL, C));
241 PetscCall(MatMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DETERMINE, &C));
242 PetscCall(MatAXPY(C, -1.0, B, SAME_NONZERO_PATTERN));
243 PetscCall(MatNorm(C, NORM_1, &norm));
244 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: |A*X - B| for LU %g\n", (double)norm));
245
246 /* Out-place LU */
247 PetscCall(MatGetFactor(A, MATSOLVERELEMENTAL, MAT_FACTOR_LU, &F));
248 PetscCall(MatLUFactorSymbolic(F, A, 0, 0, &finfo));
249 PetscCall(MatLUFactorNumeric(F, A, &finfo));
250 if (mats_view) PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD));
251 PetscCall(MatSolve(F, b, x));
252 PetscCall(MatMatSolve(F, B, X));
253 PetscCall(MatDestroy(&F));
254
255 /* Free space */
256 PetscCall(MatDestroy(&A));
257 PetscCall(MatDestroy(&Aher));
258 PetscCall(MatDestroy(&B));
259 PetscCall(MatDestroy(&C));
260 PetscCall(MatDestroy(&X));
261 PetscCall(VecDestroy(&b));
262 PetscCall(VecDestroy(&c));
263 PetscCall(VecDestroy(&d));
264 PetscCall(VecDestroy(&e));
265 PetscCall(VecDestroy(&x));
266 PetscCall(PetscRandomDestroy(&rand));
267 PetscCall(PetscFinalize());
268 return 0;
269 }
270
271 /*TEST
272
273 build:
274 requires: elemental
275
276 test:
277 nsize: 2
278 output_file: output/ex145.out
279
280 test:
281 suffix: 2
282 nsize: 6
283 output_file: output/ex145.out
284
285 TEST*/
286