2: #define USE_FAST_MAT_SET_VALUES

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

  7: #if defined(USE_FAST_MAT_SET_VALUES)
  8: #include <../src/mat/impls/aij/mpi/mpiaij.h>
  9: #define MatSetValues MatSetValues_MPIAIJ
 10: #else 
 11: #include <petscmat.h>
 12: #endif


 15: /*
 16:    Opens a separate file for each process and reads in ITS portion
 17:   of a large parallel matrix. Only requires enough memory to store
 18:   the processes portion of the matrix ONCE.

 20:     petsc-maint@mcs.anl.gov
 21: */
 24: int Mat_Parallel_Load(MPI_Comm comm,const char *name,Mat *newmat)
 25: {
 26:   Mat            A;
 27:   PetscScalar    *vals;
 29:   PetscMPIInt    rank,size;
 30:   PetscInt       i,j,rstart,rend;
 31:   PetscInt       header[4],M,N,m;
 32:   PetscInt       *ourlens,*offlens,jj,*mycols,maxnz;
 33:   PetscInt       cend,cstart,n,*rowners;
 34:   int            fd1,fd2;
 35:   PetscViewer    viewer1,viewer2;

 38:   MPI_Comm_size(comm,&size);
 39:   MPI_Comm_rank(comm,&rank);

 41:   /* Open the files; each process opens its own file */
 42:   PetscViewerBinaryOpen(PETSC_COMM_SELF,name,FILE_MODE_READ,&viewer1);
 43:   PetscViewerBinaryGetDescriptor(viewer1,&fd1);
 44:   PetscBinaryRead(fd1,(char *)header,4,PETSC_INT);

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

 52:   /* error checking on files */
 53:   if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object");
 54:   MPI_Allreduce(header+2,&N,1,MPI_INT,MPI_SUM,comm);
 55:   if (N != size*header[2]) SETERRQ(PETSC_COMM_SELF,1,"All files must have matrices with the same number of total columns");
 56: 
 57:   /* number of rows in matrix is sum of rows in all files */
 58:   m = header[1]; N = header[2];
 59:   MPI_Allreduce(&m,&M,1,MPI_INT,MPI_SUM,comm);

 61:   /* determine rows of matrices owned by each process */
 62:   PetscMalloc((size+1)*sizeof(PetscInt),&rowners);
 63:   MPI_Allgather(&m,1,MPIU_INT,rowners+1,1,MPIU_INT,comm);
 64:   rowners[0] = 0;
 65:   for (i=2; i<=size; i++) {
 66:     rowners[i] += rowners[i-1];
 67:   }
 68:   rstart = rowners[rank];
 69:   rend   = rowners[rank+1];
 70:   PetscFree(rowners);

 72:   /* determine column ownership if matrix is not square */
 73:   if (N != M) {
 74:     n      = N/size + ((N % size) > rank);
 75:     MPI_Scan(&n,&cend,1,MPIU_INT,MPI_SUM,comm);
 76:     cstart = cend - n;
 77:   } else {
 78:     cstart = rstart;
 79:     cend   = rend;
 80:     n      = cend - cstart;
 81:   }

 83:   /* read in local row lengths */
 84:   PetscMalloc(m*sizeof(PetscInt),&ourlens);
 85:   PetscMalloc(m*sizeof(PetscInt),&offlens);
 86:   PetscBinaryRead(fd1,ourlens,m,PETSC_INT);
 87:   PetscBinaryRead(fd2,ourlens,m,PETSC_INT);

 89:   /* determine buffer space needed for column indices of any one row*/
 90:   maxnz = 0;
 91:   for (i=0; i<m; i++) {
 92:     maxnz = PetscMax(maxnz,ourlens[i]);
 93:   }

 95:   /* allocate enough memory to hold a single row of column indices */
 96:   PetscMalloc(maxnz*sizeof(PetscInt),&mycols);

 98:   /* loop over local rows, determining number of off diagonal entries */
 99:   PetscMemzero(offlens,m*sizeof(PetscInt));
100:   for (i=0; i<m; i++) {
101:     PetscBinaryRead(fd1,mycols,ourlens[i],PETSC_INT);
102:     for (j=0; j<ourlens[i]; j++) {
103:       if (mycols[j] < cstart || mycols[j] >= cend) offlens[i]++;
104:     }
105:   }

107:   /* on diagonal entries are all that were not counted as off-diagonal */
108:   for (i=0; i<m; i++) {
109:     ourlens[i] -= offlens[i];
110:   }

112:   /* create our matrix */
113:   MatCreate(comm,&A);
114:   MatSetSizes(A,m,n,M,N);
115:   MatSetType(A,MATMPIAIJ);
116:   MatMPIAIJSetPreallocation(A,0,ourlens,0,offlens);

118:   for (i=0; i<m; i++) {
119:     ourlens[i] += offlens[i];
120:   }
121:   PetscFree(offlens);

123:   /* allocate enough memory to hold a single row of matrix values */
124:   PetscMalloc(maxnz*sizeof(PetscScalar),&vals);

126:   /* read in my part of the matrix numerical values and columns 1 row at a time and put in matrix  */
127:   jj = rstart;
128:   for (i=0; i<m; i++) {
129:     PetscBinaryRead(fd1,vals,ourlens[i],PETSC_SCALAR);
130:     PetscBinaryRead(fd2,mycols,ourlens[i],PETSC_INT);
131:     MatSetValues(A,1,&jj,ourlens[i],mycols,vals,INSERT_VALUES);
132:     jj++;
133:   }
134:   PetscFree(ourlens);
135:   PetscFree(vals);
136:   PetscFree(mycols);

138:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
139:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
140:   *newmat = A;
141:   PetscViewerDestroy(&viewer1);
142:   PetscViewerDestroy(&viewer2);
143:   return(0);
144: }

146: int main(int argc,char **args)
147: {
149:   Mat            A;
150:   char           name[1024];
151:   PetscBool      flg;

153:   PetscInitialize(&argc,&args,0,0);
154:   PetscOptionsGetString(PETSC_NULL,"-f",name,1024,&flg);
155:   if (!flg) SETERRQ(PETSC_COMM_SELF,1,"Must pass in filename with -f option");
156:   Mat_Parallel_Load(PETSC_COMM_WORLD,name,&A);
157:   MatDestroy(&A);
158:   PetscFinalize();
159:   return 0;
160: }