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;

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

 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;
 28:   j = 3;
 29:   MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES);
 30:   v = 2.0;
 31:   i = 1;
 32:   j = 2;
 33:   MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES);
 34:   v = 3.0;
 35:   i = 2;
 36:   j = 1;
 37:   MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES);
 38:   v = 4.0;
 39:   i = 3;
 40:   j = 0;
 41:   MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES);
 42:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
 43:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);

 45:   PetscViewerASCIIGetStdout(PETSC_COMM_WORLD, &viewer);
 46:   PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_DENSE);
 47:   PetscViewerASCIIPrintf(viewer, "Original matrix\n");
 48:   MatView(mat, viewer);

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

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

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

 67:   ISDestroy(&isrow);
 68:   ISDestroy(&iscol);

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

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

 90:   ISDestroy(&isrow);
 91:   ISDestroy(&iscol);

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

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

130:   /* Free data structures */
131:   ISDestroy(&isrow);
132:   ISDestroy(&iscol);
133:   MatDestroy(&mat);

135:   PetscFinalize();
136:   return 0;
137: }

139: /*TEST

141:    test:

143: TEST*/