1 static char help[] = "Tests LU, Cholesky, and QR factorization and MatMatSolve() for a sequential dense matrix. \n\ 2 For MATSEQDENSE matrix, the factorization is just a thin wrapper to LAPACK. \n\ 3 For MATSEQDENSECUDA, it uses cusolverDn routines \n\n"; 4 5 #include <petscmat.h> 6 7 static PetscErrorCode createMatsAndVecs(PetscInt m, PetscInt n, PetscInt nrhs, PetscBool full, Mat *_mat, Mat *_RHS, Mat *_SOLU, Vec *_x, Vec *_y, Vec *_b) 8 { 9 PetscRandom rand; 10 Mat mat, RHS, SOLU; 11 PetscInt rstart, rend; 12 PetscInt cstart, cend; 13 PetscScalar value = 1.0; 14 Vec x, y, b; 15 16 PetscFunctionBegin; 17 /* create multiple vectors RHS and SOLU */ 18 PetscCall(MatCreate(PETSC_COMM_WORLD, &RHS)); 19 PetscCall(MatSetSizes(RHS, PETSC_DECIDE, PETSC_DECIDE, m, nrhs)); 20 PetscCall(MatSetType(RHS, MATDENSE)); 21 PetscCall(MatSetOptionsPrefix(RHS, "rhs_")); 22 PetscCall(MatSetFromOptions(RHS)); 23 PetscCall(MatSeqDenseSetPreallocation(RHS, NULL)); 24 25 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand)); 26 PetscCall(PetscRandomSetFromOptions(rand)); 27 PetscCall(MatSetRandom(RHS, rand)); 28 29 if (m == n) { 30 PetscCall(MatDuplicate(RHS, MAT_DO_NOT_COPY_VALUES, &SOLU)); 31 } else { 32 PetscCall(MatCreate(PETSC_COMM_WORLD, &SOLU)); 33 PetscCall(MatSetSizes(SOLU, PETSC_DECIDE, PETSC_DECIDE, n, nrhs)); 34 PetscCall(MatSetType(SOLU, MATDENSE)); 35 PetscCall(MatSeqDenseSetPreallocation(SOLU, NULL)); 36 } 37 PetscCall(MatSetRandom(SOLU, rand)); 38 39 /* create matrix */ 40 PetscCall(MatCreate(PETSC_COMM_WORLD, &mat)); 41 PetscCall(MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, m, n)); 42 PetscCall(MatSetType(mat, MATDENSE)); 43 PetscCall(MatSetFromOptions(mat)); 44 PetscCall(MatSetUp(mat)); 45 PetscCall(MatGetOwnershipRange(mat, &rstart, &rend)); 46 PetscCall(MatGetOwnershipRangeColumn(mat, &cstart, &cend)); 47 if (!full) { 48 for (PetscInt i = rstart; i < rend; i++) { 49 if (m == n) { 50 value = (PetscReal)i + 1; 51 PetscCall(MatSetValues(mat, 1, &i, 1, &i, &value, INSERT_VALUES)); 52 } else { 53 for (PetscInt j = cstart; j < cend; j++) { 54 value = ((PetscScalar)i + 1.) / (PetscSqr(i - j) + 1.); 55 PetscCall(MatSetValues(mat, 1, &i, 1, &j, &value, INSERT_VALUES)); 56 } 57 } 58 } 59 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 60 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 61 } else { 62 PetscCall(MatSetRandom(mat, rand)); 63 if (m == n) { 64 Mat T; 65 66 PetscCall(MatMatTransposeMult(mat, mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T)); 67 PetscCall(MatDestroy(&mat)); 68 mat = T; 69 } 70 } 71 72 /* create single vectors */ 73 PetscCall(MatCreateVecs(mat, &x, &b)); 74 PetscCall(VecDuplicate(x, &y)); 75 PetscCall(VecSet(x, value)); 76 PetscCall(PetscRandomDestroy(&rand)); 77 *_mat = mat; 78 *_RHS = RHS; 79 *_SOLU = SOLU; 80 *_x = x; 81 *_y = y; 82 *_b = b; 83 PetscFunctionReturn(PETSC_SUCCESS); 84 } 85 86 int main(int argc, char **argv) 87 { 88 Mat mat, F, RHS, SOLU; 89 MatInfo info; 90 PetscInt m = 15, n = 10, i, j, nrhs = 2; 91 Vec x, y, b, ytmp; 92 IS perm; 93 PetscReal norm, tol = PETSC_SMALL; 94 PetscMPIInt size; 95 char solver[64]; 96 PetscBool inplace, full = PETSC_FALSE, ldl = PETSC_TRUE, qr = PETSC_TRUE; 97 98 PetscFunctionBeginUser; 99 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 100 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 101 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 102 PetscCall(PetscStrncpy(solver, MATSOLVERPETSC, sizeof(solver))); 103 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 104 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 105 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL)); 106 PetscCall(PetscOptionsGetBool(NULL, NULL, "-ldl", &ldl, NULL)); 107 PetscCall(PetscOptionsGetBool(NULL, NULL, "-qr", &qr, NULL)); 108 PetscCall(PetscOptionsGetBool(NULL, NULL, "-full", &full, NULL)); 109 PetscCall(PetscOptionsGetString(NULL, NULL, "-solver_type", solver, sizeof(solver), NULL)); 110 111 PetscCall(createMatsAndVecs(n, n, nrhs, full, &mat, &RHS, &SOLU, &x, &y, &b)); 112 PetscCall(VecDuplicate(y, &ytmp)); 113 114 /* Only SeqDense* support in-place factorizations and NULL permutations */ 115 PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQDENSE, &inplace)); 116 PetscCall(MatGetLocalSize(mat, &i, NULL)); 117 PetscCall(MatGetOwnershipRange(mat, &j, NULL)); 118 PetscCall(ISCreateStride(PETSC_COMM_WORLD, i, j, 1, &perm)); 119 120 PetscCall(MatGetInfo(mat, MAT_LOCAL, &info)); 121 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "matrix nonzeros = %" PetscInt_FMT ", allocated nonzeros = %" PetscInt_FMT "\n", (PetscInt)info.nz_used, (PetscInt)info.nz_allocated)); 122 PetscCall(MatMult(mat, x, b)); 123 124 /* Cholesky factorization - perm and factinfo are ignored by LAPACK */ 125 /* in-place Cholesky */ 126 if (inplace) { 127 Mat RHS2; 128 129 if (!ldl) PetscCall(MatSetOption(mat, MAT_SPD, PETSC_TRUE)); 130 PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &F)); 131 PetscCall(MatCholeskyFactor(F, perm, 0)); 132 PetscCall(MatSolve(F, b, y)); 133 PetscCall(VecAXPY(y, -1.0, x)); 134 PetscCall(VecNorm(y, NORM_2, &norm)); 135 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for in-place Cholesky %g\n", (double)norm)); 136 137 PetscCall(MatMatSolve(F, RHS, SOLU)); 138 PetscCall(MatMatMult(mat, SOLU, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &RHS2)); 139 PetscCall(MatAXPY(RHS, -1.0, RHS2, SAME_NONZERO_PATTERN)); 140 PetscCall(MatNorm(RHS, NORM_FROBENIUS, &norm)); 141 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: Norm of residual for in-place Cholesky (MatMatSolve) %g\n", (double)norm)); 142 PetscCall(MatDestroy(&F)); 143 PetscCall(MatDestroy(&RHS2)); 144 } 145 146 /* out-of-place Cholesky */ 147 if (!ldl) PetscCall(MatSetOption(mat, MAT_SPD, PETSC_TRUE)); 148 PetscCall(MatGetFactor(mat, solver, MAT_FACTOR_CHOLESKY, &F)); 149 PetscCall(MatCholeskyFactorSymbolic(F, mat, perm, 0)); 150 PetscCall(MatCholeskyFactorNumeric(F, mat, 0)); 151 PetscCall(MatSolve(F, b, y)); 152 PetscCall(VecAXPY(y, -1.0, x)); 153 PetscCall(VecNorm(y, NORM_2, &norm)); 154 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for out-of-place Cholesky %g\n", (double)norm)); 155 PetscCall(MatDestroy(&F)); 156 157 /* LU factorization - perms and factinfo are ignored by LAPACK */ 158 i = n - 1; 159 PetscCall(MatZeroRows(mat, 1, &i, -1.0, NULL, NULL)); 160 PetscCall(MatMult(mat, x, b)); 161 162 /* in-place LU */ 163 if (inplace) { 164 Mat RHS2; 165 166 PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &F)); 167 PetscCall(MatLUFactor(F, perm, perm, 0)); 168 PetscCall(MatSolve(F, b, y)); 169 PetscCall(VecAXPY(y, -1.0, x)); 170 PetscCall(VecNorm(y, NORM_2, &norm)); 171 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for in-place LU %g\n", (double)norm)); 172 PetscCall(MatMatSolve(F, RHS, SOLU)); 173 PetscCall(MatMatMult(mat, SOLU, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &RHS2)); 174 PetscCall(MatAXPY(RHS, -1.0, RHS2, SAME_NONZERO_PATTERN)); 175 PetscCall(MatNorm(RHS, NORM_FROBENIUS, &norm)); 176 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: Norm of residual for in-place LU (MatMatSolve) %g\n", (double)norm)); 177 PetscCall(MatDestroy(&F)); 178 PetscCall(MatDestroy(&RHS2)); 179 } 180 181 /* out-of-place LU */ 182 PetscCall(MatGetFactor(mat, solver, MAT_FACTOR_LU, &F)); 183 PetscCall(MatLUFactorSymbolic(F, mat, perm, perm, 0)); 184 PetscCall(MatLUFactorNumeric(F, mat, 0)); 185 PetscCall(MatSolve(F, b, y)); 186 PetscCall(VecAXPY(y, -1.0, x)); 187 PetscCall(VecNorm(y, NORM_2, &norm)); 188 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for out-of-place LU %g\n", (double)norm)); 189 190 /* free space */ 191 PetscCall(ISDestroy(&perm)); 192 PetscCall(MatDestroy(&F)); 193 PetscCall(MatDestroy(&mat)); 194 PetscCall(MatDestroy(&RHS)); 195 PetscCall(MatDestroy(&SOLU)); 196 PetscCall(VecDestroy(&x)); 197 PetscCall(VecDestroy(&b)); 198 PetscCall(VecDestroy(&y)); 199 PetscCall(VecDestroy(&ytmp)); 200 201 if (qr) { 202 /* setup rectangular */ 203 PetscCall(createMatsAndVecs(m, n, nrhs, full, &mat, &RHS, &SOLU, &x, &y, &b)); 204 PetscCall(VecDuplicate(y, &ytmp)); 205 206 /* QR factorization - perms and factinfo are ignored by LAPACK */ 207 PetscCall(MatMult(mat, x, b)); 208 209 /* in-place QR */ 210 if (inplace) { 211 Mat SOLU2; 212 213 PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &F)); 214 PetscCall(MatQRFactor(F, NULL, 0)); 215 PetscCall(MatSolve(F, b, y)); 216 PetscCall(VecAXPY(y, -1.0, x)); 217 PetscCall(VecNorm(y, NORM_2, &norm)); 218 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for in-place QR %g\n", (double)norm)); 219 PetscCall(MatMatMult(mat, SOLU, MAT_REUSE_MATRIX, PETSC_DEFAULT, &RHS)); 220 PetscCall(MatDuplicate(SOLU, MAT_DO_NOT_COPY_VALUES, &SOLU2)); 221 PetscCall(MatMatSolve(F, RHS, SOLU2)); 222 PetscCall(MatAXPY(SOLU2, -1.0, SOLU, SAME_NONZERO_PATTERN)); 223 PetscCall(MatNorm(SOLU2, NORM_FROBENIUS, &norm)); 224 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: Norm of error for in-place QR (MatMatSolve) %g\n", (double)norm)); 225 PetscCall(MatDestroy(&F)); 226 PetscCall(MatDestroy(&SOLU2)); 227 } 228 229 /* out-of-place QR */ 230 PetscCall(MatGetFactor(mat, solver, MAT_FACTOR_QR, &F)); 231 PetscCall(MatQRFactorSymbolic(F, mat, NULL, NULL)); 232 PetscCall(MatQRFactorNumeric(F, mat, NULL)); 233 PetscCall(MatSolve(F, b, y)); 234 PetscCall(VecAXPY(y, -1.0, x)); 235 PetscCall(VecNorm(y, NORM_2, &norm)); 236 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for out-of-place QR %g\n", (double)norm)); 237 238 if (m == n) { 239 /* out-of-place MatSolveTranspose */ 240 PetscCall(MatMultTranspose(mat, x, b)); 241 PetscCall(MatSolveTranspose(F, b, y)); 242 PetscCall(VecAXPY(y, -1.0, x)); 243 PetscCall(VecNorm(y, NORM_2, &norm)); 244 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for out-of-place QR %g\n", (double)norm)); 245 } 246 247 /* free space */ 248 PetscCall(MatDestroy(&F)); 249 PetscCall(MatDestroy(&mat)); 250 PetscCall(MatDestroy(&RHS)); 251 PetscCall(MatDestroy(&SOLU)); 252 PetscCall(VecDestroy(&x)); 253 PetscCall(VecDestroy(&b)); 254 PetscCall(VecDestroy(&y)); 255 PetscCall(VecDestroy(&ytmp)); 256 } 257 PetscCall(PetscFinalize()); 258 return 0; 259 } 260 261 /*TEST 262 263 test: 264 265 test: 266 requires: cuda 267 suffix: seqdensecuda 268 args: -mat_type seqdensecuda -rhs_mat_type seqdensecuda -ldl 0 -solver_type {{petsc cuda}} 269 output_file: output/ex1_1.out 270 271 test: 272 requires: cuda 273 suffix: seqdensecuda_2 274 args: -ldl 0 -solver_type cuda 275 output_file: output/ex1_1.out 276 277 test: 278 requires: cuda 279 suffix: seqdensecuda_seqaijcusparse 280 args: -mat_type seqaijcusparse -rhs_mat_type seqdensecuda -qr 0 281 output_file: output/ex1_2.out 282 283 test: 284 requires: cuda viennacl 285 suffix: seqdensecuda_seqaijviennacl 286 args: -mat_type seqaijviennacl -rhs_mat_type seqdensecuda -qr 0 287 output_file: output/ex1_2.out 288 289 test: 290 suffix: 4 291 args: -m 10 -n 10 292 output_file: output/ex1_1.out 293 294 TEST*/ 295