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;
21: PetscInitialize(&argc, &argv, (char *)0, help);
22: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
23: MPI_Comm_size(PETSC_COMM_WORLD, &size);
25: PetscRandomCreate(PETSC_COMM_WORLD, &rand);
26: PetscRandomSetFromOptions(rand);
28: /* Get local dimensions of matrices */
29: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
30: n = m;
31: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
32: p = m / 2;
33: PetscOptionsGetInt(NULL, NULL, "-p", &p, NULL);
34: PetscOptionsHasName(NULL, NULL, "-mats_view", &mats_view);
36: /* Create matrix A */
37: PetscPrintf(PETSC_COMM_WORLD, " Create Elemental matrix A\n");
38: MatCreate(PETSC_COMM_WORLD, &A);
39: MatSetSizes(A, m, n, PETSC_DECIDE, PETSC_DECIDE);
40: MatSetType(A, MATELEMENTAL);
41: MatSetFromOptions(A);
42: MatSetUp(A);
43: /* Set local matrix entries */
44: MatGetOwnershipIS(A, &isrows, &iscols);
45: ISGetLocalSize(isrows, &nrows);
46: ISGetIndices(isrows, &rows);
47: ISGetLocalSize(iscols, &ncols);
48: ISGetIndices(iscols, &cols);
49: PetscMalloc1(nrows * ncols, &v);
50: for (i = 0; i < nrows; i++) {
51: for (j = 0; j < ncols; j++) {
52: PetscRandomGetValue(rand, &rval);
53: v[i * ncols + j] = rval;
54: }
55: }
56: MatSetValues(A, nrows, rows, ncols, cols, v, INSERT_VALUES);
57: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
58: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
59: ISRestoreIndices(isrows, &rows);
60: ISRestoreIndices(iscols, &cols);
61: ISDestroy(&isrows);
62: ISDestroy(&iscols);
63: PetscFree(v);
64: if (mats_view) {
65: PetscPrintf(PETSC_COMM_WORLD, "A: nrows %" PetscInt_FMT ", m %" PetscInt_FMT "; ncols %" PetscInt_FMT ", n %" PetscInt_FMT "\n", nrows, m, ncols, n);
66: MatView(A, PETSC_VIEWER_STDOUT_WORLD);
67: }
69: /* Create rhs matrix B */
70: PetscPrintf(PETSC_COMM_WORLD, " Create rhs matrix B\n");
71: MatCreate(PETSC_COMM_WORLD, &B);
72: MatSetSizes(B, m, p, PETSC_DECIDE, PETSC_DECIDE);
73: MatSetType(B, MATELEMENTAL);
74: MatSetFromOptions(B);
75: MatSetUp(B);
76: MatGetOwnershipIS(B, &isrows, &iscols);
77: ISGetLocalSize(isrows, &nrows);
78: ISGetIndices(isrows, &rows);
79: ISGetLocalSize(iscols, &ncols);
80: ISGetIndices(iscols, &cols);
81: PetscMalloc1(nrows * ncols, &v);
82: for (i = 0; i < nrows; i++) {
83: for (j = 0; j < ncols; j++) {
84: PetscRandomGetValue(rand, &rval);
85: v[i * ncols + j] = rval;
86: }
87: }
88: MatSetValues(B, nrows, rows, ncols, cols, v, INSERT_VALUES);
89: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
90: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
91: ISRestoreIndices(isrows, &rows);
92: ISRestoreIndices(iscols, &cols);
93: ISDestroy(&isrows);
94: ISDestroy(&iscols);
95: PetscFree(v);
96: if (mats_view) {
97: PetscPrintf(PETSC_COMM_WORLD, "B: nrows %" PetscInt_FMT ", m %" PetscInt_FMT "; ncols %" PetscInt_FMT ", p %" PetscInt_FMT "\n", nrows, m, ncols, p);
98: MatView(B, PETSC_VIEWER_STDOUT_WORLD);
99: }
101: /* Create rhs vector b and solution x (same size as b) */
102: VecCreate(PETSC_COMM_WORLD, &b);
103: VecSetSizes(b, m, PETSC_DECIDE);
104: VecSetFromOptions(b);
105: VecGetArray(b, &barray);
106: for (j = 0; j < m; j++) {
107: PetscRandomGetValue(rand, &rval);
108: barray[j] = rval;
109: }
110: VecRestoreArray(b, &barray);
111: VecAssemblyBegin(b);
112: VecAssemblyEnd(b);
113: if (mats_view) {
114: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] b: m %" PetscInt_FMT "\n", rank, m);
115: PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
116: VecView(b, PETSC_VIEWER_STDOUT_WORLD);
117: }
118: VecDuplicate(b, &x);
120: /* Create matrix X - same size as B */
121: PetscPrintf(PETSC_COMM_WORLD, " Create solution matrix X\n");
122: MatCreate(PETSC_COMM_WORLD, &X);
123: MatSetSizes(X, m, p, PETSC_DECIDE, PETSC_DECIDE);
124: MatSetType(X, MATELEMENTAL);
125: MatSetFromOptions(X);
126: MatSetUp(X);
127: MatAssemblyBegin(X, MAT_FINAL_ASSEMBLY);
128: MatAssemblyEnd(X, MAT_FINAL_ASSEMBLY);
130: /* Cholesky factorization */
131: /*------------------------*/
132: PetscPrintf(PETSC_COMM_WORLD, " Create Elemental matrix Aher\n");
133: MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &Aher);
134: 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: MatGetSize(Aher, &M, &N);
141: for (i = 0; i < M; i++) {
142: rval = 100.0;
143: MatSetValues(Aher, 1, &i, 1, &i, &rval, ADD_VALUES);
144: }
145: }
146: MatAssemblyBegin(Aher, MAT_FINAL_ASSEMBLY);
147: MatAssemblyEnd(Aher, MAT_FINAL_ASSEMBLY);
148: if (mats_view) {
149: PetscPrintf(PETSC_COMM_WORLD, "Aher:\n");
150: MatView(Aher, PETSC_VIEWER_STDOUT_WORLD);
151: }
153: /* Cholesky factorization */
154: /*------------------------*/
155: PetscPrintf(PETSC_COMM_WORLD, " Test Cholesky Solver \n");
156: /* In-place Cholesky */
157: /* Create matrix factor G, then copy Aher to G */
158: MatCreate(PETSC_COMM_WORLD, &G);
159: MatSetSizes(G, m, n, PETSC_DECIDE, PETSC_DECIDE);
160: MatSetType(G, MATELEMENTAL);
161: MatSetFromOptions(G);
162: MatSetUp(G);
163: MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY);
164: MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY);
165: MatCopy(Aher, G, SAME_NONZERO_PATTERN);
167: /* Only G = U^T * U is implemented for now */
168: MatCholeskyFactor(G, 0, 0);
169: if (mats_view) {
170: PetscPrintf(PETSC_COMM_WORLD, "Cholesky Factor G:\n");
171: MatView(G, PETSC_VIEWER_STDOUT_WORLD);
172: }
174: /* Solve U^T * U x = b and U^T * U X = B */
175: MatSolve(G, b, x);
176: MatMatSolve(G, B, X);
177: MatDestroy(&G);
179: /* Out-place Cholesky */
180: MatGetFactor(Aher, MATSOLVERELEMENTAL, MAT_FACTOR_CHOLESKY, &G);
181: MatCholeskyFactorSymbolic(G, Aher, 0, &finfo);
182: MatCholeskyFactorNumeric(G, Aher, &finfo);
183: if (mats_view) MatView(G, PETSC_VIEWER_STDOUT_WORLD);
184: MatSolve(G, b, x);
185: MatMatSolve(G, B, X);
186: MatDestroy(&G);
188: /* Check norm(Aher*x - b) */
189: VecCreate(PETSC_COMM_WORLD, &c);
190: VecSetSizes(c, m, PETSC_DECIDE);
191: VecSetFromOptions(c);
192: MatMult(Aher, x, c);
193: VecAXPY(c, -1.0, b);
194: VecNorm(c, NORM_1, &norm);
195: if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "Warning: |Aher*x - b| for Cholesky %g\n", (double)norm);
197: /* Check norm(Aher*X - B) */
198: MatMatMult(Aher, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C);
199: MatAXPY(C, -1.0, B, SAME_NONZERO_PATTERN);
200: MatNorm(C, NORM_1, &norm);
201: if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "Warning: |Aher*X - B| for Cholesky %g\n", (double)norm);
203: /* LU factorization */
204: /*------------------*/
205: PetscPrintf(PETSC_COMM_WORLD, " Test LU Solver \n");
206: /* In-place LU */
207: /* Create matrix factor F, then copy A to F */
208: MatCreate(PETSC_COMM_WORLD, &F);
209: MatSetSizes(F, m, n, PETSC_DECIDE, PETSC_DECIDE);
210: MatSetType(F, MATELEMENTAL);
211: MatSetFromOptions(F);
212: MatSetUp(F);
213: MatAssemblyBegin(F, MAT_FINAL_ASSEMBLY);
214: MatAssemblyEnd(F, MAT_FINAL_ASSEMBLY);
215: MatCopy(A, F, SAME_NONZERO_PATTERN);
216: /* Create vector d to test MatSolveAdd() */
217: VecDuplicate(x, &d);
218: 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: MatLUFactor(F, 0, 0, &finfo);
225: /* Solve LUX = PB or LUX = B */
226: MatSolveAdd(F, b, d, x);
227: MatMatSolve(F, B, X);
228: MatDestroy(&F);
230: /* Check norm(A*X - B) */
231: VecCreate(PETSC_COMM_WORLD, &e);
232: VecSetSizes(e, m, PETSC_DECIDE);
233: VecSetFromOptions(e);
234: MatMult(A, x, c);
235: MatMult(A, d, e);
236: VecAXPY(c, -1.0, e);
237: VecAXPY(c, -1.0, b);
238: VecNorm(c, NORM_1, &norm);
239: if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "Warning: |A*x - b| for LU %g\n", (double)norm);
240: /* Reuse product C; replace Aher with A */
241: MatProductReplaceMats(A, NULL, NULL, C);
242: MatMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C);
243: MatAXPY(C, -1.0, B, SAME_NONZERO_PATTERN);
244: MatNorm(C, NORM_1, &norm);
245: if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "Warning: |A*X - B| for LU %g\n", (double)norm);
247: /* Out-place LU */
248: MatGetFactor(A, MATSOLVERELEMENTAL, MAT_FACTOR_LU, &F);
249: MatLUFactorSymbolic(F, A, 0, 0, &finfo);
250: MatLUFactorNumeric(F, A, &finfo);
251: if (mats_view) MatView(F, PETSC_VIEWER_STDOUT_WORLD);
252: MatSolve(F, b, x);
253: MatMatSolve(F, B, X);
254: MatDestroy(&F);
256: /* Free space */
257: MatDestroy(&A);
258: MatDestroy(&Aher);
259: MatDestroy(&B);
260: MatDestroy(&C);
261: MatDestroy(&X);
262: VecDestroy(&b);
263: VecDestroy(&c);
264: VecDestroy(&d);
265: VecDestroy(&e);
266: VecDestroy(&x);
267: PetscRandomDestroy(&rand);
268: 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*/