Actual source code: ex7.c
petsc-3.4.2 2013-07-02
2: static char help[] = "Tests matrix factorization. Note that most users should\n\
3: employ the KSP interface to the linear solvers instead of using the factorization\n\
4: routines directly.\n\n";
6: #include <petscmat.h>
10: int main(int argc,char **args)
11: {
12: Mat C,LU;
13: MatInfo info;
14: PetscInt i,j,m = 3,n = 3,Ii,J;
16: PetscScalar v,one = 1.0;
17: IS perm,iperm;
18: Vec x,u,b;
19: PetscReal norm;
20: MatFactorInfo luinfo;
22: PetscInitialize(&argc,&args,(char*)0,help);
24: /* Create the matrix for the five point stencil, YET AGAIN */
25: MatCreate(PETSC_COMM_SELF,&C);
26: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
27: MatSetFromOptions(C);
28: MatSetUp(C);
29: for (i=0; i<m; i++) {
30: for (j=0; j<n; j++) {
31: v = -1.0; Ii = j + n*i;
32: if (i>0) {J = Ii - n; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
33: if (i<m-1) {J = Ii + n; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
34: if (j>0) {J = Ii - 1; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
35: if (j<n-1) {J = Ii + 1; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
36: v = 4.0; MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES);
37: }
38: }
39: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
40: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
41: MatGetOrdering(C,MATORDERINGRCM,&perm,&iperm);
42: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
43: ISView(perm,PETSC_VIEWER_STDOUT_SELF);
45: MatFactorInfoInitialize(&luinfo);
47: luinfo.fill = 2.0;
48: luinfo.dtcol = 0.0;
49: luinfo.zeropivot = 1.e-14;
50: luinfo.pivotinblocks = 1.0;
52: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_LU,&LU);
53: MatLUFactorSymbolic(LU,C,perm,iperm,&luinfo);
54: MatLUFactorNumeric(LU,C,&luinfo);
55: MatView(LU,PETSC_VIEWER_STDOUT_WORLD);
57: VecCreateSeq(PETSC_COMM_SELF,m*n,&u);
58: VecSet(u,one);
59: VecDuplicate(u,&x);
60: VecDuplicate(u,&b);
62: MatMult(C,u,b);
63: MatSolve(LU,b,x);
65: VecView(b,PETSC_VIEWER_STDOUT_SELF);
66: VecView(x,PETSC_VIEWER_STDOUT_SELF);
68: VecAXPY(x,-1.0,u);
69: VecNorm(x,NORM_2,&norm);
70: PetscPrintf(PETSC_COMM_SELF,"Norm of error %G\n",norm);
72: MatGetInfo(C,MAT_LOCAL,&info);
73: PetscPrintf(PETSC_COMM_SELF,"original matrix nonzeros = %D\n",(PetscInt)info.nz_used);
74: MatGetInfo(LU,MAT_LOCAL,&info);
75: PetscPrintf(PETSC_COMM_SELF,"factored matrix nonzeros = %D\n",(PetscInt)info.nz_used);
77: VecDestroy(&u);
78: VecDestroy(&b);
79: VecDestroy(&x);
80: ISDestroy(&perm);
81: ISDestroy(&iperm);
82: MatDestroy(&C);
83: MatDestroy(&LU);
85: PetscFinalize();
86: return 0;
87: }