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