Actual source code: ex68.c
2: static char help[] = "Tests MatReorderForNonzeroDiagonal().\n\n";
4: #include <petscmat.h>
6: int main(int argc,char **argv)
7: {
8: Mat mat,B,C;
9: PetscInt i,j;
10: PetscMPIInt size;
11: PetscScalar v;
12: IS isrow,iscol,identity;
13: PetscViewer viewer;
15: PetscInitialize(&argc,&argv,(char*)0,help);
17: /* ------- Assemble matrix, --------- */
19: MatCreate(PETSC_COMM_WORLD,&mat);
20: MatSetSizes(mat,PETSC_DECIDE,PETSC_DECIDE,4,4);
21: MatSetFromOptions(mat);
22: MatSetUp(mat);
24: /* set anti-diagonal of matrix */
25: v = 1.0;
26: i = 0; j = 3;
27: MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
28: v = 2.0;
29: i = 1; j = 2;
30: MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
31: v = 3.0;
32: i = 2; j = 1;
33: MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
34: v = 4.0;
35: i = 3; j = 0;
36: MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
37: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
38: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
40: PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer);
41: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_DENSE);
42: PetscViewerASCIIPrintf(viewer,"Original matrix\n");
43: MatView(mat,viewer);
45: MatGetOrdering(mat,MATORDERINGNATURAL,&isrow,&iscol);
47: MatPermute(mat,isrow,iscol,&B);
48: PetscViewerASCIIPrintf(viewer,"Original matrix permuted by identity\n");
49: MatView(B,viewer);
50: MatDestroy(&B);
52: MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol);
53: MatPermute(mat,isrow,iscol,&B);
54: PetscViewerASCIIPrintf(viewer,"Original matrix permuted by identity + NonzeroDiagonal()\n");
55: MatView(B,viewer);
56: PetscViewerASCIIPrintf(viewer,"Row permutation\n");
57: ISView(isrow,viewer);
58: PetscViewerASCIIPrintf(viewer,"Column permutation\n");
59: ISView(iscol,viewer);
60: MatDestroy(&B);
62: ISDestroy(&isrow);
63: ISDestroy(&iscol);
65: MatGetOrdering(mat,MATORDERINGND,&isrow,&iscol);
66: MatPermute(mat,isrow,iscol,&B);
67: PetscViewerASCIIPrintf(viewer,"Original matrix permuted by ND\n");
68: MatView(B,viewer);
69: MatDestroy(&B);
70: PetscViewerASCIIPrintf(viewer,"ND row permutation\n");
71: ISView(isrow,viewer);
72: PetscViewerASCIIPrintf(viewer,"ND column permutation\n");
73: ISView(iscol,viewer);
75: MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol);
76: MatPermute(mat,isrow,iscol,&B);
77: PetscViewerASCIIPrintf(viewer,"Original matrix permuted by ND + NonzeroDiagonal()\n");
78: MatView(B,viewer);
79: MatDestroy(&B);
80: PetscViewerASCIIPrintf(viewer,"ND + NonzeroDiagonal() row permutation\n");
81: ISView(isrow,viewer);
82: PetscViewerASCIIPrintf(viewer,"ND + NonzeroDiagonal() column permutation\n");
83: ISView(iscol,viewer);
85: ISDestroy(&isrow);
86: ISDestroy(&iscol);
88: MatGetOrdering(mat,MATORDERINGRCM,&isrow,&iscol);
89: MatPermute(mat,isrow,iscol,&B);
90: PetscViewerASCIIPrintf(viewer,"Original matrix permuted by RCM\n");
91: MatView(B,viewer);
92: MatDestroy(&B);
93: PetscViewerASCIIPrintf(viewer,"RCM row permutation\n");
94: ISView(isrow,viewer);
95: PetscViewerASCIIPrintf(viewer,"RCM column permutation\n");
96: ISView(iscol,viewer);
98: MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol);
99: MatPermute(mat,isrow,iscol,&B);
100: PetscViewerASCIIPrintf(viewer,"Original matrix permuted by RCM + NonzeroDiagonal()\n");
101: MatView(B,viewer);
102: PetscViewerASCIIPrintf(viewer,"RCM + NonzeroDiagonal() row permutation\n");
103: ISView(isrow,viewer);
104: PetscViewerASCIIPrintf(viewer,"RCM + NonzeroDiagonal() column permutation\n");
105: ISView(iscol,viewer);
106: MPI_Comm_size(PETSC_COMM_WORLD,&size);
107: if (size == 1) {
108: MatSetOption(B,MAT_SYMMETRIC,PETSC_TRUE);
109: ISCreateStride(PETSC_COMM_SELF,4,0,1,&identity);
110: MatPermute(B,identity,identity,&C);
111: MatConvert(C,MATSEQSBAIJ,MAT_INPLACE_MATRIX,&C);
112: MatDestroy(&C);
113: ISDestroy(&identity);
114: }
115: MatDestroy(&B);
116: /* Test MatLUFactor(); set diagonal as zeros as requested by PETSc matrix factorization */
117: for (i=0; i<4; i++) {
118: v = 0.0;
119: MatSetValues(mat,1,&i,1,&i,&v,INSERT_VALUES);
120: }
121: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
122: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
123: MatLUFactor(mat,isrow,iscol,NULL);
125: /* Free data structures */
126: ISDestroy(&isrow);
127: ISDestroy(&iscol);
128: MatDestroy(&mat);
130: PetscFinalize();
131: return 0;
132: }
134: /*TEST
136: test:
138: TEST*/