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;
  9:   PetscInt       i,j;
 10:   PetscMPIInt    size;
 11:   PetscScalar    v;
 12:   IS             isrow,iscol,identity;
 13:   PetscViewer    viewer;

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

 17:   /* ------- Assemble matrix, --------- */

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

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

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

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

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

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

 62:   ISDestroy(&isrow);
 63:   ISDestroy(&iscol);

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

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

 85:   ISDestroy(&isrow);
 86:   ISDestroy(&iscol);

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

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

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

130:   PetscFinalize();
131:   return 0;
132: }

134: /*TEST

136:    test:

138: TEST*/