1 static char help[] = "Tests MatReorderForNonzeroDiagonal().\n\n"; 2 3 #include <petscmat.h> 4 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; 13 14 PetscFunctionBeginUser; 15 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 16 17 /* ------- Assemble matrix, --------- */ 18 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)); 23 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)); 43 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)); 48 49 PetscCall(MatGetOrdering(mat, MATORDERINGNATURAL, &isrow, &iscol)); 50 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)); 55 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)); 65 66 PetscCall(ISDestroy(&isrow)); 67 PetscCall(ISDestroy(&iscol)); 68 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)); 78 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)); 88 89 PetscCall(ISDestroy(&isrow)); 90 PetscCall(ISDestroy(&iscol)); 91 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)); 101 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)); 128 129 /* Free data structures */ 130 PetscCall(ISDestroy(&isrow)); 131 PetscCall(ISDestroy(&iscol)); 132 PetscCall(MatDestroy(&mat)); 133 134 PetscCall(PetscFinalize()); 135 return 0; 136 } 137 138 /*TEST 139 140 test: 141 142 TEST*/ 143