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;
 12:   PetscMPIInt    rank,size;
 13:   PetscInt       m,n,nfact,ipack=0;
 14:   PetscReal      norm,tol=1.e-12,Anorm;
 15:   IS             perm,iperm;
 16:   MatFactorInfo  info;
 17:   PetscBool      flg,testMatSolve=PETSC_TRUE;
 18:   PetscViewer    fd;              /* viewer */
 19:   char           file[PETSC_MAX_PATH_LEN]; /* input file name */

 21:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 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);
 27:   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option");

 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);
 37:   if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", 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) 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) 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) 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) 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) {
 88:           PetscPrintf(PETSC_COMM_SELF,"MatSolve: rel residual %g/%g = %g, LU numfact %d\n",norm,Anorm,norm/Anorm,nfact);
 89:         }
 90:       }
 91:     }
 92:   }

 94:   /* Free data structures */
 95:   MatDestroy(&A);
 96:   MatDestroy(&F);
 97:   ISDestroy(&perm);
 98:   ISDestroy(&iperm);
 99:   VecDestroy(&x);
100:   VecDestroy(&b);
101:   VecDestroy(&u);
102:   PetscFinalize();
103:   return ierr;
104: }