Actual source code: ex28.c
petsc-3.12.0 2019-09-29
1: static char help[] = "Illustrate how to do one symbolic factorization and multiple numeric factorizations using same matrix structure. \n\n";
3: #include <petscmat.h>
5: int main(int argc,char **args)
6: {
7: PetscInt ipack,i,rstart,rend,N=10,num_numfac=5,col[3],k;
8: Mat A[5],F;
9: Vec u,x,b;
11: PetscMPIInt rank;
12: PetscScalar value[3];
13: PetscReal norm,tol=100*PETSC_MACHINE_EPSILON;
14: IS perm,iperm;
15: MatFactorInfo info;
16: char solvertype[64]="petsc";
17: PetscBool flg,flg_superlu,flg_mumps;
19: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
20: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
22: /* Create and assemble matrices, all have same data structure */
23: for (k=0; k<num_numfac; k++) {
24: MatCreate(PETSC_COMM_WORLD,&A[k]);
25: MatSetSizes(A[k],PETSC_DECIDE,PETSC_DECIDE,N,N);
26: MatSetFromOptions(A[k]);
27: MatSetUp(A[k]);
28: MatGetOwnershipRange(A[k],&rstart,&rend);
30: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
31: for (i=rstart; i<rend; i++) {
32: col[0] = i-1; col[1] = i; col[2] = i+1;
33: if (i == 0) {
34: MatSetValues(A[k],1,&i,2,col+1,value+1,INSERT_VALUES);
35: } else if (i == N-1) {
36: MatSetValues(A[k],1,&i,2,col,value,INSERT_VALUES);
37: } else {
38: MatSetValues(A[k],1,&i,3,col,value,INSERT_VALUES);
39: }
40: }
41: MatAssemblyBegin(A[k],MAT_FINAL_ASSEMBLY);
42: MatAssemblyEnd(A[k],MAT_FINAL_ASSEMBLY);
43: MatSetOption(A[k],MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
44: }
46: /* Create vectors */
47: MatCreateVecs(A[0],&x,&b);
48: VecDuplicate(x,&u);
50: /* Set rhs vector b */
51: VecSet(b,1.0);
53: /* Get a symbolic factor F from A[0] */
54: PetscOptionsGetString(NULL, NULL, "-mat_solver_type",solvertype,64,&flg);
55: PetscStrcmp(solvertype,"superlu",&flg_superlu);
56: PetscStrcmp(solvertype,"mumps",&flg_mumps);
58: ipack = 0;
59: if (flg_superlu) ipack = 1;
60: else {
61: if (flg_mumps) ipack = 2;
62: }
64: switch (ipack) {
65: case 1:
66: #if defined(PETSC_HAVE_SUPERLU)
67: PetscPrintf(PETSC_COMM_WORLD," SUPERLU LU:\n");
68: MatGetFactor(A[0],MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);
69: break;
70: #else
71: PetscPrintf(PETSC_COMM_WORLD," SUPERLU is not installed, using PETSC LU\n");
72: #endif
73: case 2:
74: #if defined(PETSC_HAVE_MUMPS)
75: PetscPrintf(PETSC_COMM_WORLD," MUMPS LU:\n");
76: MatGetFactor(A[0],MATSOLVERMUMPS,MAT_FACTOR_LU,&F);
77: {
78: /* test mumps options */
79: PetscInt icntl_7 = 5;
80: MatMumpsSetIcntl(F,7,icntl_7);
81: }
82: break;
83: #else
84: PetscPrintf(PETSC_COMM_WORLD," MUMPS is not installed, use PETSC LU\n");
85: #endif
86: default:
87: PetscPrintf(PETSC_COMM_WORLD," PETSC LU:\n");
88: MatGetFactor(A[0],MATSOLVERPETSC,MAT_FACTOR_LU,&F);
89: }
91: info.fill = 5.0;
92: MatGetOrdering(A[0],MATORDERINGNATURAL,&perm,&iperm);
93: MatLUFactorSymbolic(F,A[0],perm,iperm,&info);
95: /* Compute numeric factors using same F, then solve */
96: for (k = 0; k < num_numfac; k++) {
97: /* Get numeric factor of A[k] */
98: MatLUFactorNumeric(F,A[k],&info);
100: /* Solve A[k] * x = b */
101: MatSolve(F,b,x);
103: /* Check the residual */
104: MatMult(A[k],x,u);
105: VecAXPY(u,-1.0,b);
106: VecNorm(u,NORM_INFINITY,&norm);
107: if (norm > tol) {
108: PetscPrintf(PETSC_COMM_WORLD,"%D-the LU numfact and solve: residual %g\n",k,(double)norm);
109: }
110: }
112: /* Free data structures */
113: for (k=0; k<num_numfac; k++) {
114: MatDestroy(&A[k]);
115: }
116: MatDestroy(&F);
117: ISDestroy(&perm);
118: ISDestroy(&iperm);
119: VecDestroy(&x);
120: VecDestroy(&b);
121: VecDestroy(&u);
122: PetscFinalize();
123: return ierr;
124: }
126: /*TEST
128: test:
130: test:
131: suffix: 2
132: args: -mat_solver_type superlu
133: requires: superlu
135: test:
136: suffix: 3
137: nsize: 2
138: requires: mumps
139: args: -mat_solver_type mumps
141: TEST*/