1 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"; 4 5 #include <petscmat.h> 6 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 */ 19 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)); 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) { 87 PetscCall(PetscPrintf(PETSC_COMM_SELF,"MatSolve: rel residual %g/%g = %g, LU numfact %d\n",norm,Anorm,norm/Anorm,nfact)); 88 } 89 } 90 } 91 } 92 93 /* Free data structures */ 94 PetscCall(MatDestroy(&A)); 95 PetscCall(MatDestroy(&F)); 96 PetscCall(ISDestroy(&perm)); 97 PetscCall(ISDestroy(&iperm)); 98 PetscCall(VecDestroy(&x)); 99 PetscCall(VecDestroy(&b)); 100 PetscCall(VecDestroy(&u)); 101 PetscCall(PetscFinalize()); 102 return 0; 103 } 104