Actual source code: ex28.c
2: static char help[] = "Tests MatReorderForNonzeroDiagonal()\n\n";
4: #include petscmat.h
8: int main(int argc,char **args)
9: {
10: Mat A,LU;
11: Vec x,y;
12: PetscInt nnz[4]={2,1,1,1},col[4],i;
14: PetscScalar values[4];
15: IS rowperm,colperm;
17: PetscInitialize(&argc,&args,(char *)0,help);
19: MatCreateSeqAIJ(PETSC_COMM_WORLD,4,4,2,nnz,&A);
21: /* build test matrix */
22: values[0]=1.0;values[1]=-1.0;
23: col[0]=0;col[1]=2; i=0;
24: MatSetValues(A,1,&i,2,col,values,INSERT_VALUES);
25: values[0]=1.0;
26: col[0]=1;i=1;
27: MatSetValues(A,1,&i,1,col,values,INSERT_VALUES);
28: values[0]=-1.0;
29: col[0]=3;i=2;
30: MatSetValues(A,1,&i,1,col,values,INSERT_VALUES);
31: values[0]=1.0;
32: col[0]=2;i=3;
33: MatSetValues(A,1,&i,1,col,values,INSERT_VALUES);
34: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
35: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
36: MatView(A,PETSC_VIEWER_STDOUT_SELF);
38: MatGetOrdering(A,MATORDERING_NATURAL,&rowperm,&colperm);
39: MatReorderForNonzeroDiagonal(A,1.e-12,rowperm,colperm);
40: PetscPrintf(PETSC_COMM_SELF,"column and row perms\n");
41: ISView(rowperm,0);
42: ISView(colperm,0);
43: MatGetFactor(A,MAT_SOLVER_PETSC,MAT_FACTOR_LU,&LU);
44: MatLUFactorSymbolic(LU,A,rowperm,colperm,PETSC_NULL);
45: MatLUFactorNumeric(LU,A,PETSC_NULL);
46: MatView(LU,PETSC_VIEWER_STDOUT_SELF);
47: VecCreate(PETSC_COMM_WORLD,&x);
48: VecSetSizes(x,PETSC_DECIDE,4);
49: VecSetFromOptions(x);
50: VecDuplicate(x,&y);
51: values[0]=0;values[1]=1.0;values[2]=-1.0;values[3]=1.0;
52: for (i=0; i<4; i++) col[i]=i;
53: VecSetValues(x,4,col,values,INSERT_VALUES);
54: VecAssemblyBegin(x);
55: VecAssemblyEnd(x);
56: VecView(x,PETSC_VIEWER_STDOUT_SELF);
58: MatSolve(LU,x,y);
59: VecView(y,PETSC_VIEWER_STDOUT_SELF);
61: ISDestroy(rowperm);
62: ISDestroy(colperm);
63: MatDestroy(LU);
64: MatDestroy(A);
65: VecDestroy(x);
66: VecDestroy(y);
67: PetscFinalize();
68: return 0;
69: }