2: static char help[] = "Reads a PETSc matrix from a file partitions it\n\n";

  4: /*T
  5:    Concepts: partitioning
  6:    Processors: n
  7: T*/

  9: /* 
 10:   Include "petscmat.h" so that we can use matrices.  Note that this file
 11:   automatically includes:
 12:      petsc.h       - base PETSc routines   petscvec.h - vectors
 13:      petscsys.h    - system routines       petscmat.h - matrices
 14:      petscis.h     - index sets            
 15:      petscviewer.h - viewers    

 17:   Example of usage:  
 18:     mpirun -np 3 ex73 -f <matfile> -mat_partitioning_type parmetis/scotch -viewer_binary_skip_info -nox
 19: */
 20:  #include petscksp.h

 24: int main(int argc,char **args)
 25: {
 26:   MatType         mtype = MATMPIAIJ; /* matrix format */
 27:   Mat             A,B;               /* matrix */
 28:   PetscViewer     fd;                /* viewer */
 29:   char            file[PETSC_MAX_PATH_LEN];         /* input file name */
 30:   PetscTruth      flg;
 31:   PetscInt        ierr,*nlocal;
 32:   PetscMPIInt     rank,size;
 33:   MatPartitioning part;
 34:   IS              is,isn;

 36:   PetscInitialize(&argc,&args,(char *)0,help);
 37:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 38:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 40:   /* 
 41:      Determine file from which we read the matrix
 42:   */
 43:   PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN-1,&flg);

 45:   /* 
 46:        Open binary file.  Note that we use FILE_MODE_READ to indicate
 47:        reading from this file.
 48:   */
 49:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);

 51:   /*
 52:       Load the matrix and vector; then destroy the viewer.
 53:   */
 54:   MatLoad(fd,mtype,&A);
 55:   PetscViewerDestroy(fd);

 57:   MatView(A,PETSC_VIEWER_DRAW_WORLD);

 59:   /*
 60:        Partition the graph of the matrix 
 61:   */
 62:   MatPartitioningCreate(PETSC_COMM_WORLD,&part);
 63:   MatPartitioningSetAdjacency(part,A);
 64:   MatPartitioningSetFromOptions(part);
 65:   /* get new processor owner number of each vertex */
 66:   MatPartitioningApply(part,&is);
 67:   /* get new global number of each old global number */
 68:   ISPartitioningToNumbering(is,&isn);
 69:   PetscMalloc(size*sizeof(PetscInt),&nlocal);
 70:   /* get number of new vertices for each processor */
 71:   ISPartitioningCount(is,nlocal);
 72:   ISDestroy(is);

 74:   /* get old global number of each new global number */
 75:   ISInvertPermutation(isn,nlocal[rank],&is);
 76:   PetscFree(nlocal);
 77:   ISDestroy(isn);
 78:   MatPartitioningDestroy(part);

 80:   ISSort(is);
 81:   ISAllGather(is,&isn);


 84:   MatGetSubMatrix(A,is,isn,PETSC_DECIDE,MAT_INITIAL_MATRIX,&B);
 85:   ISDestroy(is);
 86:   ISDestroy(isn);

 88:   MatView(B,PETSC_VIEWER_DRAW_WORLD);

 90:   /*
 91:        Free work space.  All PETSc objects should be destroyed when they
 92:        are no longer needed.
 93:   */
 94:   MatDestroy(A);
 95:   MatDestroy(B);


 98:   PetscFinalize();
 99:   return 0;
100: }