2: static char help[] = "Each process opens the file and reads its part. Not scalable do not copy\n";

  4:  #include <petscsys.h>
  5:  #include <petscviewer.h>
  6:  #include <petscmat.h>

  8: /*
  9:    Opens a separate file for each process and reads in ITS portion
 10:   of a large parallel matrix. Only requires enough memory to store
 11:   the processes portion of the matrix ONCE.

 13:     petsc-maint@mcs.anl.gov
 14: */
 15: int Mat_Parallel_Load(MPI_Comm comm,const char *name,Mat *newmat)
 16: {
 17:   Mat            A;
 18:   PetscScalar    *vals;
 20:   PetscMPIInt    rank,size;
 21:   PetscInt       i,j,rstart,rend;
 22:   PetscInt       header[4],M,N,m;
 23:   PetscInt       *ourlens,*offlens,jj,*mycols,maxnz;
 24:   PetscInt       cend,cstart,n,*rowners;
 25:   int            fd1,fd2;
 26:   PetscViewer    viewer1,viewer2;

 29:   MPI_Comm_size(comm,&size);
 30:   MPI_Comm_rank(comm,&rank);

 32:   /* Open the files; each process opens its own file */
 33:   PetscViewerBinaryOpen(PETSC_COMM_SELF,name,FILE_MODE_READ,&viewer1);
 34:   PetscViewerBinaryGetDescriptor(viewer1,&fd1);
 35:   PetscBinaryRead(fd1,(char*)header,4,PETSC_INT);

 37:   /* open the file twice so that later we can read entries from two different parts of the
 38:      file at the same time. Note that due to file caching this should not impact performance */
 39:   PetscViewerBinaryOpen(PETSC_COMM_SELF,name,FILE_MODE_READ,&viewer2);
 40:   PetscViewerBinaryGetDescriptor(viewer2,&fd2);
 41:   PetscBinaryRead(fd2,(char*)header,4,PETSC_INT);

 43:   /* error checking on files */
 44:   if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object");
 45:   MPI_Allreduce(header+2,&N,1,MPIU_INT,MPI_SUM,comm);
 46:   if (N != size*header[2]) SETERRQ(PETSC_COMM_SELF,1,"All files must have matrices with the same number of total columns");

 48:   /* number of rows in matrix is sum of rows in all files */
 49:   m    = header[1]; N = header[2];
 50:   MPI_Allreduce(&m,&M,1,MPIU_INT,MPI_SUM,comm);

 52:   /* determine rows of matrices owned by each process */
 53:   PetscMalloc1(size+1,&rowners);
 54:   MPI_Allgather(&m,1,MPIU_INT,rowners+1,1,MPIU_INT,comm);
 55:   rowners[0] = 0;
 56:   for (i=2; i<=size; i++) {
 57:     rowners[i] += rowners[i-1];
 58:   }
 59:   rstart = rowners[rank];
 60:   rend   = rowners[rank+1];
 61:   PetscFree(rowners);

 63:   /* determine column ownership if matrix is not square */
 64:   if (N != M) {
 65:     n      = N/size + ((N % size) > rank);
 66:     MPI_Scan(&n,&cend,1,MPIU_INT,MPI_SUM,comm);
 67:     cstart = cend - n;
 68:   } else {
 69:     cstart = rstart;
 70:     cend   = rend;
 71:     n      = cend - cstart;
 72:   }

 74:   /* read in local row lengths */
 75:   PetscMalloc1(m,&ourlens);
 76:   PetscMalloc1(m,&offlens);
 77:   PetscBinaryRead(fd1,ourlens,m,PETSC_INT);
 78:   PetscBinaryRead(fd2,ourlens,m,PETSC_INT);

 80:   /* determine buffer space needed for column indices of any one row*/
 81:   maxnz = 0;
 82:   for (i=0; i<m; i++) {
 83:     maxnz = PetscMax(maxnz,ourlens[i]);
 84:   }

 86:   /* allocate enough memory to hold a single row of column indices */
 87:   PetscMalloc1(maxnz,&mycols);

 89:   /* loop over local rows, determining number of off diagonal entries */
 90:   PetscMemzero(offlens,m*sizeof(PetscInt));
 91:   for (i=0; i<m; i++) {
 92:     PetscBinaryRead(fd1,mycols,ourlens[i],PETSC_INT);
 93:     for (j=0; j<ourlens[i]; j++) {
 94:       if (mycols[j] < cstart || mycols[j] >= cend) offlens[i]++;
 95:     }
 96:   }

 98:   /* on diagonal entries are all that were not counted as off-diagonal */
 99:   for (i=0; i<m; i++) {
100:     ourlens[i] -= offlens[i];
101:   }

103:   /* create our matrix */
104:   MatCreate(comm,&A);
105:   MatSetSizes(A,m,n,M,N);
106:   MatSetType(A,MATMPIAIJ);
107:   MatMPIAIJSetPreallocation(A,0,ourlens,0,offlens);

109:   for (i=0; i<m; i++) {
110:     ourlens[i] += offlens[i];
111:   }
112:   PetscFree(offlens);

114:   /* allocate enough memory to hold a single row of matrix values */
115:   PetscMalloc1(maxnz,&vals);

117:   /* read in my part of the matrix numerical values and columns 1 row at a time and put in matrix  */
118:   jj = rstart;
119:   for (i=0; i<m; i++) {
120:     PetscBinaryRead(fd1,vals,ourlens[i],PETSC_SCALAR);
121:     PetscBinaryRead(fd2,mycols,ourlens[i],PETSC_INT);
122:     MatSetValues(A,1,&jj,ourlens[i],mycols,vals,INSERT_VALUES);
123:     jj++;
124:   }
125:   PetscFree(ourlens);
126:   PetscFree(vals);
127:   PetscFree(mycols);

129:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
130:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
131:   *newmat = A;
132:   PetscViewerDestroy(&viewer1);
133:   PetscViewerDestroy(&viewer2);
134:   return(0);
135: }

137: int main(int argc,char **args)
138: {
140:   Mat            A;
141:   char           name[1024];
142:   PetscBool      flg;

144:   PetscInitialize(&argc,&args,0,help);if (ierr) return ierr;
145:   PetscOptionsGetString(NULL,NULL,"-f",name,1024,&flg);
146:   if (!flg) SETERRQ(PETSC_COMM_SELF,1,"Must pass in filename with -f option");
147:   Mat_Parallel_Load(PETSC_COMM_WORLD,name,&A);
148:   MatDestroy(&A);
149:   PetscFinalize();
150:   return ierr;
151: }