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