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;
10: PetscInt i,j;
11: PetscMPIInt size;
12: PetscScalar v;
13: IS isrow,iscol,identity;
14: PetscViewer viewer;
16: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
18: /* ------- Assemble matrix, --------- */
20: MatCreate(PETSC_COMM_WORLD,&mat);
21: MatSetSizes(mat,PETSC_DECIDE,PETSC_DECIDE,4,4);
22: MatSetFromOptions(mat);
23: MatSetUp(mat);
25: /* set anti-diagonal of matrix */
26: v = 1.0;
27: i = 0; j = 3;
28: MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
29: v = 2.0;
30: i = 1; j = 2;
31: MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
32: v = 3.0;
33: i = 2; j = 1;
34: MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
35: v = 4.0;
36: i = 3; j = 0;
37: MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
38: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
39: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
41: PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer);
42: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_DENSE);
43: PetscViewerASCIIPrintf(viewer,"Original matrix\n");
44: MatView(mat,viewer);
46: MatGetOrdering(mat,MATORDERINGNATURAL,&isrow,&iscol);
48: MatPermute(mat,isrow,iscol,&B);
49: PetscViewerASCIIPrintf(viewer,"Original matrix permuted by identity\n");
50: MatView(B,viewer);
51: MatDestroy(&B);
53: MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol);
54: MatPermute(mat,isrow,iscol,&B);
55: PetscViewerASCIIPrintf(viewer,"Original matrix permuted by identity + NonzeroDiagonal()\n");
56: MatView(B,viewer);
57: PetscViewerASCIIPrintf(viewer,"Row permutation\n");
58: ISView(isrow,viewer);
59: PetscViewerASCIIPrintf(viewer,"Column permutation\n");
60: ISView(iscol,viewer);
61: MatDestroy(&B);
63: ISDestroy(&isrow);
64: ISDestroy(&iscol);
66: MatGetOrdering(mat,MATORDERINGND,&isrow,&iscol);
67: MatPermute(mat,isrow,iscol,&B);
68: PetscViewerASCIIPrintf(viewer,"Original matrix permuted by ND\n");
69: MatView(B,viewer);
70: MatDestroy(&B);
71: PetscViewerASCIIPrintf(viewer,"ND row permutation\n");
72: ISView(isrow,viewer);
73: PetscViewerASCIIPrintf(viewer,"ND column permutation\n");
74: ISView(iscol,viewer);
76: MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol);
77: MatPermute(mat,isrow,iscol,&B);
78: PetscViewerASCIIPrintf(viewer,"Original matrix permuted by ND + NonzeroDiagonal()\n");
79: MatView(B,viewer);
80: MatDestroy(&B);
81: PetscViewerASCIIPrintf(viewer,"ND + NonzeroDiagonal() row permutation\n");
82: ISView(isrow,viewer);
83: PetscViewerASCIIPrintf(viewer,"ND + NonzeroDiagonal() column permutation\n");
84: ISView(iscol,viewer);
86: ISDestroy(&isrow);
87: ISDestroy(&iscol);
89: MatGetOrdering(mat,MATORDERINGRCM,&isrow,&iscol);
90: MatPermute(mat,isrow,iscol,&B);
91: PetscViewerASCIIPrintf(viewer,"Original matrix permuted by RCM\n");
92: MatView(B,viewer);
93: MatDestroy(&B);
94: PetscViewerASCIIPrintf(viewer,"RCM row permutation\n");
95: ISView(isrow,viewer);
96: PetscViewerASCIIPrintf(viewer,"RCM column permutation\n");
97: ISView(iscol,viewer);
99: MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol);
100: MatPermute(mat,isrow,iscol,&B);
101: PetscViewerASCIIPrintf(viewer,"Original matrix permuted by RCM + NonzeroDiagonal()\n");
102: MatView(B,viewer);
103: PetscViewerASCIIPrintf(viewer,"RCM + NonzeroDiagonal() row permutation\n");
104: ISView(isrow,viewer);
105: PetscViewerASCIIPrintf(viewer,"RCM + NonzeroDiagonal() column permutation\n");
106: ISView(iscol,viewer);
107: MPI_Comm_size(PETSC_COMM_WORLD,&size);
108: if (size == 1) {
109: MatSetOption(B,MAT_SYMMETRIC,PETSC_TRUE);
110: ISCreateStride(PETSC_COMM_SELF,4,0,1,&identity);
111: MatPermute(B,identity,identity,&C);
112: MatConvert(C,MATSEQSBAIJ,MAT_INPLACE_MATRIX,&C);
113: MatDestroy(&C);
114: ISDestroy(&identity);
115: }
116: MatDestroy(&B);
117: /* Test MatLUFactor(); set diagonal as zeros as requested by PETSc matrix factorization */
118: for (i=0; i<4; i++) {
119: v = 0.0;
120: MatSetValues(mat,1,&i,1,&i,&v,INSERT_VALUES);
121: }
122: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
123: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
124: MatLUFactor(mat,isrow,iscol,NULL);
126: /* Free data structures */
127: ISDestroy(&isrow);
128: ISDestroy(&iscol);
129: MatDestroy(&mat);
131: PetscFinalize();
132: return ierr;
133: }
135: /*TEST
137: test:
139: TEST*/