Actual source code: ex28.c
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 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: MatFactorType facttype = MAT_FACTOR_LU;
17: char solvertype[64];
18: char factortype[64];
20: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
21: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
23: /* Create and assemble matrices, all have same data structure */
24: for (k=0; k<num_numfac; k++) {
25: MatCreate(PETSC_COMM_WORLD,&A[k]);
26: MatSetSizes(A[k],PETSC_DECIDE,PETSC_DECIDE,N,N);
27: MatSetFromOptions(A[k]);
28: MatSetUp(A[k]);
29: MatGetOwnershipRange(A[k],&rstart,&rend);
31: value[0] = -1.0*(k+1);
32: value[1] = 2.0*(k+1);
33: value[2] = -1.0*(k+1);
34: for (i=rstart; i<rend; i++) {
35: col[0] = i-1; col[1] = i; col[2] = i+1;
36: if (i == 0) {
37: MatSetValues(A[k],1,&i,2,col+1,value+1,INSERT_VALUES);
38: } else if (i == N-1) {
39: MatSetValues(A[k],1,&i,2,col,value,INSERT_VALUES);
40: } else {
41: MatSetValues(A[k],1,&i,3,col,value,INSERT_VALUES);
42: }
43: }
44: MatAssemblyBegin(A[k],MAT_FINAL_ASSEMBLY);
45: MatAssemblyEnd(A[k],MAT_FINAL_ASSEMBLY);
46: MatSetOption(A[k],MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
47: }
49: /* Create vectors */
50: MatCreateVecs(A[0],&x,&b);
51: VecDuplicate(x,&u);
53: /* Set rhs vector b */
54: VecSet(b,1.0);
56: /* Get a symbolic factor F from A[0] */
57: PetscStrncpy(solvertype,"petsc",sizeof(solvertype));
58: PetscOptionsGetString(NULL, NULL, "-mat_solver_type",solvertype,sizeof(solvertype),NULL);
59: PetscOptionsGetEnum(NULL,NULL,"-mat_factor_type",MatFactorTypes,(PetscEnum*)&facttype,NULL);
61: MatGetFactor(A[0],solvertype,facttype,&F);
62: /* test mumps options */
63: #if defined(PETSC_HAVE_MUMPS)
64: MatMumpsSetIcntl(F,7,5);
65: #endif
66: PetscStrncpy(factortype,MatFactorTypes[facttype],sizeof(factortype));
67: PetscStrtoupper(solvertype);
68: PetscStrtoupper(factortype);
69: PetscPrintf(PETSC_COMM_WORLD," %s %s:\n",solvertype,factortype);
71: MatFactorInfoInitialize(&info);
72: info.fill = 5.0;
73: MatGetOrdering(A[0],MATORDERINGNATURAL,&perm,&iperm);
74: switch (facttype) {
75: case MAT_FACTOR_LU:
76: MatLUFactorSymbolic(F,A[0],perm,iperm,&info);
77: break;
78: case MAT_FACTOR_ILU:
79: MatILUFactorSymbolic(F,A[0],perm,iperm,&info);
80: break;
81: case MAT_FACTOR_ICC:
82: MatICCFactorSymbolic(F,A[0],perm,&info);
83: break;
84: case MAT_FACTOR_CHOLESKY:
85: MatCholeskyFactorSymbolic(F,A[0],perm,&info);
86: break;
87: default:
88: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not for factor type %s\n",factortype);
89: }
91: /* Compute numeric factors using same F, then solve */
92: for (k = 0; k < num_numfac; k++) {
93: switch (facttype) {
94: case MAT_FACTOR_LU:
95: case MAT_FACTOR_ILU:
96: MatLUFactorNumeric(F,A[k],&info);
97: break;
98: case MAT_FACTOR_ICC:
99: case MAT_FACTOR_CHOLESKY:
100: MatCholeskyFactorNumeric(F,A[k],&info);
101: break;
102: default:
103: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not for factor type %s\n",factortype);
104: }
106: /* Solve A[k] * x = b */
107: MatSolve(F,b,x);
109: /* Check the residual */
110: MatMult(A[k],x,u);
111: VecAXPY(u,-1.0,b);
112: VecNorm(u,NORM_INFINITY,&norm);
113: if (norm > tol) {
114: PetscPrintf(PETSC_COMM_WORLD,"%D-the %s numfact and solve: residual %g\n",k,factortype,(double)norm);
115: }
116: }
118: /* Free data structures */
119: for (k=0; k<num_numfac; k++) {
120: MatDestroy(&A[k]);
121: }
122: MatDestroy(&F);
123: ISDestroy(&perm);
124: ISDestroy(&iperm);
125: VecDestroy(&x);
126: VecDestroy(&b);
127: VecDestroy(&u);
128: PetscFinalize();
129: return ierr;
130: }
132: /*TEST
134: test:
136: test:
137: suffix: 2
138: args: -mat_solver_type superlu
139: requires: superlu
141: test:
142: suffix: 3
143: nsize: 2
144: requires: mumps
145: args: -mat_solver_type mumps
147: test:
148: suffix: 4
149: args: -mat_solver_type cusparse -mat_type aijcusparse -mat_factor_type {{lu cholesky ilu icc}separate output}
150: requires: cuda
152: TEST*/