2: static char help[] = "Tests MatReorderForNonzeroDiagonal().\n\n";

  4: #include <petscmat.h>

  8: int main(int argc,char **argv)
  9: {
 10:   Mat            mat,B;
 12:   PetscInt       i,j;
 13:   PetscScalar    v;
 14:   IS             isrow,iscol;
 15:   PetscViewer    viewer;

 17:   PetscInitialize(&argc,&argv,(char*)0,help);


 20:   /* ------- Assemble matrix, --------- */

 22:   MatCreate(PETSC_COMM_WORLD,&mat);
 23:   MatSetSizes(mat,PETSC_DECIDE,PETSC_DECIDE,4,4);
 24:   MatSetFromOptions(mat);
 25:   MatSetUp(mat);

 27:   /* set anti-diagonal of matrix */
 28:   v    = 1.0;
 29:   i    = 0; j = 3;
 30:   MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
 31:   v    = 2.0;
 32:   i    = 1; j = 2;
 33:   MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
 34:   v    = 3.0;
 35:   i    = 2; j = 1;
 36:   MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
 37:   v    = 4.0;
 38:   i    = 3; j = 0;
 39:   MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
 40:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
 41:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

 43:   PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer);
 44:   PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_DENSE);
 45:   PetscViewerASCIIPrintf(viewer,"Original matrix\n");
 46:   MatView(mat,viewer);

 48:   MatGetOrdering(mat,MATORDERINGNATURAL,&isrow,&iscol);

 50:   MatPermute(mat,isrow,iscol,&B);
 51:   PetscViewerASCIIPrintf(viewer,"Original matrix permuted by identity\n");
 52:   MatView(B,viewer);
 53:   MatDestroy(&B);

 55:   MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol);
 56:   MatPermute(mat,isrow,iscol,&B);
 57:   PetscViewerASCIIPrintf(viewer,"Original matrix permuted by identity + NonzeroDiagonal()\n");
 58:   MatView(B,viewer);
 59:   PetscViewerASCIIPrintf(viewer,"Row permutation\n");
 60:   ISView(isrow,viewer);
 61:   PetscViewerASCIIPrintf(viewer,"Column permutation\n");
 62:   ISView(iscol,viewer);
 63:   MatDestroy(&B);

 65:   ISDestroy(&isrow);
 66:   ISDestroy(&iscol);

 68:   MatGetOrdering(mat,MATORDERINGND,&isrow,&iscol);
 69:   MatPermute(mat,isrow,iscol,&B);
 70:   PetscViewerASCIIPrintf(viewer,"Original matrix permuted by ND\n");
 71:   MatView(B,viewer);
 72:   MatDestroy(&B);
 73:   PetscViewerASCIIPrintf(viewer,"ND row permutation\n");
 74:   ISView(isrow,viewer);
 75:   PetscViewerASCIIPrintf(viewer,"ND column permutation\n");
 76:   ISView(iscol,viewer);

 78:   MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol);
 79:   MatPermute(mat,isrow,iscol,&B);
 80:   PetscViewerASCIIPrintf(viewer,"Original matrix permuted by ND + NonzeroDiagonal()\n");
 81:   MatView(B,viewer);
 82:   MatDestroy(&B);
 83:   PetscViewerASCIIPrintf(viewer,"ND + NonzeroDiagonal() row permutation\n");
 84:   ISView(isrow,viewer);
 85:   PetscViewerASCIIPrintf(viewer,"ND + NonzeroDiagonal() column permutation\n");
 86:   ISView(iscol,viewer);

 88:   ISDestroy(&isrow);
 89:   ISDestroy(&iscol);

 91:   MatGetOrdering(mat,MATORDERINGRCM,&isrow,&iscol);
 92:   MatPermute(mat,isrow,iscol,&B);
 93:   PetscViewerASCIIPrintf(viewer,"Original matrix permuted by RCM\n");
 94:   MatView(B,viewer);
 95:   MatDestroy(&B);
 96:   PetscViewerASCIIPrintf(viewer,"RCM row permutation\n");
 97:   ISView(isrow,viewer);
 98:   PetscViewerASCIIPrintf(viewer,"RCM column permutation\n");
 99:   ISView(iscol,viewer);

101:   MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol);
102:   MatPermute(mat,isrow,iscol,&B);
103:   PetscViewerASCIIPrintf(viewer,"Original matrix permuted by RCM + NonzeroDiagonal()\n");
104:   MatView(B,viewer);
105:   MatDestroy(&B);
106:   PetscViewerASCIIPrintf(viewer,"RCM + NonzeroDiagonal() row permutation\n");
107:   ISView(isrow,viewer);
108:   PetscViewerASCIIPrintf(viewer,"RCM + NonzeroDiagonal() column permutation\n");
109:   ISView(iscol,viewer);

111:   /* Test MatLUFactor(); set diagonal as zeros as requested by PETSc matrix factorization */
112:   for (i=0; i<4; i++) {
113:     v = 0.0;
114:     MatSetValues(mat,1,&i,1,&i,&v,INSERT_VALUES);
115:   }
116:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
117:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
118:   MatLUFactor(mat,isrow,iscol,NULL);

120:   /* Free data structures */
121:   ISDestroy(&isrow);
122:   ISDestroy(&iscol);
123:   MatDestroy(&mat);

125:   PetscFinalize();
126:   return 0;
127: }