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