Actual source code: ex245.c
1: static char help[] = "Tests LU, Cholesky factorization and MatMatSolve() for a ScaLAPACK dense matrix.\n\n";
3: #include <petscmat.h>
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;
18: PetscFunctionBeginUser;
19: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
20: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
21: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
23: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
24: PetscCall(PetscRandomSetFromOptions(rand));
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));
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: }
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: }
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));
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));
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: }
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));
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: }
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));
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));
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));
170: /* Check norm(Aher*X - B) */
171: PetscCall(MatMatMult(Aher, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &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));
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));
186: /* PF=LU factorization */
187: PetscCall(MatLUFactor(F, 0, 0, NULL));
189: /* Solve LUX = PB */
190: PetscCall(MatSolveAdd(F, b, d, x));
191: PetscCall(MatMatSolve(F, B, X));
192: PetscCall(MatDestroy(&F));
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_DEFAULT, &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));
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));
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: }
236: /*TEST
238: build:
239: requires: scalapack
241: test:
242: nsize: 2
243: output_file: output/ex245.out
245: test:
246: suffix: 2
247: nsize: 6
248: output_file: output/ex245.out
250: TEST*/