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