Actual source code: ex68.c

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

  3: #include <petscmat.h>

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

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

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

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

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

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

 49:   PetscCall(MatGetOrdering(mat, MATORDERINGNATURAL, &isrow, &iscol));

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

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

 66:   PetscCall(ISDestroy(&isrow));
 67:   PetscCall(ISDestroy(&iscol));

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

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

 89:   PetscCall(ISDestroy(&isrow));
 90:   PetscCall(ISDestroy(&iscol));

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

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

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

134:   PetscCall(PetscFinalize());
135:   return 0;
136: }

138: /*TEST

140:    test:

142: TEST*/