1 2 static char help[] = "Tests LU, Cholesky factorization and MatMatSolve() for a ScaLAPACK 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.e5 * PETSC_MACHINE_EPSILON; 13 PetscMPIInt size, rank; 14 PetscRandom rand; 15 const PetscInt *rows, *cols; 16 IS isrows, iscols; 17 PetscBool mats_view = PETSC_FALSE; 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 ScaLAPACK matrix A\n")); 37 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 38 PetscCall(MatSetSizes(A, m, n, PETSC_DECIDE, PETSC_DECIDE)); 39 PetscCall(MatSetType(A, MATSCALAPACK)); 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, MATSCALAPACK)); 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(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &X)); 122 123 /* Cholesky factorization */ 124 /*------------------------*/ 125 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Create ScaLAPACK matrix Aher\n")); 126 PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &Aher)); 127 PetscCall(MatAXPY(Aher, 1.0, A, SAME_NONZERO_PATTERN)); /* Aher = A + A^T */ 128 PetscCall(MatShift(Aher, 100.0)); /* add 100.0 to diagonals of Aher to make it spd */ 129 if (mats_view) { 130 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Aher:\n")); 131 PetscCall(MatView(Aher, PETSC_VIEWER_STDOUT_WORLD)); 132 } 133 134 /* Cholesky factorization */ 135 /*------------------------*/ 136 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Test Cholesky Solver \n")); 137 /* In-place Cholesky */ 138 /* Create matrix factor G, with a copy of Aher */ 139 PetscCall(MatDuplicate(Aher, MAT_COPY_VALUES, &G)); 140 141 /* G = L * L^T */ 142 PetscCall(MatCholeskyFactor(G, 0, 0)); 143 if (mats_view) { 144 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Cholesky Factor G:\n")); 145 PetscCall(MatView(G, PETSC_VIEWER_STDOUT_WORLD)); 146 } 147 148 /* Solve L * L^T x = b and L * L^T * X = B */ 149 PetscCall(MatSolve(G, b, x)); 150 PetscCall(MatMatSolve(G, B, X)); 151 PetscCall(MatDestroy(&G)); 152 153 /* Out-place Cholesky */ 154 PetscCall(MatGetFactor(Aher, MATSOLVERSCALAPACK, MAT_FACTOR_CHOLESKY, &G)); 155 PetscCall(MatCholeskyFactorSymbolic(G, Aher, 0, NULL)); 156 PetscCall(MatCholeskyFactorNumeric(G, Aher, NULL)); 157 if (mats_view) PetscCall(MatView(G, PETSC_VIEWER_STDOUT_WORLD)); 158 PetscCall(MatSolve(G, b, x)); 159 PetscCall(MatMatSolve(G, B, X)); 160 PetscCall(MatDestroy(&G)); 161 162 /* Check norm(Aher*x - b) */ 163 PetscCall(VecCreate(PETSC_COMM_WORLD, &c)); 164 PetscCall(VecSetSizes(c, m, PETSC_DECIDE)); 165 PetscCall(VecSetFromOptions(c)); 166 PetscCall(MatMult(Aher, x, c)); 167 PetscCall(VecAXPY(c, -1.0, b)); 168 PetscCall(VecNorm(c, NORM_1, &norm)); 169 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: ||Aher*x - b||=%g for Cholesky\n", (double)norm)); 170 171 /* Check norm(Aher*X - B) */ 172 PetscCall(MatMatMult(Aher, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C)); 173 PetscCall(MatAXPY(C, -1.0, B, SAME_NONZERO_PATTERN)); 174 PetscCall(MatNorm(C, NORM_1, &norm)); 175 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: ||Aher*X - B||=%g for Cholesky\n", (double)norm)); 176 177 /* LU factorization */ 178 /*------------------*/ 179 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Test LU Solver \n")); 180 /* In-place LU */ 181 /* Create matrix factor F, with a copy of A */ 182 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &F)); 183 /* Create vector d to test MatSolveAdd() */ 184 PetscCall(VecDuplicate(x, &d)); 185 PetscCall(VecCopy(x, d)); 186 187 /* PF=LU factorization */ 188 PetscCall(MatLUFactor(F, 0, 0, NULL)); 189 190 /* Solve LUX = PB */ 191 PetscCall(MatSolveAdd(F, b, d, x)); 192 PetscCall(MatMatSolve(F, B, X)); 193 PetscCall(MatDestroy(&F)); 194 195 /* Check norm(A*X - B) */ 196 PetscCall(VecCreate(PETSC_COMM_WORLD, &e)); 197 PetscCall(VecSetSizes(e, m, PETSC_DECIDE)); 198 PetscCall(VecSetFromOptions(e)); 199 PetscCall(MatMult(A, x, c)); 200 PetscCall(MatMult(A, d, e)); 201 PetscCall(VecAXPY(c, -1.0, e)); 202 PetscCall(VecAXPY(c, -1.0, b)); 203 PetscCall(VecNorm(c, NORM_1, &norm)); 204 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: ||A*x - b||=%g for LU\n", (double)norm)); 205 /* Reuse product C; replace Aher with A */ 206 PetscCall(MatProductReplaceMats(A, NULL, NULL, C)); 207 PetscCall(MatMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C)); 208 PetscCall(MatAXPY(C, -1.0, B, SAME_NONZERO_PATTERN)); 209 PetscCall(MatNorm(C, NORM_1, &norm)); 210 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: ||A*X - B||=%g for LU\n", (double)norm)); 211 212 /* Out-place LU */ 213 PetscCall(MatGetFactor(A, MATSOLVERSCALAPACK, MAT_FACTOR_LU, &F)); 214 PetscCall(MatLUFactorSymbolic(F, A, 0, 0, NULL)); 215 PetscCall(MatLUFactorNumeric(F, A, NULL)); 216 if (mats_view) PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD)); 217 PetscCall(MatSolve(F, b, x)); 218 PetscCall(MatMatSolve(F, B, X)); 219 PetscCall(MatDestroy(&F)); 220 221 /* Free space */ 222 PetscCall(MatDestroy(&A)); 223 PetscCall(MatDestroy(&Aher)); 224 PetscCall(MatDestroy(&B)); 225 PetscCall(MatDestroy(&C)); 226 PetscCall(MatDestroy(&X)); 227 PetscCall(VecDestroy(&b)); 228 PetscCall(VecDestroy(&c)); 229 PetscCall(VecDestroy(&d)); 230 PetscCall(VecDestroy(&e)); 231 PetscCall(VecDestroy(&x)); 232 PetscCall(PetscRandomDestroy(&rand)); 233 PetscCall(PetscFinalize()); 234 return 0; 235 } 236 237 /*TEST 238 239 build: 240 requires: scalapack 241 242 test: 243 nsize: 2 244 output_file: output/ex245.out 245 246 test: 247 suffix: 2 248 nsize: 6 249 output_file: output/ex245.out 250 251 TEST*/ 252