1 static char help[] = "Tests LU, Cholesky factorization and MatMatSolve() for an Elemental dense matrix.\n\n"; 2 3 #include <petscmat.h> 4 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.e-11; 12 PetscMPIInt size, rank; 13 PetscRandom rand; 14 const PetscInt *rows, *cols; 15 IS isrows, iscols; 16 PetscBool mats_view = PETSC_FALSE; 17 MatFactorInfo finfo; 18 19 PetscFunctionBeginUser; 20 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 21 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 22 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 23 24 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand)); 25 PetscCall(PetscRandomSetFromOptions(rand)); 26 27 /* Get local dimensions of matrices */ 28 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 29 n = m; 30 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 31 p = m / 2; 32 PetscCall(PetscOptionsGetInt(NULL, NULL, "-p", &p, NULL)); 33 PetscCall(PetscOptionsHasName(NULL, NULL, "-mats_view", &mats_view)); 34 35 /* Create matrix A */ 36 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Create Elemental matrix A\n")); 37 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 38 PetscCall(MatSetSizes(A, m, n, PETSC_DECIDE, PETSC_DECIDE)); 39 PetscCall(MatSetType(A, MATELEMENTAL)); 40 PetscCall(MatSetFromOptions(A)); 41 PetscCall(MatSetUp(A)); 42 /* Set local matrix entries */ 43 PetscCall(MatGetOwnershipIS(A, &isrows, &iscols)); 44 PetscCall(ISGetLocalSize(isrows, &nrows)); 45 PetscCall(ISGetIndices(isrows, &rows)); 46 PetscCall(ISGetLocalSize(iscols, &ncols)); 47 PetscCall(ISGetIndices(iscols, &cols)); 48 PetscCall(PetscMalloc1(nrows * ncols, &v)); 49 for (i = 0; i < nrows; i++) { 50 for (j = 0; j < ncols; j++) { 51 PetscCall(PetscRandomGetValue(rand, &rval)); 52 v[i * ncols + j] = rval; 53 } 54 } 55 PetscCall(MatSetValues(A, nrows, rows, ncols, cols, v, INSERT_VALUES)); 56 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 57 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 58 PetscCall(ISRestoreIndices(isrows, &rows)); 59 PetscCall(ISRestoreIndices(iscols, &cols)); 60 PetscCall(ISDestroy(&isrows)); 61 PetscCall(ISDestroy(&iscols)); 62 PetscCall(PetscFree(v)); 63 if (mats_view) { 64 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "A: nrows %" PetscInt_FMT ", m %" PetscInt_FMT "; ncols %" PetscInt_FMT ", n %" PetscInt_FMT "\n", nrows, m, ncols, n)); 65 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 66 } 67 68 /* Create rhs matrix B */ 69 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Create rhs matrix B\n")); 70 PetscCall(MatCreate(PETSC_COMM_WORLD, &B)); 71 PetscCall(MatSetSizes(B, m, p, PETSC_DECIDE, PETSC_DECIDE)); 72 PetscCall(MatSetType(B, MATELEMENTAL)); 73 PetscCall(MatSetFromOptions(B)); 74 PetscCall(MatSetUp(B)); 75 PetscCall(MatGetOwnershipIS(B, &isrows, &iscols)); 76 PetscCall(ISGetLocalSize(isrows, &nrows)); 77 PetscCall(ISGetIndices(isrows, &rows)); 78 PetscCall(ISGetLocalSize(iscols, &ncols)); 79 PetscCall(ISGetIndices(iscols, &cols)); 80 PetscCall(PetscMalloc1(nrows * ncols, &v)); 81 for (i = 0; i < nrows; i++) { 82 for (j = 0; j < ncols; j++) { 83 PetscCall(PetscRandomGetValue(rand, &rval)); 84 v[i * ncols + j] = rval; 85 } 86 } 87 PetscCall(MatSetValues(B, nrows, rows, ncols, cols, v, INSERT_VALUES)); 88 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 89 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 90 PetscCall(ISRestoreIndices(isrows, &rows)); 91 PetscCall(ISRestoreIndices(iscols, &cols)); 92 PetscCall(ISDestroy(&isrows)); 93 PetscCall(ISDestroy(&iscols)); 94 PetscCall(PetscFree(v)); 95 if (mats_view) { 96 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B: nrows %" PetscInt_FMT ", m %" PetscInt_FMT "; ncols %" PetscInt_FMT ", p %" PetscInt_FMT "\n", nrows, m, ncols, p)); 97 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD)); 98 } 99 100 /* Create rhs vector b and solution x (same size as b) */ 101 PetscCall(VecCreate(PETSC_COMM_WORLD, &b)); 102 PetscCall(VecSetSizes(b, m, PETSC_DECIDE)); 103 PetscCall(VecSetFromOptions(b)); 104 PetscCall(VecGetArray(b, &barray)); 105 for (j = 0; j < m; j++) { 106 PetscCall(PetscRandomGetValue(rand, &rval)); 107 barray[j] = rval; 108 } 109 PetscCall(VecRestoreArray(b, &barray)); 110 PetscCall(VecAssemblyBegin(b)); 111 PetscCall(VecAssemblyEnd(b)); 112 if (mats_view) { 113 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] b: m %" PetscInt_FMT "\n", rank, m)); 114 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 115 PetscCall(VecView(b, PETSC_VIEWER_STDOUT_WORLD)); 116 } 117 PetscCall(VecDuplicate(b, &x)); 118 119 /* Create matrix X - same size as B */ 120 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Create solution matrix X\n")); 121 PetscCall(MatCreate(PETSC_COMM_WORLD, &X)); 122 PetscCall(MatSetSizes(X, m, p, PETSC_DECIDE, PETSC_DECIDE)); 123 PetscCall(MatSetType(X, MATELEMENTAL)); 124 PetscCall(MatSetFromOptions(X)); 125 PetscCall(MatSetUp(X)); 126 PetscCall(MatAssemblyBegin(X, MAT_FINAL_ASSEMBLY)); 127 PetscCall(MatAssemblyEnd(X, MAT_FINAL_ASSEMBLY)); 128 129 /* Cholesky factorization */ 130 /*------------------------*/ 131 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Create Elemental matrix Aher\n")); 132 PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &Aher)); 133 PetscCall(MatAXPY(Aher, 1.0, A, SAME_NONZERO_PATTERN)); /* Aher = A + A^T */ 134 if (rank == 0) { /* add 100.0 to diagonals of Aher to make it spd */ 135 136 /* TODO: Replace this with a call to El::ShiftDiagonal( A, 100.), 137 or at least pre-allocate the right amount of space */ 138 PetscInt M, N; 139 PetscCall(MatGetSize(Aher, &M, &N)); 140 for (i = 0; i < M; i++) { 141 rval = 100.0; 142 PetscCall(MatSetValues(Aher, 1, &i, 1, &i, &rval, ADD_VALUES)); 143 } 144 } 145 PetscCall(MatAssemblyBegin(Aher, MAT_FINAL_ASSEMBLY)); 146 PetscCall(MatAssemblyEnd(Aher, MAT_FINAL_ASSEMBLY)); 147 if (mats_view) { 148 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Aher:\n")); 149 PetscCall(MatView(Aher, PETSC_VIEWER_STDOUT_WORLD)); 150 } 151 152 /* Cholesky factorization */ 153 /*------------------------*/ 154 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Test Cholesky Solver \n")); 155 /* In-place Cholesky */ 156 /* Create matrix factor G, then copy Aher to G */ 157 PetscCall(MatCreate(PETSC_COMM_WORLD, &G)); 158 PetscCall(MatSetSizes(G, m, n, PETSC_DECIDE, PETSC_DECIDE)); 159 PetscCall(MatSetType(G, MATELEMENTAL)); 160 PetscCall(MatSetFromOptions(G)); 161 PetscCall(MatSetUp(G)); 162 PetscCall(MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY)); 163 PetscCall(MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY)); 164 PetscCall(MatCopy(Aher, G, SAME_NONZERO_PATTERN)); 165 166 /* Only G = U^T * U is implemented for now */ 167 PetscCall(MatCholeskyFactor(G, 0, 0)); 168 if (mats_view) { 169 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Cholesky Factor G:\n")); 170 PetscCall(MatView(G, PETSC_VIEWER_STDOUT_WORLD)); 171 } 172 173 /* Solve U^T * U x = b and U^T * U X = B */ 174 PetscCall(MatSolve(G, b, x)); 175 PetscCall(MatMatSolve(G, B, X)); 176 PetscCall(MatDestroy(&G)); 177 178 /* Out-place Cholesky */ 179 PetscCall(MatGetFactor(Aher, MATSOLVERELEMENTAL, MAT_FACTOR_CHOLESKY, &G)); 180 PetscCall(MatCholeskyFactorSymbolic(G, Aher, 0, &finfo)); 181 PetscCall(MatCholeskyFactorNumeric(G, Aher, &finfo)); 182 if (mats_view) PetscCall(MatView(G, PETSC_VIEWER_STDOUT_WORLD)); 183 PetscCall(MatSolve(G, b, x)); 184 PetscCall(MatMatSolve(G, B, X)); 185 PetscCall(MatDestroy(&G)); 186 187 /* Check norm(Aher*x - b) */ 188 PetscCall(VecCreate(PETSC_COMM_WORLD, &c)); 189 PetscCall(VecSetSizes(c, m, PETSC_DECIDE)); 190 PetscCall(VecSetFromOptions(c)); 191 PetscCall(MatMult(Aher, x, c)); 192 PetscCall(VecAXPY(c, -1.0, b)); 193 PetscCall(VecNorm(c, NORM_1, &norm)); 194 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: |Aher*x - b| for Cholesky %g\n", (double)norm)); 195 196 /* Check norm(Aher*X - B) */ 197 PetscCall(MatMatMult(Aher, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C)); 198 PetscCall(MatAXPY(C, -1.0, B, SAME_NONZERO_PATTERN)); 199 PetscCall(MatNorm(C, NORM_1, &norm)); 200 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: |Aher*X - B| for Cholesky %g\n", (double)norm)); 201 202 /* LU factorization */ 203 /*------------------*/ 204 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Test LU Solver \n")); 205 /* In-place LU */ 206 /* Create matrix factor F, then copy A to F */ 207 PetscCall(MatCreate(PETSC_COMM_WORLD, &F)); 208 PetscCall(MatSetSizes(F, m, n, PETSC_DECIDE, PETSC_DECIDE)); 209 PetscCall(MatSetType(F, MATELEMENTAL)); 210 PetscCall(MatSetFromOptions(F)); 211 PetscCall(MatSetUp(F)); 212 PetscCall(MatAssemblyBegin(F, MAT_FINAL_ASSEMBLY)); 213 PetscCall(MatAssemblyEnd(F, MAT_FINAL_ASSEMBLY)); 214 PetscCall(MatCopy(A, F, SAME_NONZERO_PATTERN)); 215 /* Create vector d to test MatSolveAdd() */ 216 PetscCall(VecDuplicate(x, &d)); 217 PetscCall(VecCopy(x, d)); 218 219 /* PF=LU or F=LU factorization - perms is ignored by Elemental; 220 set finfo.dtcol !0 or 0 to enable/disable partial pivoting */ 221 finfo.dtcol = 0.1; 222 PetscCall(MatLUFactor(F, 0, 0, &finfo)); 223 224 /* Solve LUX = PB or LUX = B */ 225 PetscCall(MatSolveAdd(F, b, d, x)); 226 PetscCall(MatMatSolve(F, B, X)); 227 PetscCall(MatDestroy(&F)); 228 229 /* Check norm(A*X - B) */ 230 PetscCall(VecCreate(PETSC_COMM_WORLD, &e)); 231 PetscCall(VecSetSizes(e, m, PETSC_DECIDE)); 232 PetscCall(VecSetFromOptions(e)); 233 PetscCall(MatMult(A, x, c)); 234 PetscCall(MatMult(A, d, e)); 235 PetscCall(VecAXPY(c, -1.0, e)); 236 PetscCall(VecAXPY(c, -1.0, b)); 237 PetscCall(VecNorm(c, NORM_1, &norm)); 238 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: |A*x - b| for LU %g\n", (double)norm)); 239 /* Reuse product C; replace Aher with A */ 240 PetscCall(MatProductReplaceMats(A, NULL, NULL, C)); 241 PetscCall(MatMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C)); 242 PetscCall(MatAXPY(C, -1.0, B, SAME_NONZERO_PATTERN)); 243 PetscCall(MatNorm(C, NORM_1, &norm)); 244 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: |A*X - B| for LU %g\n", (double)norm)); 245 246 /* Out-place LU */ 247 PetscCall(MatGetFactor(A, MATSOLVERELEMENTAL, MAT_FACTOR_LU, &F)); 248 PetscCall(MatLUFactorSymbolic(F, A, 0, 0, &finfo)); 249 PetscCall(MatLUFactorNumeric(F, A, &finfo)); 250 if (mats_view) PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD)); 251 PetscCall(MatSolve(F, b, x)); 252 PetscCall(MatMatSolve(F, B, X)); 253 PetscCall(MatDestroy(&F)); 254 255 /* Free space */ 256 PetscCall(MatDestroy(&A)); 257 PetscCall(MatDestroy(&Aher)); 258 PetscCall(MatDestroy(&B)); 259 PetscCall(MatDestroy(&C)); 260 PetscCall(MatDestroy(&X)); 261 PetscCall(VecDestroy(&b)); 262 PetscCall(VecDestroy(&c)); 263 PetscCall(VecDestroy(&d)); 264 PetscCall(VecDestroy(&e)); 265 PetscCall(VecDestroy(&x)); 266 PetscCall(PetscRandomDestroy(&rand)); 267 PetscCall(PetscFinalize()); 268 return 0; 269 } 270 271 /*TEST 272 273 build: 274 requires: elemental 275 276 test: 277 nsize: 2 278 output_file: output/ex145.out 279 280 test: 281 suffix: 2 282 nsize: 6 283 output_file: output/ex145.out 284 285 TEST*/ 286