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