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: break;
90: }
92: /* Compute numeric factors using same F, then solve */
93: for (k = 0; k < num_numfac; k++) {
94: switch (facttype) {
95: case MAT_FACTOR_LU:
96: case MAT_FACTOR_ILU:
97: MatLUFactorNumeric(F,A[k],&info);
98: break;
99: case MAT_FACTOR_ICC:
100: case MAT_FACTOR_CHOLESKY:
101: MatCholeskyFactorNumeric(F,A[k],&info);
102: break;
103: default:
104: SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not for factor type %s\n",factortype);
105: break;
106: }
108: /* Solve A[k] * x = b */
109: MatSolve(F,b,x);
111: /* Check the residual */
112: MatMult(A[k],x,u);
113: VecAXPY(u,-1.0,b);
114: VecNorm(u,NORM_INFINITY,&norm);
115: if (norm > tol) {
116: PetscPrintf(PETSC_COMM_WORLD,"%D-the %s numfact and solve: residual %g\n",k,factortype,(double)norm);
117: }
118: }
120: /* Free data structures */
121: for (k=0; k<num_numfac; k++) {
122: MatDestroy(&A[k]);
123: }
124: MatDestroy(&F);
125: ISDestroy(&perm);
126: ISDestroy(&iperm);
127: VecDestroy(&x);
128: VecDestroy(&b);
129: VecDestroy(&u);
130: PetscFinalize();
131: return ierr;
132: }
134: /*TEST
136: test:
138: test:
139: suffix: 2
140: args: -mat_solver_type superlu
141: requires: superlu
143: test:
144: suffix: 3
145: nsize: 2
146: requires: mumps
147: args: -mat_solver_type mumps
149: test:
150: suffix: 4
151: args: -mat_solver_type cusparse -mat_type aijcusparse -mat_factor_type {{lu cholesky ilu icc}separate output}
152: requires: cuda
154: TEST*/