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