1 2 static char help[] = "Tests MatReorderForNonzeroDiagonal().\n\n"; 3 4 #include <petscmat.h> 5 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; 14 15 PetscFunctionBeginUser; 16 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 17 18 /* ------- Assemble matrix, --------- */ 19 20 PetscCall(MatCreate(PETSC_COMM_WORLD, &mat)); 21 PetscCall(MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, 4, 4)); 22 PetscCall(MatSetFromOptions(mat)); 23 PetscCall(MatSetUp(mat)); 24 25 /* set anti-diagonal of matrix */ 26 v = 1.0; 27 i = 0; 28 j = 3; 29 PetscCall(MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES)); 30 v = 2.0; 31 i = 1; 32 j = 2; 33 PetscCall(MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES)); 34 v = 3.0; 35 i = 2; 36 j = 1; 37 PetscCall(MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES)); 38 v = 4.0; 39 i = 3; 40 j = 0; 41 PetscCall(MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES)); 42 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 43 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 44 45 PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD, &viewer)); 46 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_DENSE)); 47 PetscCall(PetscViewerASCIIPrintf(viewer, "Original matrix\n")); 48 PetscCall(MatView(mat, viewer)); 49 50 PetscCall(MatGetOrdering(mat, MATORDERINGNATURAL, &isrow, &iscol)); 51 52 PetscCall(MatPermute(mat, isrow, iscol, &B)); 53 PetscCall(PetscViewerASCIIPrintf(viewer, "Original matrix permuted by identity\n")); 54 PetscCall(MatView(B, viewer)); 55 PetscCall(MatDestroy(&B)); 56 57 PetscCall(MatReorderForNonzeroDiagonal(mat, 1.e-8, isrow, iscol)); 58 PetscCall(MatPermute(mat, isrow, iscol, &B)); 59 PetscCall(PetscViewerASCIIPrintf(viewer, "Original matrix permuted by identity + NonzeroDiagonal()\n")); 60 PetscCall(MatView(B, viewer)); 61 PetscCall(PetscViewerASCIIPrintf(viewer, "Row permutation\n")); 62 PetscCall(ISView(isrow, viewer)); 63 PetscCall(PetscViewerASCIIPrintf(viewer, "Column permutation\n")); 64 PetscCall(ISView(iscol, viewer)); 65 PetscCall(MatDestroy(&B)); 66 67 PetscCall(ISDestroy(&isrow)); 68 PetscCall(ISDestroy(&iscol)); 69 70 PetscCall(MatGetOrdering(mat, MATORDERINGND, &isrow, &iscol)); 71 PetscCall(MatPermute(mat, isrow, iscol, &B)); 72 PetscCall(PetscViewerASCIIPrintf(viewer, "Original matrix permuted by ND\n")); 73 PetscCall(MatView(B, viewer)); 74 PetscCall(MatDestroy(&B)); 75 PetscCall(PetscViewerASCIIPrintf(viewer, "ND row permutation\n")); 76 PetscCall(ISView(isrow, viewer)); 77 PetscCall(PetscViewerASCIIPrintf(viewer, "ND column permutation\n")); 78 PetscCall(ISView(iscol, viewer)); 79 80 PetscCall(MatReorderForNonzeroDiagonal(mat, 1.e-8, isrow, iscol)); 81 PetscCall(MatPermute(mat, isrow, iscol, &B)); 82 PetscCall(PetscViewerASCIIPrintf(viewer, "Original matrix permuted by ND + NonzeroDiagonal()\n")); 83 PetscCall(MatView(B, viewer)); 84 PetscCall(MatDestroy(&B)); 85 PetscCall(PetscViewerASCIIPrintf(viewer, "ND + NonzeroDiagonal() row permutation\n")); 86 PetscCall(ISView(isrow, viewer)); 87 PetscCall(PetscViewerASCIIPrintf(viewer, "ND + NonzeroDiagonal() column permutation\n")); 88 PetscCall(ISView(iscol, viewer)); 89 90 PetscCall(ISDestroy(&isrow)); 91 PetscCall(ISDestroy(&iscol)); 92 93 PetscCall(MatGetOrdering(mat, MATORDERINGRCM, &isrow, &iscol)); 94 PetscCall(MatPermute(mat, isrow, iscol, &B)); 95 PetscCall(PetscViewerASCIIPrintf(viewer, "Original matrix permuted by RCM\n")); 96 PetscCall(MatView(B, viewer)); 97 PetscCall(MatDestroy(&B)); 98 PetscCall(PetscViewerASCIIPrintf(viewer, "RCM row permutation\n")); 99 PetscCall(ISView(isrow, viewer)); 100 PetscCall(PetscViewerASCIIPrintf(viewer, "RCM column permutation\n")); 101 PetscCall(ISView(iscol, viewer)); 102 103 PetscCall(MatReorderForNonzeroDiagonal(mat, 1.e-8, isrow, iscol)); 104 PetscCall(MatPermute(mat, isrow, iscol, &B)); 105 PetscCall(PetscViewerASCIIPrintf(viewer, "Original matrix permuted by RCM + NonzeroDiagonal()\n")); 106 PetscCall(MatView(B, viewer)); 107 PetscCall(PetscViewerASCIIPrintf(viewer, "RCM + NonzeroDiagonal() row permutation\n")); 108 PetscCall(ISView(isrow, viewer)); 109 PetscCall(PetscViewerASCIIPrintf(viewer, "RCM + NonzeroDiagonal() column permutation\n")); 110 PetscCall(ISView(iscol, viewer)); 111 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 112 if (size == 1) { 113 PetscCall(MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE)); 114 PetscCall(ISCreateStride(PETSC_COMM_SELF, 4, 0, 1, &identity)); 115 PetscCall(MatPermute(B, identity, identity, &C)); 116 PetscCall(MatConvert(C, MATSEQSBAIJ, MAT_INPLACE_MATRIX, &C)); 117 PetscCall(MatDestroy(&C)); 118 PetscCall(ISDestroy(&identity)); 119 } 120 PetscCall(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 PetscCall(MatSetValues(mat, 1, &i, 1, &i, &v, INSERT_VALUES)); 125 } 126 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 127 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 128 PetscCall(MatLUFactor(mat, isrow, iscol, NULL)); 129 130 /* Free data structures */ 131 PetscCall(ISDestroy(&isrow)); 132 PetscCall(ISDestroy(&iscol)); 133 PetscCall(MatDestroy(&mat)); 134 135 PetscCall(PetscFinalize()); 136 return 0; 137 } 138 139 /*TEST 140 141 test: 142 143 TEST*/ 144