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