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

  4:  #include <petscmat.h>

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

 15:   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:   MatDestroy(&B);
104:   PetscViewerASCIIPrintf(viewer,"RCM + NonzeroDiagonal() row permutation\n");
105:   ISView(isrow,viewer);
106:   PetscViewerASCIIPrintf(viewer,"RCM + NonzeroDiagonal() column permutation\n");
107:   ISView(iscol,viewer);

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

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

123:   PetscFinalize();
124:   return ierr;
125: }