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;

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


 19:   /* ------- Assemble matrix, --------- */

 21:   MatCreateSeqAIJ(PETSC_COMM_WORLD,4,4,0,0,&mat);

 23:   /* set anti-diagonal of matrix */
 24:   v = 1.0;
 25:   i = 0; j = 3;
 26:   MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
 27:   v = 2.0;
 28:   i = 1; j = 2;
 29:   MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
 30:   v = 3.0;
 31:   i = 2; j = 1;
 32:   MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
 33:   v = 4.0;
 34:   i = 3; j = 0;
 35:   MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);

 37:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
 38:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

 40:   printf("Original matrix\n");
 41:   PetscViewerSetFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_DENSE);
 42:   MatView(mat,PETSC_VIEWER_STDOUT_SELF);

 44:   MatGetOrdering(mat,MATORDERING_NATURAL,&isrow,&iscol);

 46:   MatPermute(mat,isrow,iscol,&B);
 47:   printf("Original matrix permuted by identity\n");
 48:   MatView(B,PETSC_VIEWER_STDOUT_SELF);
 49:   MatDestroy(B);

 51:   MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol);
 52:   MatPermute(mat,isrow,iscol,&B);
 53:   printf("Original matrix permuted by identity + NonzeroDiagonal()\n");
 54:   MatView(B,PETSC_VIEWER_STDOUT_SELF);
 55:   printf("Row permutation\n");
 56:   ISView(isrow,PETSC_VIEWER_STDOUT_SELF);
 57:   printf("Column permutation\n");
 58:   ISView(iscol,PETSC_VIEWER_STDOUT_SELF);
 59:   MatDestroy(B);

 61:   ISDestroy(isrow);
 62:   ISDestroy(iscol);

 64:   MatGetOrdering(mat,MATORDERING_ND,&isrow,&iscol);
 65:   MatPermute(mat,isrow,iscol,&B);
 66:   printf("Original matrix permuted by ND\n");
 67:   MatView(B,PETSC_VIEWER_STDOUT_SELF);
 68:   MatDestroy(B);
 69:   printf("ND row permutation\n");
 70:   ISView(isrow,PETSC_VIEWER_STDOUT_SELF);
 71:   printf("ND column permutation\n");
 72:   ISView(iscol,PETSC_VIEWER_STDOUT_SELF);

 74:   MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol);
 75:   MatPermute(mat,isrow,iscol,&B);
 76:   printf("Original matrix permuted by ND + NonzeroDiagonal()\n");
 77:   MatView(B,PETSC_VIEWER_STDOUT_SELF);
 78:   MatDestroy(B);
 79:   printf("ND + NonzeroDiagonal() row permutation\n");
 80:   ISView(isrow,PETSC_VIEWER_STDOUT_SELF);
 81:   printf("ND + NonzeroDiagonal() column permutation\n");
 82:   ISView(iscol,PETSC_VIEWER_STDOUT_SELF);

 84:   ISDestroy(isrow);
 85:   ISDestroy(iscol);

 87:   MatGetOrdering(mat,MATORDERING_RCM,&isrow,&iscol);
 88:   MatPermute(mat,isrow,iscol,&B);
 89:   printf("Original matrix permuted by RCM\n");
 90:   MatView(B,PETSC_VIEWER_STDOUT_SELF);
 91:   MatDestroy(B);
 92:   printf("RCM row permutation\n");
 93:   ISView(isrow,PETSC_VIEWER_STDOUT_SELF);
 94:   printf("RCM column permutation\n");
 95:   ISView(iscol,PETSC_VIEWER_STDOUT_SELF);

 97:   MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol);
 98:   MatPermute(mat,isrow,iscol,&B);
 99:   printf("Original matrix permuted by RCM + NonzeroDiagonal()\n");
100:   MatView(B,PETSC_VIEWER_STDOUT_SELF);
101:   MatDestroy(B);
102:   printf("RCM + NonzeroDiagonal() row permutation\n");
103:   ISView(isrow,PETSC_VIEWER_STDOUT_SELF);
104:   printf("RCM + NonzeroDiagonal() column permutation\n");
105:   ISView(iscol,PETSC_VIEWER_STDOUT_SELF);

107:    MatLUFactor(mat,isrow,iscol,PETSC_NULL);
108:   printf("Factored matrix permuted by RCM + NonzeroDiagonal()\n");
109:   MatView(mat,PETSC_VIEWER_STDOUT_SELF);

111:   /* Free data structures */
112:   ISDestroy(isrow);
113:   ISDestroy(iscol);
114:   MatDestroy(mat);

116:   PetscFinalize();
117:   return 0;
118: }
119: