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: MatFactorInfoInitialize(&factinfo);
48: factinfo.fill = 1.0;
49: MatMult(mat,x,b);
50: MatConvert(mat,MATSAME,MAT_INITIAL_MATRIX,&fact);
51: MatCholeskyFactor(fact,perm,&factinfo);
52: MatSolve(fact,b,y);
53: MatDestroy(fact);
54: value = -1.0; VecAXPY(y,value,x);
55: VecNorm(y,NORM_2,&norm);
56: PetscPrintf(PETSC_COMM_WORLD,"Norm of error for Cholesky %A\n",norm);
57: MatGetFactor(mat,MAT_SOLVER_PETSC,MAT_FACTOR_CHOLESKY,&fact);
58: MatCholeskyFactorSymbolic(fact,mat,perm,&factinfo);
59: MatCholeskyFactorNumeric(fact,mat,&factinfo);
60: MatSolve(fact,b,y);
61: value = -1.0; VecAXPY(y,value,x);
62: VecNorm(y,NORM_2,&norm);
63: PetscPrintf(PETSC_COMM_WORLD,"Norm of error for Cholesky %A\n",norm);
64: MatDestroy(fact);
66: i = m-1; value = 1.0;
67: MatSetValues(mat,1,&i,1,&i,&value,INSERT_VALUES);
68: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
69: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
70: MatMult(mat,x,b);
71: MatConvert(mat,MATSAME,MAT_INITIAL_MATRIX,&fact);
73: MatFactorInfoInitialize(&luinfo);
74: luinfo.fill = 1.0;
75: luinfo.dtcol = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */
76: luinfo.zeropivot = 1.e-12;
77: luinfo.shifttype = (PetscReal)MAT_SHIFT_INBLOCKS;
78: luinfo.shiftamount = 1.e-12;
79:
80: MatLUFactor(fact,perm,perm,&luinfo);
81: MatSolve(fact,b,y);
82: value = -1.0; VecAXPY(y,value,x);
83: VecNorm(y,NORM_2,&norm);
84: PetscPrintf(PETSC_COMM_WORLD,"Norm of error for LU %A\n",norm);
85: MatDestroy(fact);
87: luinfo.fill = 1.0;
88: luinfo.dtcol = 0.0;
89: luinfo.zeropivot = 1.e-14;
90: luinfo.shifttype = (PetscReal)MAT_SHIFT_INBLOCKS;
91: luinfo.shiftamount = 1.e-12;
92: MatGetFactor(mat,MAT_SOLVER_PETSC,MAT_FACTOR_LU,&fact);
93: MatLUFactorSymbolic(fact,mat,perm,perm,&luinfo);
94: MatLUFactorNumeric(fact,mat,&luinfo);
95: MatSolve(fact,b,y);
96: value = -1.0; VecAXPY(y,value,x);
97: VecNorm(y,NORM_2,&norm);
98: PetscPrintf(PETSC_COMM_WORLD,"Norm of error for LU %A\n",norm);
99: MatDestroy(fact);
100: MatDestroy(mat);
101: VecDestroy(x);
102: VecDestroy(b);
103: VecDestroy(y);
104: ISDestroy(perm);
105: PetscFinalize();
106: return 0;
107: }
108: