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;


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

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

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

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

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

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

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

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

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

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

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

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

100:   MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol);
101:   MatPermute(mat,isrow,iscol,&B);
102:   PetscViewerASCIIPrintf(viewer,"Original matrix permuted by RCM + NonzeroDiagonal()\n");
103:   MatView(B,viewer);
104:   PetscViewerASCIIPrintf(viewer,"RCM + NonzeroDiagonal() row permutation\n");
105:   ISView(isrow,viewer);
106:   PetscViewerASCIIPrintf(viewer,"RCM + NonzeroDiagonal() column permutation\n");
107:   ISView(iscol,viewer);
108:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
109:   if (size == 1) {
110:     MatSetOption(B,MAT_SYMMETRIC,PETSC_TRUE);
111:     ISCreateStride(PETSC_COMM_SELF,4,0,1,&identity);
112:     MatPermute(B,identity,identity,&C);
113:     MatConvert(C,MATSEQSBAIJ,MAT_INPLACE_MATRIX,&C);
114:     MatDestroy(&C);
115:     ISDestroy(&identity);
116:   }
117:   MatDestroy(&B);
118:   /* Test MatLUFactor(); set diagonal as zeros as requested by PETSc matrix factorization */
119:   for (i=0; i<4; i++) {
120:     v = 0.0;
121:     MatSetValues(mat,1,&i,1,&i,&v,INSERT_VALUES);
122:   }
123:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
124:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
125:   MatLUFactor(mat,isrow,iscol,NULL);

127:   /* Free data structures */
128:   ISDestroy(&isrow);
129:   ISDestroy(&iscol);
130:   MatDestroy(&mat);

132:   PetscFinalize();
133:   return ierr;
134: }



138: /*TEST

140:    test:

142: TEST*/