1 static char help[] = "Tests MatSolve() and MatMatSolve() (interface to superlu_dist, mumps and mkl_pardiso).\n\ 2 Example: mpiexec -n <np> ./ex125 -f <matrix binary file> -nrhs 4 \n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc, char **args) 7 { 8 Mat A, RHS = NULL, RHS1 = NULL, C, F, X; 9 Vec u, x, b; 10 PetscMPIInt size; 11 PetscInt m, n, nfact, nsolve, nrhs, ipack = 0; 12 PetscReal norm, tol = 1.e-10; 13 IS perm, iperm; 14 MatFactorInfo info; 15 PetscRandom rand; 16 PetscBool flg, testMatSolve = PETSC_TRUE, testMatMatSolve = PETSC_TRUE, testMatMatSolveTranspose = PETSC_TRUE; 17 PetscBool chol = PETSC_FALSE, view = PETSC_FALSE, matsolvexx = PETSC_FALSE; 18 #if defined(PETSC_HAVE_MUMPS) 19 PetscBool test_mumps_opts = PETSC_FALSE; 20 #endif 21 PetscViewer fd; /* viewer */ 22 char file[PETSC_MAX_PATH_LEN]; /* input file name */ 23 24 PetscFunctionBeginUser; 25 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 26 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 27 28 /* Determine file from which we read the matrix A */ 29 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg)); 30 if (flg) { /* Load matrix A */ 31 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd)); 32 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 33 PetscCall(MatSetFromOptions(A)); 34 PetscCall(MatLoad(A, fd)); 35 PetscCall(PetscViewerDestroy(&fd)); 36 } else { 37 n = 13; 38 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 39 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 40 PetscCall(MatSetType(A, MATAIJ)); 41 PetscCall(MatSetFromOptions(A)); 42 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n)); 43 PetscCall(MatSetUp(A)); 44 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 45 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 46 PetscCall(MatShift(A, 1.0)); 47 } 48 PetscCall(MatGetLocalSize(A, &m, &n)); 49 PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n); 50 51 /* if A is symmetric, set its flag -- required by MatGetInertia() */ 52 PetscCall(MatIsSymmetric(A, 0.0, &flg)); 53 54 PetscCall(MatViewFromOptions(A, NULL, "-A_view")); 55 56 /* Create dense matrix C and X; C holds true solution with identical columns */ 57 nrhs = 2; 58 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL)); 59 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ex125: nrhs %" PetscInt_FMT "\n", nrhs)); 60 PetscCall(MatCreate(PETSC_COMM_WORLD, &C)); 61 PetscCall(MatSetOptionsPrefix(C, "rhs_")); 62 PetscCall(MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs)); 63 PetscCall(MatSetType(C, MATDENSE)); 64 PetscCall(MatSetFromOptions(C)); 65 PetscCall(MatSetUp(C)); 66 67 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_factor", &view, NULL)); 68 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmatsolve", &testMatMatSolve, NULL)); 69 PetscCall(PetscOptionsGetBool(NULL, NULL, "-cholesky", &chol, NULL)); 70 #if defined(PETSC_HAVE_MUMPS) 71 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_mumps_opts", &test_mumps_opts, NULL)); 72 #endif 73 74 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand)); 75 PetscCall(PetscRandomSetFromOptions(rand)); 76 PetscCall(MatSetRandom(C, rand)); 77 PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &X)); 78 79 /* Create vectors */ 80 PetscCall(MatCreateVecs(A, &x, &b)); 81 PetscCall(VecDuplicate(x, &u)); /* save the true solution */ 82 83 /* Test Factorization */ 84 PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm)); 85 86 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_solver_type", &ipack, NULL)); 87 switch (ipack) { 88 #if defined(PETSC_HAVE_SUPERLU) 89 case 0: 90 PetscCheck(!chol, PETSC_COMM_WORLD, PETSC_ERR_SUP, "SuperLU does not provide Cholesky!"); 91 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SUPERLU LU:\n")); 92 PetscCall(MatGetFactor(A, MATSOLVERSUPERLU, MAT_FACTOR_LU, &F)); 93 matsolvexx = PETSC_TRUE; 94 break; 95 #endif 96 #if defined(PETSC_HAVE_SUPERLU_DIST) 97 case 1: 98 PetscCheck(!chol, PETSC_COMM_WORLD, PETSC_ERR_SUP, "SuperLU does not provide Cholesky!"); 99 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SUPERLU_DIST LU:\n")); 100 PetscCall(MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_LU, &F)); 101 matsolvexx = PETSC_TRUE; 102 break; 103 #endif 104 #if defined(PETSC_HAVE_MUMPS) 105 case 2: 106 if (chol) { 107 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MUMPS CHOLESKY:\n")); 108 PetscCall(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_CHOLESKY, &F)); 109 } else { 110 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MUMPS LU:\n")); 111 PetscCall(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_LU, &F)); 112 } 113 matsolvexx = PETSC_TRUE; 114 if (test_mumps_opts) { 115 /* test mumps options */ 116 PetscInt icntl; 117 PetscReal cntl; 118 119 icntl = 2; /* sequential matrix ordering */ 120 PetscCall(MatMumpsSetIcntl(F, 7, icntl)); 121 122 cntl = 1.e-6; /* threshold for row pivot detection */ 123 PetscCall(MatMumpsSetIcntl(F, 24, 1)); 124 PetscCall(MatMumpsSetCntl(F, 3, cntl)); 125 } 126 break; 127 #endif 128 #if defined(PETSC_HAVE_MKL_PARDISO) 129 case 3: 130 if (chol) { 131 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO CHOLESKY:\n")); 132 PetscCall(MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_CHOLESKY, &F)); 133 } else { 134 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO LU:\n")); 135 PetscCall(MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_LU, &F)); 136 } 137 break; 138 #endif 139 #if defined(PETSC_HAVE_CUDA) 140 case 4: 141 if (chol) { 142 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE CHOLESKY:\n")); 143 PetscCall(MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_CHOLESKY, &F)); 144 } else { 145 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE LU:\n")); 146 PetscCall(MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_LU, &F)); 147 } 148 break; 149 #endif 150 default: 151 if (chol) { 152 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PETSC CHOLESKY:\n")); 153 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &F)); 154 } else { 155 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PETSC LU:\n")); 156 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F)); 157 } 158 matsolvexx = PETSC_TRUE; 159 } 160 161 PetscCall(MatFactorInfoInitialize(&info)); 162 info.fill = 5.0; 163 info.shifttype = (PetscReal)MAT_SHIFT_NONE; 164 if (chol) { 165 PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info)); 166 } else { 167 PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info)); 168 } 169 170 for (nfact = 0; nfact < 2; nfact++) { 171 if (chol) { 172 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the CHOLESKY numfactorization \n", nfact)); 173 PetscCall(MatCholeskyFactorNumeric(F, A, &info)); 174 } else { 175 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the LU numfactorization \n", nfact)); 176 PetscCall(MatLUFactorNumeric(F, A, &info)); 177 } 178 if (view) { 179 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO)); 180 PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD)); 181 PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD)); 182 view = PETSC_FALSE; 183 } 184 185 #if defined(PETSC_HAVE_SUPERLU_DIST) 186 if (ipack == 1) { /* Test MatSuperluDistGetDiagU() 187 -- input: matrix factor F; output: main diagonal of matrix U on all processes */ 188 PetscInt M; 189 PetscScalar *diag; 190 #if !defined(PETSC_USE_COMPLEX) 191 PetscInt nneg, nzero, npos; 192 #endif 193 194 PetscCall(MatGetSize(F, &M, NULL)); 195 PetscCall(PetscMalloc1(M, &diag)); 196 PetscCall(MatSuperluDistGetDiagU(F, diag)); 197 PetscCall(PetscFree(diag)); 198 199 #if !defined(PETSC_USE_COMPLEX) 200 /* Test MatGetInertia() */ 201 PetscCall(MatGetInertia(F, &nneg, &nzero, &npos)); 202 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos)); 203 #endif 204 } 205 #endif 206 207 #if defined(PETSC_HAVE_MUMPS) 208 /* mumps interface allows repeated call of MatCholeskyFactorSymbolic(), while the succession calls do nothing */ 209 if (ipack == 2) { 210 if (chol) { 211 PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info)); 212 PetscCall(MatCholeskyFactorNumeric(F, A, &info)); 213 } else { 214 PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info)); 215 PetscCall(MatLUFactorNumeric(F, A, &info)); 216 } 217 } 218 #endif 219 220 /* Test MatMatSolve(), A X = B, where B can be dense or sparse */ 221 if (testMatMatSolve) { 222 if (!nfact) { 223 PetscCall(MatMatMult(A, C, MAT_INITIAL_MATRIX, 2.0, &RHS)); 224 } else { 225 PetscCall(MatMatMult(A, C, MAT_REUSE_MATRIX, 2.0, &RHS)); 226 } 227 for (nsolve = 0; nsolve < 2; nsolve++) { 228 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatMatSolve \n", nsolve)); 229 PetscCall(MatMatSolve(F, RHS, X)); 230 231 /* Check the error */ 232 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN)); 233 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm)); 234 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve)); 235 } 236 237 if (matsolvexx) { 238 /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */ 239 PetscCall(MatCopy(RHS, X, SAME_NONZERO_PATTERN)); 240 PetscCall(MatMatSolve(F, X, X)); 241 /* Check the error */ 242 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN)); 243 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm)); 244 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMatSolve(F,RHS,RHS): Norm of error %g\n", (double)norm)); 245 } 246 247 if (ipack == 2 && size == 1) { 248 Mat spRHS, spRHST, RHST; 249 250 PetscCall(MatTranspose(RHS, MAT_INITIAL_MATRIX, &RHST)); 251 PetscCall(MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST)); 252 PetscCall(MatCreateTranspose(spRHST, &spRHS)); 253 for (nsolve = 0; nsolve < 2; nsolve++) { 254 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the sparse MatMatSolve \n", nsolve)); 255 PetscCall(MatMatSolve(F, spRHS, X)); 256 257 /* Check the error */ 258 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN)); 259 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm)); 260 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the sparse MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve)); 261 } 262 PetscCall(MatDestroy(&spRHST)); 263 PetscCall(MatDestroy(&spRHS)); 264 PetscCall(MatDestroy(&RHST)); 265 } 266 } 267 268 /* Test testMatMatSolveTranspose(), A^T X = B, where B can be dense or sparse */ 269 if (testMatMatSolveTranspose) { 270 if (!nfact) { 271 PetscCall(MatTransposeMatMult(A, C, MAT_INITIAL_MATRIX, 2.0, &RHS1)); 272 } else { 273 PetscCall(MatTransposeMatMult(A, C, MAT_REUSE_MATRIX, 2.0, &RHS1)); 274 } 275 276 for (nsolve = 0; nsolve < 2; nsolve++) { 277 PetscCall(MatMatSolveTranspose(F, RHS1, X)); 278 279 /* Check the error */ 280 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN)); 281 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm)); 282 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve)); 283 } 284 285 if (ipack == 2 && size == 1) { 286 Mat spRHS, spRHST, RHST; 287 288 PetscCall(MatTranspose(RHS1, MAT_INITIAL_MATRIX, &RHST)); 289 PetscCall(MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST)); 290 PetscCall(MatCreateTranspose(spRHST, &spRHS)); 291 for (nsolve = 0; nsolve < 2; nsolve++) { 292 PetscCall(MatMatSolveTranspose(F, spRHS, X)); 293 294 /* Check the error */ 295 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN)); 296 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm)); 297 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the sparse MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve)); 298 } 299 PetscCall(MatDestroy(&spRHST)); 300 PetscCall(MatDestroy(&spRHS)); 301 PetscCall(MatDestroy(&RHST)); 302 } 303 } 304 305 /* Test MatSolve() */ 306 if (testMatSolve) { 307 for (nsolve = 0; nsolve < 2; nsolve++) { 308 PetscCall(VecSetRandom(x, rand)); 309 PetscCall(VecCopy(x, u)); 310 PetscCall(MatMult(A, x, b)); 311 312 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatSolve \n", nsolve)); 313 PetscCall(MatSolve(F, b, x)); 314 315 /* Check the error */ 316 PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */ 317 PetscCall(VecNorm(u, NORM_2, &norm)); 318 if (norm > tol) { 319 PetscReal resi; 320 PetscCall(MatMult(A, x, u)); /* u = A*x */ 321 PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */ 322 PetscCall(VecNorm(u, NORM_2, &resi)); 323 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolve: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact)); 324 } 325 } 326 } 327 } 328 329 /* Free data structures */ 330 PetscCall(MatDestroy(&A)); 331 PetscCall(MatDestroy(&C)); 332 PetscCall(MatDestroy(&F)); 333 PetscCall(MatDestroy(&X)); 334 PetscCall(MatDestroy(&RHS)); 335 PetscCall(MatDestroy(&RHS1)); 336 337 PetscCall(PetscRandomDestroy(&rand)); 338 PetscCall(ISDestroy(&perm)); 339 PetscCall(ISDestroy(&iperm)); 340 PetscCall(VecDestroy(&x)); 341 PetscCall(VecDestroy(&b)); 342 PetscCall(VecDestroy(&u)); 343 PetscCall(PetscFinalize()); 344 return 0; 345 } 346 347 /*TEST 348 349 test: 350 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 351 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 10 352 output_file: output/ex125.out 353 354 test: 355 suffix: 2 356 args: -mat_solver_type 10 357 output_file: output/ex125.out 358 359 test: 360 suffix: mkl_pardiso 361 requires: mkl_pardiso datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 362 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 3 363 364 test: 365 suffix: mkl_pardiso_2 366 requires: mkl_pardiso 367 args: -mat_solver_type 3 368 output_file: output/ex125_mkl_pardiso.out 369 370 test: 371 suffix: mumps 372 requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 373 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2 374 output_file: output/ex125_mumps_seq.out 375 376 test: 377 suffix: mumps_2 378 nsize: 3 379 requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 380 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2 381 output_file: output/ex125_mumps_par.out 382 383 test: 384 suffix: mumps_3 385 requires: mumps 386 args: -mat_solver_type 2 387 output_file: output/ex125_mumps_seq.out 388 389 test: 390 suffix: mumps_4 391 nsize: 3 392 requires: mumps 393 args: -mat_solver_type 2 394 output_file: output/ex125_mumps_par.out 395 396 test: 397 suffix: mumps_5 398 nsize: 3 399 requires: mumps 400 args: -mat_solver_type 2 -cholesky 401 output_file: output/ex125_mumps_par_cholesky.out 402 403 test: 404 suffix: superlu_dist 405 nsize: {{1 3}} 406 requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist 407 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM 408 409 test: 410 suffix: superlu_dist_2 411 nsize: {{1 3}} 412 requires: superlu_dist !complex 413 args: -n 36 -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM 414 output_file: output/ex125_superlu_dist.out 415 416 test: 417 suffix: superlu_dist_complex 418 nsize: 3 419 requires: datafilespath double superlu_dist complex !defined(PETSC_USE_64BIT_INDICES) 420 args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type 1 421 output_file: output/ex125_superlu_dist_complex.out 422 423 test: 424 suffix: superlu_dist_complex_2 425 nsize: 3 426 requires: superlu_dist complex 427 args: -mat_solver_type 1 428 output_file: output/ex125_superlu_dist_complex.out 429 430 test: 431 suffix: cusparse 432 requires: cuda datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 433 #todo: fix the bug with cholesky 434 #args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type 4 -cholesky {{0 1}separate output} 435 args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type 4 -cholesky {{0}separate output} 436 437 test: 438 suffix: cusparse_2 439 requires: cuda 440 args: -mat_type aijcusparse -mat_solver_type 4 -cholesky {{0 1}separate output} 441 442 TEST*/ 443