1 2 static char help[] = "Tests LU, Cholesky factorization and MatMatSolve() for an Elemental dense matrix.\n\n"; 3 4 #include <petscmat.h> 5 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; 19 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)); 24 25 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand)); 26 PetscCall(PetscRandomSetFromOptions(rand)); 27 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)); 35 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 } 68 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 } 100 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)); 119 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)); 129 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 */ 136 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 } 152 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)); 166 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 } 173 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)); 178 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)); 187 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)); 196 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)); 202 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)); 219 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)); 224 225 /* Solve LUX = PB or LUX = B */ 226 PetscCall(MatSolveAdd(F, b, d, x)); 227 PetscCall(MatMatSolve(F, B, X)); 228 PetscCall(MatDestroy(&F)); 229 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)); 246 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)); 255 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 } 271 272 /*TEST 273 274 build: 275 requires: elemental 276 277 test: 278 nsize: 2 279 output_file: output/ex145.out 280 281 test: 282 suffix: 2 283 nsize: 6 284 output_file: output/ex145.out 285 286 TEST*/ 287