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