1 static char help[] = "Tests MatPermute() for a square matrix in parallel.\n\n"; 2 #include <petscmat.h> 3 /* Results: 4 Sequential: 5 - seqaij: correct permutation 6 - seqbaij: permutation not supported for this MATTYPE 7 - seqsbaij: permutation not supported for this MATTYPE 8 Parallel: 9 - mpiaij: incorrect permutation (disable this op for now) 10 - seqbaij: correct permutation, but all manner of memory leaks 11 (disable this op for now, because it appears broken for nonsquare matrices; see ex151) 12 - seqsbaij: permutation not supported for this MATTYPE 13 14 */ 15 int main(int argc, char **argv) 16 { 17 const struct { 18 PetscInt i, j; 19 PetscScalar v; 20 } entries[] = { 21 {0, 3, 1.}, 22 {1, 2, 2.}, 23 {2, 1, 3.}, 24 {2, 4, 4.}, 25 {3, 0, 5.}, 26 {3, 3, 6.}, 27 {4, 1, 7.}, 28 {4, 4, 8.} 29 }; 30 const PetscInt ixrow[5] = {4, 2, 1, 3, 0}, ixcol[5] = {3, 2, 1, 4, 0}; 31 Mat A, B; 32 PetscInt i, rstart, rend, cstart, cend; 33 IS isrow, iscol; 34 PetscViewer viewer, sviewer; 35 PetscBool view_sparse; 36 37 PetscFunctionBeginUser; 38 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 39 /* ------- Assemble matrix, --------- */ 40 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 41 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 5, 5)); 42 PetscCall(MatSetFromOptions(A)); 43 PetscCall(MatSetUp(A)); 44 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 45 PetscCall(MatGetOwnershipRangeColumn(A, &cstart, &cend)); 46 47 for (i = 0; i < (PetscInt)PETSC_STATIC_ARRAY_LENGTH(entries); i++) PetscCall(MatSetValue(A, entries[i].i, entries[i].j, entries[i].v, INSERT_VALUES)); 48 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 49 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 50 51 /* ------ Prepare index sets ------ */ 52 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, rend - rstart, ixrow + rstart, PETSC_USE_POINTER, &isrow)); 53 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 5, ixcol, PETSC_USE_POINTER, &iscol)); 54 PetscCall(ISSetPermutation(isrow)); 55 PetscCall(ISSetPermutation(iscol)); 56 57 PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD, &viewer)); 58 view_sparse = PETSC_FALSE; 59 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_sparse", &view_sparse, NULL)); 60 if (!view_sparse) PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_DENSE)); 61 PetscCall(PetscViewerASCIIPrintf(viewer, "Original matrix\n")); 62 PetscCall(MatView(A, viewer)); 63 64 PetscCall(MatPermute(A, isrow, iscol, &B)); 65 PetscCall(PetscViewerASCIIPrintf(viewer, "Permuted matrix\n")); 66 PetscCall(MatView(B, viewer)); 67 if (!view_sparse) PetscCall(PetscViewerPopFormat(viewer)); 68 69 PetscCall(PetscViewerASCIIPrintf(viewer, "Row permutation\n")); 70 PetscCall(ISView(isrow, viewer)); 71 PetscCall(PetscViewerASCIIPrintf(viewer, "Column permutation\n")); 72 PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer)); 73 PetscCall(ISView(iscol, sviewer)); 74 PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer)); 75 76 /* Free data structures */ 77 PetscCall(ISDestroy(&isrow)); 78 PetscCall(ISDestroy(&iscol)); 79 PetscCall(MatDestroy(&A)); 80 PetscCall(MatDestroy(&B)); 81 82 PetscCall(PetscFinalize()); 83 return 0; 84 } 85