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, NULL, 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: }