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 -mat_solver_type <>\n\n"; 3 4 /* 5 -mat_solver_type: 6 superlu 7 superlu_dist 8 mumps 9 mkl_pardiso 10 cusparse 11 petsc 12 */ 13 14 #include <petscmat.h> 15 16 PetscErrorCode CreateRandom(PetscInt n, PetscInt m, Mat *A) 17 { 18 PetscFunctionBeginUser; 19 PetscCall(MatCreate(PETSC_COMM_WORLD, A)); 20 PetscCall(MatSetType(*A, MATAIJ)); 21 PetscCall(MatSetFromOptions(*A)); 22 PetscCall(MatSetSizes(*A, PETSC_DECIDE, PETSC_DECIDE, n, m)); 23 PetscCall(MatSeqAIJSetPreallocation(*A, 5, NULL)); 24 PetscCall(MatMPIAIJSetPreallocation(*A, 5, NULL, 5, NULL)); 25 PetscCall(MatSetRandom(*A, NULL)); 26 PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY)); 27 PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY)); 28 PetscFunctionReturn(PETSC_SUCCESS); 29 } 30 31 PetscErrorCode CreateIdentity(PetscInt n, Mat *A) 32 { 33 PetscFunctionBeginUser; 34 PetscCall(MatCreate(PETSC_COMM_WORLD, A)); 35 PetscCall(MatSetType(*A, MATAIJ)); 36 PetscCall(MatSetFromOptions(*A)); 37 PetscCall(MatSetSizes(*A, PETSC_DECIDE, PETSC_DECIDE, n, n)); 38 PetscCall(MatSetUp(*A)); 39 PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY)); 40 PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY)); 41 PetscCall(MatShift(*A, 1.0)); 42 PetscFunctionReturn(PETSC_SUCCESS); 43 } 44 45 int main(int argc, char **args) 46 { 47 Mat A, Ae, RHS = NULL, RHS1 = NULL, C, F, X; 48 Vec u, x, b; 49 PetscMPIInt size; 50 PetscInt m, n, nfact, nsolve, nrhs, ipack = 5; 51 PetscReal norm, tol = 10 * PETSC_SQRT_MACHINE_EPSILON; 52 IS perm = NULL, iperm = NULL; 53 MatFactorInfo info; 54 PetscRandom rand; 55 PetscBool flg, symm, testMatSolve = PETSC_TRUE, testMatMatSolve = PETSC_TRUE, testMatMatSolveTranspose = PETSC_TRUE, testMatSolveTranspose = PETSC_TRUE, match = PETSC_FALSE; 56 PetscBool chol = PETSC_FALSE, view = PETSC_FALSE, matsolvexx = PETSC_FALSE; 57 #if defined(PETSC_HAVE_MUMPS) 58 PetscBool test_mumps_opts = PETSC_FALSE; 59 #endif 60 PetscViewer fd; /* viewer */ 61 char file[PETSC_MAX_PATH_LEN]; /* input file name */ 62 char pack[PETSC_MAX_PATH_LEN]; 63 64 PetscFunctionBeginUser; 65 PetscCall(PetscInitialize(&argc, &args, NULL, help)); 66 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 67 68 /* Determine file from which we read the matrix A */ 69 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg)); 70 if (flg) { /* Load matrix A */ 71 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd)); 72 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 73 PetscCall(MatSetFromOptions(A)); 74 PetscCall(MatLoad(A, fd)); 75 PetscCall(PetscViewerDestroy(&fd)); 76 } else { 77 n = 13; 78 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 79 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 80 PetscCall(MatSetType(A, MATAIJ)); 81 PetscCall(MatSetFromOptions(A)); 82 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n)); 83 PetscCall(MatSetUp(A)); 84 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 85 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 86 PetscCall(MatShift(A, 1.0)); 87 } 88 89 /* if A is symmetric, set its flag -- required by MatGetInertia() */ 90 PetscCall(MatIsSymmetric(A, 0.0, &symm)); 91 PetscCall(MatSetOption(A, MAT_SYMMETRIC, symm)); 92 93 PetscCall(PetscOptionsGetBool(NULL, NULL, "-cholesky", &chol, NULL)); 94 95 /* test MATNEST support */ 96 flg = PETSC_FALSE; 97 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nest", &flg, NULL)); 98 if (flg) { 99 Mat B; 100 101 flg = PETSC_FALSE; 102 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nest_bordered", &flg, NULL)); 103 if (!flg) { 104 Mat mats[9] = {NULL, NULL, A, NULL, A, NULL, A, NULL, NULL}; 105 106 /* Create a nested matrix representing 107 | 0 0 A | 108 | 0 A 0 | 109 | A 0 0 | 110 */ 111 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 3, NULL, 3, NULL, mats, &B)); 112 flg = PETSC_TRUE; 113 } else { 114 Mat mats[4]; 115 116 /* Create a nested matrix representing 117 | Id R | 118 | R^t A | 119 */ 120 PetscCall(MatGetSize(A, NULL, &n)); 121 m = n + 12; 122 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 123 PetscCall(CreateIdentity(m, &mats[0])); 124 PetscCall(CreateRandom(m, n, &mats[1])); 125 mats[3] = A; 126 127 /* use CreateTranspose/CreateHermitianTranspose or explicit matrix for debugging purposes */ 128 flg = PETSC_FALSE; 129 PetscCall(PetscOptionsGetBool(NULL, NULL, "-expl", &flg, NULL)); 130 #if PetscDefined(USE_COMPLEX) 131 if (chol) { /* Hermitian transpose not supported by MUMPS Cholesky factor */ 132 if (!flg) PetscCall(MatCreateTranspose(mats[1], &mats[2])); 133 else PetscCall(MatTranspose(mats[1], MAT_INITIAL_MATRIX, &mats[2])); 134 flg = PETSC_TRUE; 135 } else { 136 if (!flg) { 137 Mat B; 138 139 PetscCall(MatDuplicate(mats[1], MAT_COPY_VALUES, &B)); 140 PetscCall(MatCreateHermitianTranspose(B, &mats[2])); 141 PetscCall(MatDestroy(&B)); 142 if (n == m) { 143 PetscCall(MatScale(mats[2], PetscCMPLX(4.0, -2.0))); 144 PetscCall(MatShift(mats[2], PetscCMPLX(-2.0, 1.0))); // mats[2] = (4 - 2i) B* - (2 - i) I 145 PetscCall(MatCreateHermitianTranspose(mats[2], &B)); 146 PetscCall(MatDestroy(mats + 2)); 147 PetscCall(MatScale(B, 0.5)); 148 PetscCall(MatShift(B, PetscCMPLX(1.0, 0.5))); // B = 0.5 mats[2]* - (1 - 0.5i) I = (2 + i) B - (1 + 0.5i) I + (1 + 0.5i) I = (2 + i) B 149 PetscCall(MatCreateHermitianTranspose(B, &mats[2])); // mats[2] = B* = (2 - i) B* 150 PetscCall(MatDestroy(&B)); 151 PetscCall(MatScale(mats[1], PetscCMPLX(2.0, 1.0))); // mats[1] = (2 + i) B = mats[2]* 152 } else flg = PETSC_TRUE; 153 } else PetscCall(MatHermitianTranspose(mats[1], MAT_INITIAL_MATRIX, &mats[2])); 154 } 155 #else 156 if (!flg) { 157 Mat B; 158 159 PetscCall(MatDuplicate(mats[1], MAT_COPY_VALUES, &B)); 160 PetscCall(MatCreateTranspose(B, &mats[2])); 161 PetscCall(MatDestroy(&B)); 162 if (n == m) { 163 PetscCall(MatScale(mats[2], 4.0)); 164 PetscCall(MatShift(mats[2], -2.0)); // mats[2] = 4 B' - 2 I 165 PetscCall(MatCreateTranspose(mats[2], &B)); 166 PetscCall(MatDestroy(mats + 2)); 167 PetscCall(MatScale(B, 0.5)); 168 PetscCall(MatShift(B, 1.0)); // B = 0.5 mats[2]' + I = 0.5 (4 B' - 2 I)' + I = 2 B 169 PetscCall(MatCreateTranspose(B, &mats[2])); // mats[2] = B' = 2 B' 170 PetscCall(MatDestroy(&B)); 171 PetscCall(MatScale(mats[1], 2.0)); // mats[1] = 2 B = mats[2]' 172 } else flg = PETSC_TRUE; 173 } else PetscCall(MatTranspose(mats[1], MAT_INITIAL_MATRIX, &mats[2])); 174 #endif 175 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, mats, &B)); 176 PetscCall(MatDestroy(&mats[0])); 177 PetscCall(MatDestroy(&mats[1])); 178 PetscCall(MatDestroy(&mats[2])); 179 } 180 PetscCall(MatDestroy(&A)); 181 A = B; 182 PetscCall(MatSetOption(A, MAT_SYMMETRIC, symm)); 183 184 /* not all the combinations of MatMat operations are supported by MATNEST. */ 185 PetscCall(MatComputeOperator(A, MATAIJ, &Ae)); 186 } else { 187 PetscCall(PetscObjectReference((PetscObject)A)); 188 Ae = A; 189 flg = PETSC_TRUE; 190 } 191 PetscCall(MatGetLocalSize(A, &m, &n)); 192 PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n); 193 194 PetscCall(MatViewFromOptions(A, NULL, "-A_view")); 195 PetscCall(MatViewFromOptions(Ae, NULL, "-A_view_expl")); 196 197 /* Create dense matrix C and X; C holds true solution with identical columns */ 198 nrhs = 2; 199 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL)); 200 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ex125: nrhs %" PetscInt_FMT "\n", nrhs)); 201 PetscCall(MatCreate(PETSC_COMM_WORLD, &C)); 202 PetscCall(MatSetOptionsPrefix(C, "rhs_")); 203 PetscCall(MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs)); 204 PetscCall(MatSetType(C, MATDENSE)); 205 PetscCall(MatSetFromOptions(C)); 206 PetscCall(MatSetUp(C)); 207 208 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_factor", &view, NULL)); 209 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmatsolve", &testMatMatSolve, NULL)); 210 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmatsolvetranspose", &testMatMatSolveTranspose, NULL)); 211 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matsolvetranspose", &testMatSolveTranspose, NULL)); 212 #if defined(PETSC_HAVE_MUMPS) 213 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_mumps_opts", &test_mumps_opts, NULL)); 214 #endif 215 216 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand)); 217 PetscCall(PetscRandomSetFromOptions(rand)); 218 PetscCall(MatSetRandom(C, rand)); 219 PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &X)); 220 221 /* Create vectors */ 222 PetscCall(MatCreateVecs(A, &x, &b)); 223 PetscCall(VecDuplicate(x, &u)); /* save the true solution */ 224 225 /* Test Factorization */ 226 if (flg) PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm)); // TODO FIXME: MatConvert_Nest_AIJ() does not support chained MatCreate[Hermitian]Transpose() 227 228 PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_solver_type", pack, sizeof(pack), NULL)); 229 #if defined(PETSC_HAVE_SUPERLU) 230 PetscCall(PetscStrcmp(MATSOLVERSUPERLU, pack, &match)); 231 if (match) { 232 PetscCheck(!chol, PETSC_COMM_WORLD, PETSC_ERR_SUP, "SuperLU does not provide Cholesky!"); 233 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SUPERLU LU:\n")); 234 PetscCall(MatGetFactor(A, MATSOLVERSUPERLU, MAT_FACTOR_LU, &F)); 235 matsolvexx = PETSC_FALSE; /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix, need further work */ 236 ipack = 0; 237 goto skipoptions; 238 } 239 #endif 240 #if defined(PETSC_HAVE_SUPERLU_DIST) 241 PetscCall(PetscStrcmp(MATSOLVERSUPERLU_DIST, pack, &match)); 242 if (match) { 243 PetscCheck(!chol, PETSC_COMM_WORLD, PETSC_ERR_SUP, "SuperLU does not provide Cholesky!"); 244 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SUPERLU_DIST LU:\n")); 245 PetscCall(MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_LU, &F)); 246 matsolvexx = PETSC_TRUE; 247 if (symm) { /* A is symmetric */ 248 testMatMatSolveTranspose = PETSC_TRUE; 249 testMatSolveTranspose = PETSC_TRUE; 250 } else { /* superlu_dist does not support solving A^t x = rhs yet */ 251 testMatMatSolveTranspose = PETSC_FALSE; 252 testMatSolveTranspose = PETSC_FALSE; 253 } 254 ipack = 1; 255 goto skipoptions; 256 } 257 #endif 258 #if defined(PETSC_HAVE_MUMPS) 259 PetscCall(PetscStrcmp(MATSOLVERMUMPS, pack, &match)); 260 if (match) { 261 if (chol) { 262 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MUMPS CHOLESKY:\n")); 263 PetscCall(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_CHOLESKY, &F)); 264 } else { 265 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MUMPS LU:\n")); 266 PetscCall(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_LU, &F)); 267 } 268 matsolvexx = PETSC_TRUE; 269 if (test_mumps_opts) { 270 /* test mumps options */ 271 PetscInt icntl; 272 PetscReal cntl; 273 274 icntl = 2; /* sequential matrix ordering */ 275 PetscCall(MatMumpsSetIcntl(F, 7, icntl)); 276 277 cntl = 1.e-6; /* threshold for row pivot detection */ 278 PetscCall(MatMumpsSetIcntl(F, 24, 1)); 279 PetscCall(MatMumpsSetCntl(F, 3, cntl)); 280 } 281 ipack = 2; 282 goto skipoptions; 283 } 284 #endif 285 #if defined(PETSC_HAVE_MKL_PARDISO) 286 PetscCall(PetscStrcmp(MATSOLVERMKL_PARDISO, pack, &match)); 287 if (match) { 288 if (chol) { 289 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO CHOLESKY:\n")); 290 PetscCall(MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_CHOLESKY, &F)); 291 } else { 292 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO LU:\n")); 293 PetscCall(MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_LU, &F)); 294 } 295 ipack = 3; 296 goto skipoptions; 297 } 298 #endif 299 #if defined(PETSC_HAVE_CUDA) 300 PetscCall(PetscStrcmp(MATSOLVERCUSPARSE, pack, &match)); 301 if (match) { 302 if (chol) { 303 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE CHOLESKY:\n")); 304 PetscCall(MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_CHOLESKY, &F)); 305 } else { 306 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE LU:\n")); 307 PetscCall(MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_LU, &F)); 308 } 309 testMatSolveTranspose = PETSC_FALSE; 310 testMatMatSolveTranspose = PETSC_FALSE; 311 ipack = 4; 312 goto skipoptions; 313 } 314 #endif 315 /* PETSc */ 316 match = PETSC_TRUE; 317 if (match) { 318 if (chol) { 319 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PETSC CHOLESKY:\n")); 320 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &F)); 321 } else { 322 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PETSC LU:\n")); 323 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F)); 324 } 325 matsolvexx = PETSC_TRUE; 326 ipack = 5; 327 goto skipoptions; 328 } 329 330 skipoptions: 331 PetscCall(MatFactorInfoInitialize(&info)); 332 info.fill = 5.0; 333 info.shifttype = (PetscReal)MAT_SHIFT_NONE; 334 if (chol) { 335 PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info)); 336 } else { 337 PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info)); 338 } 339 340 for (nfact = 0; nfact < 2; nfact++) { 341 if (chol) { 342 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the CHOLESKY numfactorization \n", nfact)); 343 PetscCall(MatCholeskyFactorNumeric(F, A, &info)); 344 } else { 345 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the LU numfactorization \n", nfact)); 346 PetscCall(MatLUFactorNumeric(F, A, &info)); 347 } 348 if (view) { 349 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO)); 350 PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD)); 351 PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD)); 352 view = PETSC_FALSE; 353 } 354 355 #if defined(PETSC_HAVE_SUPERLU_DIST) 356 if (ipack == 1) { /* Test MatSuperluDistGetDiagU() 357 -- input: matrix factor F; output: main diagonal of matrix U on all processes */ 358 PetscInt M; 359 PetscScalar *diag; 360 #if !defined(PETSC_USE_COMPLEX) 361 PetscInt nneg, nzero, npos; 362 #endif 363 364 PetscCall(MatGetSize(F, &M, NULL)); 365 PetscCall(PetscMalloc1(M, &diag)); 366 PetscCall(MatSuperluDistGetDiagU(F, diag)); 367 PetscCall(PetscFree(diag)); 368 369 #if !defined(PETSC_USE_COMPLEX) 370 /* Test MatGetInertia() */ 371 if (symm) { /* A is symmetric */ 372 PetscCall(MatGetInertia(F, &nneg, &nzero, &npos)); 373 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos)); 374 } 375 #endif 376 } 377 #endif 378 379 #if defined(PETSC_HAVE_MUMPS) 380 /* mumps interface allows repeated call of MatCholeskyFactorSymbolic(), while the succession calls do nothing */ 381 if (ipack == 2) { 382 if (chol) { 383 PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info)); 384 PetscCall(MatCholeskyFactorNumeric(F, A, &info)); 385 } else { 386 PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info)); 387 PetscCall(MatLUFactorNumeric(F, A, &info)); 388 } 389 } 390 #endif 391 392 /* Test MatMatSolve(), A X = B, where B can be dense or sparse */ 393 if (testMatMatSolve) { 394 if (!nfact) { 395 PetscCall(MatMatMult(Ae, C, MAT_INITIAL_MATRIX, 2.0, &RHS)); 396 } else { 397 PetscCall(MatMatMult(Ae, C, MAT_REUSE_MATRIX, 2.0, &RHS)); 398 } 399 for (nsolve = 0; nsolve < 2; nsolve++) { 400 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatMatSolve \n", nsolve)); 401 PetscCall(MatMatSolve(F, RHS, X)); 402 403 /* Check the error */ 404 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN)); 405 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm)); 406 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve)); 407 } 408 409 if (matsolvexx) { 410 /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */ 411 PetscCall(MatCopy(RHS, X, SAME_NONZERO_PATTERN)); 412 PetscCall(MatMatSolve(F, X, X)); 413 /* Check the error */ 414 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN)); 415 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm)); 416 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMatSolve(F,RHS,RHS): Norm of error %g\n", (double)norm)); 417 } 418 419 if (ipack == 2 && size == 1) { 420 Mat spRHS, spRHST, RHST; 421 422 PetscCall(MatTranspose(RHS, MAT_INITIAL_MATRIX, &RHST)); 423 PetscCall(MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST)); 424 PetscCall(MatCreateTranspose(spRHST, &spRHS)); 425 for (nsolve = 0; nsolve < 2; nsolve++) { 426 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the sparse MatMatSolve \n", nsolve)); 427 PetscCall(MatMatSolve(F, spRHS, X)); 428 429 /* Check the error */ 430 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN)); 431 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm)); 432 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)); 433 } 434 PetscCall(MatDestroy(&spRHST)); 435 PetscCall(MatDestroy(&spRHS)); 436 PetscCall(MatDestroy(&RHST)); 437 } 438 } 439 440 /* Test testMatMatSolveTranspose(), A^T X = B, where B can be dense or sparse */ 441 if (testMatMatSolveTranspose) { 442 if (!nfact) { 443 PetscCall(MatTransposeMatMult(Ae, C, MAT_INITIAL_MATRIX, 2.0, &RHS1)); 444 } else { 445 PetscCall(MatTransposeMatMult(Ae, C, MAT_REUSE_MATRIX, 2.0, &RHS1)); 446 } 447 448 for (nsolve = 0; nsolve < 2; nsolve++) { 449 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatMatSolveTranspose\n", nsolve)); 450 PetscCall(MatMatSolveTranspose(F, RHS1, X)); 451 452 /* Check the error */ 453 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN)); 454 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm)); 455 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve)); 456 } 457 458 if (ipack == 2 && size == 1) { 459 Mat spRHS, spRHST, RHST; 460 461 PetscCall(MatTranspose(RHS1, MAT_INITIAL_MATRIX, &RHST)); 462 PetscCall(MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST)); 463 PetscCall(MatCreateTranspose(spRHST, &spRHS)); 464 for (nsolve = 0; nsolve < 2; nsolve++) { 465 PetscCall(MatMatSolveTranspose(F, spRHS, X)); 466 467 /* Check the error */ 468 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN)); 469 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm)); 470 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)); 471 } 472 PetscCall(MatDestroy(&spRHST)); 473 PetscCall(MatDestroy(&spRHS)); 474 PetscCall(MatDestroy(&RHST)); 475 } 476 } 477 478 /* Test MatSolve() */ 479 if (testMatSolve) { 480 for (nsolve = 0; nsolve < 2; nsolve++) { 481 PetscCall(VecSetRandom(x, rand)); 482 PetscCall(VecCopy(x, u)); 483 PetscCall(MatMult(Ae, x, b)); 484 485 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatSolve \n", nsolve)); 486 PetscCall(MatSolve(F, b, x)); 487 488 /* Check the error */ 489 PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */ 490 PetscCall(VecNorm(u, NORM_2, &norm)); 491 if (norm > tol) { 492 PetscReal resi; 493 PetscCall(MatMult(Ae, x, u)); /* u = A*x */ 494 PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */ 495 PetscCall(VecNorm(u, NORM_2, &resi)); 496 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolve: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact)); 497 } 498 } 499 } 500 501 /* Test MatSolveTranspose() */ 502 if (testMatSolveTranspose) { 503 for (nsolve = 0; nsolve < 2; nsolve++) { 504 PetscCall(VecSetRandom(x, rand)); 505 PetscCall(VecCopy(x, u)); 506 PetscCall(MatMultTranspose(Ae, x, b)); 507 508 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatSolveTranspose\n", nsolve)); 509 PetscCall(MatSolveTranspose(F, b, x)); 510 511 /* Check the error */ 512 PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */ 513 PetscCall(VecNorm(u, NORM_2, &norm)); 514 if (norm > tol) { 515 PetscReal resi; 516 PetscCall(MatMultTranspose(Ae, x, u)); /* u = A*x */ 517 PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */ 518 PetscCall(VecNorm(u, NORM_2, &resi)); 519 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveTranspose: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact)); 520 } 521 } 522 } 523 } 524 525 /* Free data structures */ 526 PetscCall(MatDestroy(&Ae)); 527 PetscCall(MatDestroy(&A)); 528 PetscCall(MatDestroy(&C)); 529 PetscCall(MatDestroy(&F)); 530 PetscCall(MatDestroy(&X)); 531 PetscCall(MatDestroy(&RHS)); 532 PetscCall(MatDestroy(&RHS1)); 533 534 PetscCall(PetscRandomDestroy(&rand)); 535 PetscCall(ISDestroy(&perm)); 536 PetscCall(ISDestroy(&iperm)); 537 PetscCall(VecDestroy(&x)); 538 PetscCall(VecDestroy(&b)); 539 PetscCall(VecDestroy(&u)); 540 PetscCall(PetscFinalize()); 541 return 0; 542 } 543 544 /*TEST 545 546 test: 547 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 548 args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type petsc 549 output_file: output/ex125.out 550 551 test: 552 suffix: 2 553 args: -mat_solver_type petsc 554 output_file: output/ex125.out 555 556 test: 557 suffix: mkl_pardiso 558 requires: mkl_pardiso datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 559 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mkl_pardiso 560 561 test: 562 suffix: mkl_pardiso_2 563 requires: mkl_pardiso 564 args: -mat_solver_type mkl_pardiso 565 output_file: output/ex125_mkl_pardiso.out 566 567 test: 568 suffix: mumps 569 requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 570 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps 571 output_file: output/ex125_mumps_seq.out 572 573 test: 574 suffix: mumps_nest 575 requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 576 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}} 577 output_file: output/ex125_mumps_seq.out 578 579 test: 580 suffix: mumps_2 581 nsize: 3 582 requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 583 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps 584 output_file: output/ex125_mumps_par.out 585 586 test: 587 suffix: mumps_2_nest 588 nsize: 3 589 requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 590 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}} 591 output_file: output/ex125_mumps_par.out 592 593 test: 594 suffix: mumps_3 595 requires: mumps 596 args: -mat_solver_type mumps 597 output_file: output/ex125_mumps_seq.out 598 599 test: 600 suffix: mumps_3_nest 601 requires: mumps 602 args: -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}} 603 output_file: output/ex125_mumps_seq.out 604 605 test: 606 suffix: mumps_4 607 nsize: 3 608 requires: mumps 609 args: -mat_solver_type mumps 610 output_file: output/ex125_mumps_par.out 611 612 test: 613 suffix: mumps_4_nest 614 nsize: 3 615 requires: mumps 616 args: -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}} 617 output_file: output/ex125_mumps_par.out 618 619 test: 620 suffix: mumps_5 621 nsize: 3 622 requires: mumps 623 args: -mat_solver_type mumps -cholesky 624 output_file: output/ex125_mumps_par_cholesky.out 625 626 test: 627 suffix: mumps_5_nest 628 nsize: 3 629 requires: mumps 630 args: -mat_solver_type mumps -cholesky -test_nest -test_nest_bordered {{0 1}} 631 output_file: output/ex125_mumps_par_cholesky.out 632 633 test: 634 suffix: mumps_6 635 nsize: 2 636 requires: mumps 637 args: -mat_solver_type mumps -test_nest -test_nest_bordered -m 13 -n 13 638 output_file: output/ex125_mumps_par.out 639 640 test: 641 suffix: superlu 642 requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu 643 args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type superlu 644 output_file: output/ex125_superlu.out 645 646 test: 647 suffix: superlu_dist 648 nsize: {{1 3}} 649 requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist 650 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM 651 output_file: output/ex125_superlu_dist.out 652 653 test: 654 suffix: superlu_dist_2 655 nsize: {{1 3}} 656 requires: superlu_dist !complex 657 args: -n 36 -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM 658 output_file: output/ex125_superlu_dist.out 659 660 test: 661 suffix: superlu_dist_3 662 nsize: {{1 3}} 663 requires: superlu_dist !complex 664 requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist 665 args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM 666 output_file: output/ex125_superlu_dist_nonsymmetric.out 667 668 test: 669 suffix: superlu_dist_complex 670 nsize: 3 671 requires: datafilespath double superlu_dist complex !defined(PETSC_USE_64BIT_INDICES) 672 args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type superlu_dist 673 output_file: output/ex125_superlu_dist_complex.out 674 675 test: 676 suffix: superlu_dist_complex_2 677 nsize: 3 678 requires: superlu_dist complex 679 args: -mat_solver_type superlu_dist 680 output_file: output/ex125_superlu_dist_complex_2.out 681 682 test: 683 suffix: cusparse 684 requires: cuda datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 685 #TODO: fix the bug with cholesky 686 #args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type cusparse -cholesky {{0 1}separate output} 687 args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type cusparse -cholesky {{0}separate output} 688 689 test: 690 suffix: cusparse_2 691 requires: cuda 692 args: -mat_type aijcusparse -mat_solver_type cusparse -cholesky {{0 1}separate output} 693 694 TEST*/ 695