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: MatLUFactorSymbolic(A,rowperm,colperm,PETSC_NULL,&LU);
44: MatLUFactorNumeric(A,PETSC_NULL,&LU);
45: MatView(LU,PETSC_VIEWER_STDOUT_SELF);
46: VecCreate(PETSC_COMM_WORLD,&x);
47: VecSetSizes(x,PETSC_DECIDE,4);
48: VecSetFromOptions(x);
49: VecDuplicate(x,&y);
50: values[0]=0;values[1]=1.0;values[2]=-1.0;values[3]=1.0;
51: for (i=0; i<4; i++) col[i]=i;
52: VecSetValues(x,4,col,values,INSERT_VALUES);
53: VecAssemblyBegin(x);
54: VecAssemblyEnd(x);
55: VecView(x,PETSC_VIEWER_STDOUT_SELF);
57: MatSolve(LU,x,y);
58: VecView(y,PETSC_VIEWER_STDOUT_SELF);
60: ISDestroy(rowperm);
61: ISDestroy(colperm);
62: MatDestroy(LU);
63: MatDestroy(A);
64: VecDestroy(x);
65: VecDestroy(y);
66: PetscFinalize();
67: return 0;
68: }