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