Actual source code: ex5.c

      2: #define USE_FAST_MAT_SET_VALUES

     4:  #include petsc.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 seperate 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,PETSC_FILE_RDONLY,&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,PETSC_FILE_RDONLY,&viewer2);
     49:   PetscViewerBinaryGetDescriptor(viewer2,&fd2);
     50:   PetscBinaryRead(fd2,(char *)header,4,PETSC_INT);

     52:   /* error checking on files */
     53:   if (header[0] != MAT_FILE_COOKIE) SETERRQ(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(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:   MatSetOption(A,MAT_COLUMNS_SORTED);
    119:   for (i=0; i<m; i++) {
    120:     ourlens[i] += offlens[i];
    121:   }
    122:   PetscFree(offlens);

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

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

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

    147: int main(int argc,char **args)
    148: {
    150:   Mat            A;
    151:   char           name[1024];
    152:   PetscTruth     flg;

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