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 */
20: PetscInitialize(&argc,&args,(char*)0,help);
21: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
22: MPI_Comm_size(PETSC_COMM_WORLD, &size);
24: /* Determine file from which we read the matrix A */
25: PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);
28: /* Load matrix A */
29: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
30: MatCreate(PETSC_COMM_WORLD,&A);
31: MatLoad(A,fd);
32: VecCreate(PETSC_COMM_WORLD,&b);
33: VecLoad(b,fd);
34: PetscViewerDestroy(&fd);
35: MatGetLocalSize(A,&m,&n);
37: MatNorm(A,NORM_INFINITY,&Anorm);
39: /* Create vectors */
40: VecDuplicate(b,&x);
41: VecDuplicate(x,&u); /* save the true solution */
43: /* Test LU Factorization */
44: MatGetOrdering(A,MATORDERINGNATURAL,&perm,&iperm);
46: 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: 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: MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);
58: {
59: /* test mumps options */
60: PetscInt icntl_7 = 5;
61: MatMumpsSetIcntl(F,7,icntl_7);
62: }
63: break;
64: #endif
65: default:
66: if (rank == 0) printf(" PETSC LU:\n");
67: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
68: }
70: info.fill = 5.0;
71: 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: MatLUFactorNumeric(F,A,&info);
77: /* Test MatSolve() */
78: if (testMatSolve) {
79: MatSolve(F,b,x);
81: /* Check the residual */
82: MatMult(A,x,u);
83: VecAXPY(u,-1.0,b);
84: VecNorm(u,NORM_INFINITY,&norm);
85: if (norm > tol) {
86: if (rank == 0) {
87: PetscPrintf(PETSC_COMM_SELF,"MatSolve: rel residual %g/%g = %g, LU numfact %d\n",norm,Anorm,norm/Anorm,nfact);
88: }
89: }
90: }
91: }
93: /* Free data structures */
94: MatDestroy(&A);
95: MatDestroy(&F);
96: ISDestroy(&perm);
97: ISDestroy(&iperm);
98: VecDestroy(&x);
99: VecDestroy(&b);
100: VecDestroy(&u);
101: PetscFinalize();
102: return 0;
103: }