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