Actual source code: ex30.c
2: static char help[] = "Tests ILU factorization and illustrates drawing of matrix sparsity structure with MatView().\n\
3: Input parameters are:\n\
4: -lf <level> : level of fill for ILU (default is 0)\n\
5: -lu : use full LU factorization\n\
6: -m <value>,-n <value> : grid dimensions\n\
7: Note that most users should employ the KSP interface to the\n\
8: linear solvers instead of using the factorization routines\n\
9: directly.\n\n";
11: #include petscmat.h
15: int main(int argc,char **args)
16: {
17: Mat C,A;
18: PetscInt i,j,m = 5,n = 5,Ii,J,lf = 0;
20: PetscTruth flg1;
21: PetscScalar v;
22: IS row,col;
23: PetscViewer viewer1,viewer2;
24: MatFactorInfo info;
26: PetscInitialize(&argc,&args,(char *)0,help);
27: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
28: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
29: PetscOptionsGetInt(PETSC_NULL,"-lf",&lf,PETSC_NULL);
31: PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,0,0,400,400,&viewer1);
32: PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,400,0,400,400,&viewer2);
34: MatCreate(PETSC_COMM_SELF,&C);
35: MatSetSizes(C,m*n,m*n,m*n,m*n);
36: MatSetFromOptions(C);
37: MatSeqBDiagSetPreallocation(C,0,1,PETSC_NULL,PETSC_NULL);
38: MatSeqAIJSetPreallocation(C,5,PETSC_NULL);
40: /* Create the matrix. (This is five-point stencil with some extra elements) */
41: for (i=0; i<m; i++) {
42: for (j=0; j<n; j++) {
43: v = -1.0; Ii = j + n*i;
44: J = Ii - n; if (J>=0) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
45: J = Ii + n; if (J<m*n) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
46: J = Ii - 1; if (J>=0) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
47: J = Ii + 1; if (J<m*n) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
48: v = 4.0; MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES);
49: }
50: }
51: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
52: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
54: MatGetOrdering(C,MATORDERING_RCM,&row,&col);
55: printf("original matrix:\n");
56: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO);
57: MatView(C,PETSC_VIEWER_STDOUT_SELF);
58: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
59: MatView(C,PETSC_VIEWER_STDOUT_SELF);
60: MatView(C,viewer1);
62: /* Compute factorization */
63: PetscOptionsHasName(PETSC_NULL,"-lu",&flg1);
64: if (flg1){
65: MatLUFactorSymbolic(C,row,col,PETSC_NULL,&A);
66: } else {
67: MatFactorInfoInitialize(&info);
68: info.levels = lf;
69: info.fill = 1.0;
70: info.diagonal_fill = 0;
71: info.shiftnz = 0;
72: info.zeropivot = 0.0;
73: MatILUFactorSymbolic(C,row,col,&info,&A);
74: }
75: MatLUFactorNumeric(C,&info,&A);
77: printf("factored matrix:\n");
78: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO);
79: MatView(A,PETSC_VIEWER_STDOUT_SELF);
80: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
81: MatView(A,PETSC_VIEWER_STDOUT_SELF);
82: MatView(A,viewer2);
84: /* Free data structures */
85: MatDestroy(C);
86: MatDestroy(A);
87: ISDestroy(row);
88: ISDestroy(col);
89: PetscViewerDestroy(viewer1);
90: PetscViewerDestroy(viewer2);
91: PetscFinalize();
92: return 0;
93: }