Actual source code: ex130.c
1: static char help[] = "Tests external direct solvers. Simplified from ex125.c\n\
2: Example: mpiexec -n <np> ./ex130 -f <matrix binary file> -mat_solver_type 1 -mat_superlu_equil \n\n";
4: #include <petscmat.h>
6: int main(int argc, char **args)
7: {
8: Mat A, F;
9: Vec u, x, b;
10: PetscMPIInt rank, size;
11: PetscInt m, n, nfact, ipack = 0;
12: PetscReal norm, tol = 1.e-12, Anorm;
13: IS perm, iperm;
14: MatFactorInfo info;
15: PetscBool flg, testMatSolve = PETSC_TRUE;
16: PetscViewer fd; /* viewer */
17: char file[PETSC_MAX_PATH_LEN]; /* input file name */
19: PetscFunctionBeginUser;
20: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
21: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
22: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
24: /* Determine file from which we read the matrix A */
25: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
26: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
28: /* Load matrix A */
29: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
30: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
31: PetscCall(MatLoad(A, fd));
32: PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
33: PetscCall(VecLoad(b, fd));
34: PetscCall(PetscViewerDestroy(&fd));
35: PetscCall(MatGetLocalSize(A, &m, &n));
36: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
37: PetscCall(MatNorm(A, NORM_INFINITY, &Anorm));
39: /* Create vectors */
40: PetscCall(VecDuplicate(b, &x));
41: PetscCall(VecDuplicate(x, &u)); /* save the true solution */
43: /* Test LU Factorization */
44: PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &perm, &iperm));
46: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_solver_type", &ipack, NULL));
47: switch (ipack) {
48: case 1:
49: #if defined(PETSC_HAVE_SUPERLU)
50: if (rank == 0) printf(" SUPERLU LU:\n");
51: PetscCall(MatGetFactor(A, MATSOLVERSUPERLU, MAT_FACTOR_LU, &F));
52: break;
53: #endif
54: case 2:
55: #if defined(PETSC_HAVE_MUMPS)
56: if (rank == 0) printf(" MUMPS LU:\n");
57: PetscCall(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_LU, &F));
58: {
59: /* test mumps options */
60: PetscInt icntl_7 = 5;
61: PetscCall(MatMumpsSetIcntl(F, 7, icntl_7));
62: }
63: break;
64: #endif
65: default:
66: if (rank == 0) printf(" PETSC LU:\n");
67: PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F));
68: }
70: info.fill = 5.0;
71: PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info));
73: for (nfact = 0; nfact < 1; nfact++) {
74: if (rank == 0) printf(" %d-the LU numfactorization \n", nfact);
75: PetscCall(MatLUFactorNumeric(F, A, &info));
77: /* Test MatSolve() */
78: if (testMatSolve) {
79: PetscCall(MatSolve(F, b, x));
81: /* Check the residual */
82: PetscCall(MatMult(A, x, u));
83: PetscCall(VecAXPY(u, -1.0, b));
84: PetscCall(VecNorm(u, NORM_INFINITY, &norm));
85: if (norm > tol) {
86: if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSolve: rel residual %g/%g = %g, LU numfact %d\n", norm, Anorm, norm / Anorm, nfact));
87: }
88: }
89: }
91: /* Free data structures */
92: PetscCall(MatDestroy(&A));
93: PetscCall(MatDestroy(&F));
94: PetscCall(ISDestroy(&perm));
95: PetscCall(ISDestroy(&iperm));
96: PetscCall(VecDestroy(&x));
97: PetscCall(VecDestroy(&b));
98: PetscCall(VecDestroy(&u));
99: PetscCall(PetscFinalize());
100: return 0;
101: }