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