Actual source code: ex166.c

  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

 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;

 38:   PetscInitialize(&argc, &argv, (char *)0, help);
 39:   /* ------- Assemble matrix, --------- */
 40:   MatCreate(PETSC_COMM_WORLD, &A);
 41:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 5, 5);
 42:   MatSetFromOptions(A);
 43:   MatSetUp(A);
 44:   MatGetOwnershipRange(A, &rstart, &rend);
 45:   MatGetOwnershipRangeColumn(A, &cstart, &cend);

 47:   for (i = 0; i < (PetscInt)PETSC_STATIC_ARRAY_LENGTH(entries); i++) MatSetValue(A, entries[i].i, entries[i].j, entries[i].v, INSERT_VALUES);
 48:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 49:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

 51:   /* ------ Prepare index sets ------ */
 52:   ISCreateGeneral(PETSC_COMM_WORLD, rend - rstart, ixrow + rstart, PETSC_USE_POINTER, &isrow);
 53:   ISCreateGeneral(PETSC_COMM_SELF, 5, ixcol, PETSC_USE_POINTER, &iscol);
 54:   ISSetPermutation(isrow);
 55:   ISSetPermutation(iscol);

 57:   PetscViewerASCIIGetStdout(PETSC_COMM_WORLD, &viewer);
 58:   view_sparse = PETSC_FALSE;
 59:   PetscOptionsGetBool(NULL, NULL, "-view_sparse", &view_sparse, NULL);
 60:   if (!view_sparse) PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_DENSE);
 61:   PetscViewerASCIIPrintf(viewer, "Original matrix\n");
 62:   MatView(A, viewer);

 64:   MatPermute(A, isrow, iscol, &B);
 65:   PetscViewerASCIIPrintf(viewer, "Permuted matrix\n");
 66:   MatView(B, viewer);
 67:   if (!view_sparse) PetscViewerPopFormat(viewer);

 69:   PetscViewerASCIIPrintf(viewer, "Row permutation\n");
 70:   ISView(isrow, viewer);
 71:   PetscViewerASCIIPrintf(viewer, "Column permutation\n");
 72:   PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
 73:   ISView(iscol, sviewer);
 74:   PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);

 76:   /* Free data structures */
 77:   ISDestroy(&isrow);
 78:   ISDestroy(&iscol);
 79:   MatDestroy(&A);
 80:   MatDestroy(&B);

 82:   PetscFinalize();
 83:   return 0;
 84: }