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*/