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; j = 3; 28 PetscCall(MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES)); 29 v = 2.0; 30 i = 1; j = 2; 31 PetscCall(MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES)); 32 v = 3.0; 33 i = 2; j = 1; 34 PetscCall(MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES)); 35 v = 4.0; 36 i = 3; j = 0; 37 PetscCall(MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES)); 38 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 39 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 40 41 PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer)); 42 PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_DENSE)); 43 PetscCall(PetscViewerASCIIPrintf(viewer,"Original matrix\n")); 44 PetscCall(MatView(mat,viewer)); 45 46 PetscCall(MatGetOrdering(mat,MATORDERINGNATURAL,&isrow,&iscol)); 47 48 PetscCall(MatPermute(mat,isrow,iscol,&B)); 49 PetscCall(PetscViewerASCIIPrintf(viewer,"Original matrix permuted by identity\n")); 50 PetscCall(MatView(B,viewer)); 51 PetscCall(MatDestroy(&B)); 52 53 PetscCall(MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol)); 54 PetscCall(MatPermute(mat,isrow,iscol,&B)); 55 PetscCall(PetscViewerASCIIPrintf(viewer,"Original matrix permuted by identity + NonzeroDiagonal()\n")); 56 PetscCall(MatView(B,viewer)); 57 PetscCall(PetscViewerASCIIPrintf(viewer,"Row permutation\n")); 58 PetscCall(ISView(isrow,viewer)); 59 PetscCall(PetscViewerASCIIPrintf(viewer,"Column permutation\n")); 60 PetscCall(ISView(iscol,viewer)); 61 PetscCall(MatDestroy(&B)); 62 63 PetscCall(ISDestroy(&isrow)); 64 PetscCall(ISDestroy(&iscol)); 65 66 PetscCall(MatGetOrdering(mat,MATORDERINGND,&isrow,&iscol)); 67 PetscCall(MatPermute(mat,isrow,iscol,&B)); 68 PetscCall(PetscViewerASCIIPrintf(viewer,"Original matrix permuted by ND\n")); 69 PetscCall(MatView(B,viewer)); 70 PetscCall(MatDestroy(&B)); 71 PetscCall(PetscViewerASCIIPrintf(viewer,"ND row permutation\n")); 72 PetscCall(ISView(isrow,viewer)); 73 PetscCall(PetscViewerASCIIPrintf(viewer,"ND column permutation\n")); 74 PetscCall(ISView(iscol,viewer)); 75 76 PetscCall(MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol)); 77 PetscCall(MatPermute(mat,isrow,iscol,&B)); 78 PetscCall(PetscViewerASCIIPrintf(viewer,"Original matrix permuted by ND + NonzeroDiagonal()\n")); 79 PetscCall(MatView(B,viewer)); 80 PetscCall(MatDestroy(&B)); 81 PetscCall(PetscViewerASCIIPrintf(viewer,"ND + NonzeroDiagonal() row permutation\n")); 82 PetscCall(ISView(isrow,viewer)); 83 PetscCall(PetscViewerASCIIPrintf(viewer,"ND + NonzeroDiagonal() column permutation\n")); 84 PetscCall(ISView(iscol,viewer)); 85 86 PetscCall(ISDestroy(&isrow)); 87 PetscCall(ISDestroy(&iscol)); 88 89 PetscCall(MatGetOrdering(mat,MATORDERINGRCM,&isrow,&iscol)); 90 PetscCall(MatPermute(mat,isrow,iscol,&B)); 91 PetscCall(PetscViewerASCIIPrintf(viewer,"Original matrix permuted by RCM\n")); 92 PetscCall(MatView(B,viewer)); 93 PetscCall(MatDestroy(&B)); 94 PetscCall(PetscViewerASCIIPrintf(viewer,"RCM row permutation\n")); 95 PetscCall(ISView(isrow,viewer)); 96 PetscCall(PetscViewerASCIIPrintf(viewer,"RCM column permutation\n")); 97 PetscCall(ISView(iscol,viewer)); 98 99 PetscCall(MatReorderForNonzeroDiagonal(mat,1.e-8,isrow,iscol)); 100 PetscCall(MatPermute(mat,isrow,iscol,&B)); 101 PetscCall(PetscViewerASCIIPrintf(viewer,"Original matrix permuted by RCM + NonzeroDiagonal()\n")); 102 PetscCall(MatView(B,viewer)); 103 PetscCall(PetscViewerASCIIPrintf(viewer,"RCM + NonzeroDiagonal() row permutation\n")); 104 PetscCall(ISView(isrow,viewer)); 105 PetscCall(PetscViewerASCIIPrintf(viewer,"RCM + NonzeroDiagonal() column permutation\n")); 106 PetscCall(ISView(iscol,viewer)); 107 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 108 if (size == 1) { 109 PetscCall(MatSetOption(B,MAT_SYMMETRIC,PETSC_TRUE)); 110 PetscCall(ISCreateStride(PETSC_COMM_SELF,4,0,1,&identity)); 111 PetscCall(MatPermute(B,identity,identity,&C)); 112 PetscCall(MatConvert(C,MATSEQSBAIJ,MAT_INPLACE_MATRIX,&C)); 113 PetscCall(MatDestroy(&C)); 114 PetscCall(ISDestroy(&identity)); 115 } 116 PetscCall(MatDestroy(&B)); 117 /* Test MatLUFactor(); set diagonal as zeros as requested by PETSc matrix factorization */ 118 for (i=0; i<4; i++) { 119 v = 0.0; 120 PetscCall(MatSetValues(mat,1,&i,1,&i,&v,INSERT_VALUES)); 121 } 122 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 123 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 124 PetscCall(MatLUFactor(mat,isrow,iscol,NULL)); 125 126 /* Free data structures */ 127 PetscCall(ISDestroy(&isrow)); 128 PetscCall(ISDestroy(&iscol)); 129 PetscCall(MatDestroy(&mat)); 130 131 PetscCall(PetscFinalize()); 132 return 0; 133 } 134 135 /*TEST 136 137 test: 138 139 TEST*/ 140