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