2: #define USE_FAST_MAT_SET_VALUES

 4:  #include petsc.h
 5:  #include petscviewer.h

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


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

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

 40:   MPI_Comm_size(comm,&size);
 41:   MPI_Comm_rank(comm,&rank);

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

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

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

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

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

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

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

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

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

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

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

120:   MatSetOption(A,MAT_COLUMNS_SORTED);
121:   for (i=0; i<m; i++) {
122:     ourlens[i] += offlens[i];
123:   }
124:   PetscFree(offlens);

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

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

141:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
142:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
143:   *newmat = A;
144:   PetscViewerDestroy(viewer1);
145:   PetscViewerDestroy(viewer2);
146:   return(0);
147: }

149: int main(int argc,char **args)
150: {
152:   Mat            A;
153:   char           name[1024];
154:   PetscTruth     flg;

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