1 2 static char help[] = "Tests LU, Cholesky, and QR factorization and MatMatSolve() for a sequential dense matrix. \n\ 3 For MATSEQDENSE matrix, the factorization is just a thin wrapper to LAPACK. \n\ 4 For MATSEQDENSECUDA, it uses cusolverDn routines \n\n"; 5 6 #include <petscmat.h> 7 8 static PetscErrorCode createMatsAndVecs(PetscInt m, PetscInt n, PetscInt nrhs, PetscBool full, Mat *_mat, Mat *_RHS, Mat *_SOLU, Vec *_x, Vec *_y, Vec *_b) { 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(0); 84 } 85 86 int main(int argc, char **argv) { 87 Mat mat, F, RHS, SOLU; 88 MatInfo info; 89 PetscInt m = 15, n = 10, i, j, nrhs = 2; 90 Vec x, y, b, ytmp; 91 IS perm; 92 PetscReal norm, tol = PETSC_SMALL; 93 PetscMPIInt size; 94 char solver[64]; 95 PetscBool inplace, full = PETSC_FALSE, ldl = PETSC_TRUE, qr = PETSC_TRUE; 96 97 PetscFunctionBeginUser; 98 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 99 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 100 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 101 PetscCall(PetscStrcpy(solver, "petsc")); 102 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 103 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 104 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL)); 105 PetscCall(PetscOptionsGetBool(NULL, NULL, "-ldl", &ldl, NULL)); 106 PetscCall(PetscOptionsGetBool(NULL, NULL, "-qr", &qr, NULL)); 107 PetscCall(PetscOptionsGetBool(NULL, NULL, "-full", &full, NULL)); 108 PetscCall(PetscOptionsGetString(NULL, NULL, "-solver_type", solver, sizeof(solver), NULL)); 109 110 PetscCall(createMatsAndVecs(n, n, nrhs, full, &mat, &RHS, &SOLU, &x, &y, &b)); 111 PetscCall(VecDuplicate(y, &ytmp)); 112 113 /* Only SeqDense* support in-place factorizations and NULL permutations */ 114 PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQDENSE, &inplace)); 115 PetscCall(MatGetLocalSize(mat, &i, NULL)); 116 PetscCall(MatGetOwnershipRange(mat, &j, NULL)); 117 PetscCall(ISCreateStride(PETSC_COMM_WORLD, i, j, 1, &perm)); 118 119 PetscCall(MatGetInfo(mat, MAT_LOCAL, &info)); 120 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "matrix nonzeros = %" PetscInt_FMT ", allocated nonzeros = %" PetscInt_FMT "\n", (PetscInt)info.nz_used, (PetscInt)info.nz_allocated)); 121 PetscCall(MatMult(mat, x, b)); 122 123 /* Cholesky factorization - perm and factinfo are ignored by LAPACK */ 124 /* in-place Cholesky */ 125 if (inplace) { 126 Mat RHS2; 127 128 PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &F)); 129 if (!ldl) PetscCall(MatSetOption(F, MAT_SPD, PETSC_TRUE)); 130 PetscCall(MatCholeskyFactor(F, perm, 0)); 131 PetscCall(MatSolve(F, b, y)); 132 PetscCall(VecAXPY(y, -1.0, x)); 133 PetscCall(VecNorm(y, NORM_2, &norm)); 134 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for in-place Cholesky %g\n", (double)norm)); } 135 136 PetscCall(MatMatSolve(F, RHS, SOLU)); 137 PetscCall(MatMatMult(mat, SOLU, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &RHS2)); 138 PetscCall(MatAXPY(RHS, -1.0, RHS2, SAME_NONZERO_PATTERN)); 139 PetscCall(MatNorm(RHS, NORM_FROBENIUS, &norm)); 140 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: Norm of residual for in-place Cholesky (MatMatSolve) %g\n", (double)norm)); } 141 PetscCall(MatDestroy(&F)); 142 PetscCall(MatDestroy(&RHS2)); 143 } 144 145 /* out-of-place Cholesky */ 146 PetscCall(MatGetFactor(mat, solver, MAT_FACTOR_CHOLESKY, &F)); 147 if (!ldl) PetscCall(MatSetOption(F, MAT_SPD, PETSC_TRUE)); 148 PetscCall(MatCholeskyFactorSymbolic(F, mat, perm, 0)); 149 PetscCall(MatCholeskyFactorNumeric(F, mat, 0)); 150 PetscCall(MatSolve(F, b, y)); 151 PetscCall(VecAXPY(y, -1.0, x)); 152 PetscCall(VecNorm(y, NORM_2, &norm)); 153 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for out-of-place Cholesky %g\n", (double)norm)); } 154 PetscCall(MatDestroy(&F)); 155 156 /* LU factorization - perms and factinfo are ignored by LAPACK */ 157 i = n - 1; 158 PetscCall(MatZeroRows(mat, 1, &i, -1.0, NULL, NULL)); 159 PetscCall(MatMult(mat, x, b)); 160 161 /* in-place LU */ 162 if (inplace) { 163 Mat RHS2; 164 165 PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &F)); 166 PetscCall(MatLUFactor(F, perm, perm, 0)); 167 PetscCall(MatSolve(F, b, y)); 168 PetscCall(VecAXPY(y, -1.0, x)); 169 PetscCall(VecNorm(y, NORM_2, &norm)); 170 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for in-place LU %g\n", (double)norm)); } 171 PetscCall(MatMatSolve(F, RHS, SOLU)); 172 PetscCall(MatMatMult(mat, SOLU, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &RHS2)); 173 PetscCall(MatAXPY(RHS, -1.0, RHS2, SAME_NONZERO_PATTERN)); 174 PetscCall(MatNorm(RHS, NORM_FROBENIUS, &norm)); 175 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: Norm of residual for in-place LU (MatMatSolve) %g\n", (double)norm)); } 176 PetscCall(MatDestroy(&F)); 177 PetscCall(MatDestroy(&RHS2)); 178 } 179 180 /* out-of-place LU */ 181 PetscCall(MatGetFactor(mat, solver, MAT_FACTOR_LU, &F)); 182 PetscCall(MatLUFactorSymbolic(F, mat, perm, perm, 0)); 183 PetscCall(MatLUFactorNumeric(F, mat, 0)); 184 PetscCall(MatSolve(F, b, y)); 185 PetscCall(VecAXPY(y, -1.0, x)); 186 PetscCall(VecNorm(y, NORM_2, &norm)); 187 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for out-of-place LU %g\n", (double)norm)); } 188 189 /* free space */ 190 PetscCall(ISDestroy(&perm)); 191 PetscCall(MatDestroy(&F)); 192 PetscCall(MatDestroy(&mat)); 193 PetscCall(MatDestroy(&RHS)); 194 PetscCall(MatDestroy(&SOLU)); 195 PetscCall(VecDestroy(&x)); 196 PetscCall(VecDestroy(&b)); 197 PetscCall(VecDestroy(&y)); 198 PetscCall(VecDestroy(&ytmp)); 199 200 if (qr) { 201 /* setup rectanglar */ 202 PetscCall(createMatsAndVecs(m, n, nrhs, full, &mat, &RHS, &SOLU, &x, &y, &b)); 203 PetscCall(VecDuplicate(y, &ytmp)); 204 205 /* QR factorization - perms and factinfo are ignored by LAPACK */ 206 PetscCall(MatMult(mat, x, b)); 207 208 /* in-place QR */ 209 if (inplace) { 210 Mat SOLU2; 211 212 PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &F)); 213 PetscCall(MatQRFactor(F, NULL, 0)); 214 PetscCall(MatSolve(F, b, y)); 215 PetscCall(VecAXPY(y, -1.0, x)); 216 PetscCall(VecNorm(y, NORM_2, &norm)); 217 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for in-place QR %g\n", (double)norm)); } 218 PetscCall(MatMatMult(mat, SOLU, MAT_REUSE_MATRIX, PETSC_DEFAULT, &RHS)); 219 PetscCall(MatDuplicate(SOLU, MAT_DO_NOT_COPY_VALUES, &SOLU2)); 220 PetscCall(MatMatSolve(F, RHS, SOLU2)); 221 PetscCall(MatAXPY(SOLU2, -1.0, SOLU, SAME_NONZERO_PATTERN)); 222 PetscCall(MatNorm(SOLU2, NORM_FROBENIUS, &norm)); 223 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: Norm of error for in-place QR (MatMatSolve) %g\n", (double)norm)); } 224 PetscCall(MatDestroy(&F)); 225 PetscCall(MatDestroy(&SOLU2)); 226 } 227 228 /* out-of-place QR */ 229 PetscCall(MatGetFactor(mat, solver, MAT_FACTOR_QR, &F)); 230 PetscCall(MatQRFactorSymbolic(F, mat, NULL, NULL)); 231 PetscCall(MatQRFactorNumeric(F, mat, NULL)); 232 PetscCall(MatSolve(F, b, y)); 233 PetscCall(VecAXPY(y, -1.0, x)); 234 PetscCall(VecNorm(y, NORM_2, &norm)); 235 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for out-of-place QR %g\n", (double)norm)); } 236 237 if (m == n) { 238 /* out-of-place MatSolveTranspose */ 239 PetscCall(MatMultTranspose(mat, x, b)); 240 PetscCall(MatSolveTranspose(F, b, y)); 241 PetscCall(VecAXPY(y, -1.0, x)); 242 PetscCall(VecNorm(y, NORM_2, &norm)); 243 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Norm of error for out-of-place QR %g\n", (double)norm)); } 244 } 245 246 /* free space */ 247 PetscCall(MatDestroy(&F)); 248 PetscCall(MatDestroy(&mat)); 249 PetscCall(MatDestroy(&RHS)); 250 PetscCall(MatDestroy(&SOLU)); 251 PetscCall(VecDestroy(&x)); 252 PetscCall(VecDestroy(&b)); 253 PetscCall(VecDestroy(&y)); 254 PetscCall(VecDestroy(&ytmp)); 255 } 256 PetscCall(PetscFinalize()); 257 return 0; 258 } 259 260 /*TEST 261 262 test: 263 264 test: 265 requires: cuda 266 suffix: seqdensecuda 267 args: -mat_type seqdensecuda -rhs_mat_type seqdensecuda -ldl 0 -solver_type {{petsc cuda}} 268 output_file: output/ex1_1.out 269 270 test: 271 requires: cuda 272 suffix: seqdensecuda_2 273 args: -ldl 0 -solver_type cuda 274 output_file: output/ex1_1.out 275 276 test: 277 requires: cuda 278 suffix: seqdensecuda_seqaijcusparse 279 args: -mat_type seqaijcusparse -rhs_mat_type seqdensecuda -qr 0 280 output_file: output/ex1_2.out 281 282 test: 283 requires: cuda viennacl 284 suffix: seqdensecuda_seqaijviennacl 285 args: -mat_type seqaijviennacl -rhs_mat_type seqdensecuda -qr 0 286 output_file: output/ex1_2.out 287 288 test: 289 suffix: 4 290 args: -m 10 -n 10 291 output_file: output/ex1_1.out 292 293 TEST*/ 294