1 static char help[] = "Tests LU, Cholesky factorization and MatMatSolve() for a ScaLAPACK 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.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; 17 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)); 22 23 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand)); 24 PetscCall(PetscRandomSetFromOptions(rand)); 25 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)); 33 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 } 66 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 } 98 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)); 117 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)); 121 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 } 132 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)); 139 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 } 146 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)); 151 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)); 160 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)); 169 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)); 175 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)); 185 186 /* PF=LU factorization */ 187 PetscCall(MatLUFactor(F, 0, 0, NULL)); 188 189 /* Solve LUX = PB */ 190 PetscCall(MatSolveAdd(F, b, d, x)); 191 PetscCall(MatMatSolve(F, B, X)); 192 PetscCall(MatDestroy(&F)); 193 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)); 210 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)); 219 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 } 235 236 /*TEST 237 238 build: 239 requires: scalapack 240 241 test: 242 nsize: 2 243 output_file: output/ex245.out 244 245 test: 246 suffix: 2 247 nsize: 6 248 output_file: output/ex245.out 249 250 TEST*/ 251