Actual source code: ex1.c
2: static char help[] = "Tests LU and Cholesky factorization for a dense matrix.\n\n";
4: #include petscmat.h
8: int main(int argc,char **argv)
9: {
10: Mat mat,fact;
11: MatInfo info;
13: PetscInt m = 10,n = 10,i = 4,rstart,rend;
14: PetscScalar value = 1.0;
15: Vec x,y,b;
16: PetscReal norm;
17: IS perm;
18: MatFactorInfo luinfo,factinfo;
20: PetscInitialize(&argc,&argv,(char*) 0,help);
22: VecCreate(PETSC_COMM_WORLD,&y);
23: VecSetSizes(y,PETSC_DECIDE,m);
24: VecSetFromOptions(y);
25: VecDuplicate(y,&x);
26: VecSet(x,value);
27: VecCreate(PETSC_COMM_WORLD,&b);
28: VecSetSizes(b,PETSC_DECIDE,n);
29: VecSetFromOptions(b);
30: ISCreateStride(PETSC_COMM_WORLD,m,0,1,&perm);
32: MatCreateSeqDense(PETSC_COMM_WORLD,m,n,PETSC_NULL,&mat);
34: MatGetOwnershipRange(mat,&rstart,&rend);
35: for (i=rstart; i<rend; i++) {
36: value = (PetscReal)i+1;
37: MatSetValues(mat,1,&i,1,&i,&value,INSERT_VALUES);
38: }
39: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
40: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
42: MatGetInfo(mat,MAT_LOCAL,&info);
43: PetscPrintf(PETSC_COMM_WORLD,"matrix nonzeros = %D, allocated nonzeros = %D\n",
44: (PetscInt)info.nz_used,(PetscInt)info.nz_allocated);
46: /* Cholesky factorization is not yet in place for this matrix format */
47: factinfo.fill = 1.0;
48: MatMult(mat,x,b);
49: MatConvert(mat,MATSAME,MAT_INITIAL_MATRIX,&fact);
50: MatCholeskyFactor(fact,perm,&factinfo);
51: MatSolve(fact,b,y);
52: MatDestroy(fact);
53: value = -1.0; VecAXPY(y,value,x);
54: VecNorm(y,NORM_2,&norm);
55: PetscPrintf(PETSC_COMM_WORLD,"Norm of error for Cholesky %A\n",norm);
56: MatGetFactor(mat,MAT_SOLVER_PETSC,MAT_FACTOR_CHOLESKY,&fact);
57: MatCholeskyFactorSymbolic(fact,mat,perm,&factinfo);
58: MatCholeskyFactorNumeric(fact,mat,&factinfo);
59: MatSolve(fact,b,y);
60: value = -1.0; VecAXPY(y,value,x);
61: VecNorm(y,NORM_2,&norm);
62: PetscPrintf(PETSC_COMM_WORLD,"Norm of error for Cholesky %A\n",norm);
63: MatDestroy(fact);
65: i = m-1; value = 1.0;
66: MatSetValues(mat,1,&i,1,&i,&value,INSERT_VALUES);
67: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
68: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
69: MatMult(mat,x,b);
70: MatConvert(mat,MATSAME,MAT_INITIAL_MATRIX,&fact);
71: luinfo.fill = 5.0;
72: luinfo.dtcol = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */
73: luinfo.shiftnz = 0.0;
74: luinfo.zeropivot = 1.e-12;
75: luinfo.pivotinblocks = 1.0;
76: luinfo.shiftpd = 0.0; /* false */
77: MatLUFactor(fact,perm,perm,&luinfo);
78: MatSolve(fact,b,y);
79: value = -1.0; VecAXPY(y,value,x);
80: VecNorm(y,NORM_2,&norm);
81: PetscPrintf(PETSC_COMM_WORLD,"Norm of error for LU %A\n",norm);
82: MatDestroy(fact);
84: luinfo.fill = 2.0;
85: luinfo.dtcol = 0.0;
86: luinfo.shiftnz = 0.0;
87: luinfo.zeropivot = 1.e-14;
88: luinfo.pivotinblocks = 1.0;
89: MatGetFactor(mat,MAT_SOLVER_PETSC,MAT_FACTOR_LU,&fact);
90: MatLUFactorSymbolic(fact,mat,perm,perm,&luinfo);
91: MatLUFactorNumeric(fact,mat,&luinfo);
92: MatSolve(fact,b,y);
93: value = -1.0; VecAXPY(y,value,x);
94: VecNorm(y,NORM_2,&norm);
95: PetscPrintf(PETSC_COMM_WORLD,"Norm of error for LU %A\n",norm);
96: MatDestroy(fact);
97: MatDestroy(mat);
98: VecDestroy(x);
99: VecDestroy(b);
100: VecDestroy(y);
101: ISDestroy(perm);
103: PetscFinalize();
104: return 0;
105: }
106: