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;
10: PetscMPIInt rank;
11: PetscScalar value[3];
12: PetscReal norm,tol=100*PETSC_MACHINE_EPSILON;
13: IS perm,iperm;
14: MatFactorInfo info;
15: MatFactorType facttype = MAT_FACTOR_LU;
16: char solvertype[64];
17: char factortype[64];
19: PetscInitialize(&argc,&args,(char*)0,help);
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*(k+1);
31: value[1] = 2.0*(k+1);
32: value[2] = -1.0*(k+1);
33: for (i=rstart; i<rend; i++) {
34: col[0] = i-1; col[1] = i; col[2] = i+1;
35: if (i == 0) {
36: MatSetValues(A[k],1,&i,2,col+1,value+1,INSERT_VALUES);
37: } else if (i == N-1) {
38: MatSetValues(A[k],1,&i,2,col,value,INSERT_VALUES);
39: } else {
40: MatSetValues(A[k],1,&i,3,col,value,INSERT_VALUES);
41: }
42: }
43: MatAssemblyBegin(A[k],MAT_FINAL_ASSEMBLY);
44: MatAssemblyEnd(A[k],MAT_FINAL_ASSEMBLY);
45: MatSetOption(A[k],MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
46: }
48: /* Create vectors */
49: MatCreateVecs(A[0],&x,&b);
50: VecDuplicate(x,&u);
52: /* Set rhs vector b */
53: VecSet(b,1.0);
55: /* Get a symbolic factor F from A[0] */
56: PetscStrncpy(solvertype,"petsc",sizeof(solvertype));
57: PetscOptionsGetString(NULL, NULL, "-mat_solver_type",solvertype,sizeof(solvertype),NULL);
58: PetscOptionsGetEnum(NULL,NULL,"-mat_factor_type",MatFactorTypes,(PetscEnum*)&facttype,NULL);
60: MatGetFactor(A[0],solvertype,facttype,&F);
61: /* test mumps options */
62: #if defined(PETSC_HAVE_MUMPS)
63: MatMumpsSetIcntl(F,7,5);
64: #endif
65: PetscStrncpy(factortype,MatFactorTypes[facttype],sizeof(factortype));
66: PetscStrtoupper(solvertype);
67: PetscStrtoupper(factortype);
68: PetscPrintf(PETSC_COMM_WORLD," %s %s:\n",solvertype,factortype);
70: MatFactorInfoInitialize(&info);
71: info.fill = 5.0;
72: MatGetOrdering(A[0],MATORDERINGNATURAL,&perm,&iperm);
73: switch (facttype) {
74: case MAT_FACTOR_LU:
75: MatLUFactorSymbolic(F,A[0],perm,iperm,&info);
76: break;
77: case MAT_FACTOR_ILU:
78: MatILUFactorSymbolic(F,A[0],perm,iperm,&info);
79: break;
80: case MAT_FACTOR_ICC:
81: MatICCFactorSymbolic(F,A[0],perm,&info);
82: break;
83: case MAT_FACTOR_CHOLESKY:
84: MatCholeskyFactorSymbolic(F,A[0],perm,&info);
85: break;
86: default:
87: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not for factor type %s",factortype);
88: }
90: /* Compute numeric factors using same F, then solve */
91: for (k = 0; k < num_numfac; k++) {
92: switch (facttype) {
93: case MAT_FACTOR_LU:
94: case MAT_FACTOR_ILU:
95: MatLUFactorNumeric(F,A[k],&info);
96: break;
97: case MAT_FACTOR_ICC:
98: case MAT_FACTOR_CHOLESKY:
99: MatCholeskyFactorNumeric(F,A[k],&info);
100: break;
101: default:
102: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not for factor type %s",factortype);
103: }
105: /* Solve A[k] * x = b */
106: MatSolve(F,b,x);
108: /* Check the residual */
109: MatMult(A[k],x,u);
110: VecAXPY(u,-1.0,b);
111: VecNorm(u,NORM_INFINITY,&norm);
112: if (norm > tol) {
113: PetscPrintf(PETSC_COMM_WORLD,"%" PetscInt_FMT "-the %s numfact and solve: residual %g\n",k,factortype,(double)norm);
114: }
115: }
117: /* Free data structures */
118: for (k=0; k<num_numfac; k++) {
119: MatDestroy(&A[k]);
120: }
121: MatDestroy(&F);
122: ISDestroy(&perm);
123: ISDestroy(&iperm);
124: VecDestroy(&x);
125: VecDestroy(&b);
126: VecDestroy(&u);
127: PetscFinalize();
128: return 0;
129: }
131: /*TEST
133: test:
135: test:
136: suffix: 2
137: args: -mat_solver_type superlu
138: requires: superlu
140: test:
141: suffix: 3
142: nsize: 2
143: requires: mumps
144: args: -mat_solver_type mumps
146: test:
147: suffix: 4
148: args: -mat_solver_type cusparse -mat_type aijcusparse -mat_factor_type {{lu cholesky ilu icc}separate output}
149: requires: cuda
151: TEST*/