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";
3
4 #include <petscmat.h>
5
main(int argc,char ** args)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 */
18
19 PetscFunctionBeginUser;
20 PetscCall(PetscInitialize(&argc, &args, NULL, help));
21 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
22 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
23
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");
27
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));
38
39 /* Create vectors */
40 PetscCall(VecDuplicate(b, &x));
41 PetscCall(VecDuplicate(x, &u)); /* save the true solution */
42
43 /* Test LU Factorization */
44 PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &perm, &iperm));
45
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 }
69
70 info.fill = 5.0;
71 PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info));
72
73 for (nfact = 0; nfact < 1; nfact++) {
74 if (rank == 0) printf(" %d-the LU numfactorization \n", nfact);
75 PetscCall(MatLUFactorNumeric(F, A, &info));
76
77 /* Test MatSolve() */
78 if (testMatSolve) {
79 PetscCall(MatSolve(F, b, x));
80
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 }
90
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 }
102