Actual source code: ex5.c
petsc-3.7.7 2017-09-25
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: */
17: int Mat_Parallel_Load(MPI_Comm comm,const char *name,Mat *newmat)
18: {
19: Mat A;
20: PetscScalar *vals;
22: PetscMPIInt rank,size;
23: PetscInt i,j,rstart,rend;
24: PetscInt header[4],M,N,m;
25: PetscInt *ourlens,*offlens,jj,*mycols,maxnz;
26: PetscInt cend,cstart,n,*rowners;
27: int fd1,fd2;
28: PetscViewer viewer1,viewer2;
31: MPI_Comm_size(comm,&size);
32: MPI_Comm_rank(comm,&rank);
34: /* Open the files; each process opens its own file */
35: PetscViewerBinaryOpen(PETSC_COMM_SELF,name,FILE_MODE_READ,&viewer1);
36: PetscViewerBinaryGetDescriptor(viewer1,&fd1);
37: PetscBinaryRead(fd1,(char*)header,4,PETSC_INT);
39: /* open the file twice so that later we can read entries from two different parts of the
40: file at the same time. Note that due to file caching this should not impact performance */
41: PetscViewerBinaryOpen(PETSC_COMM_SELF,name,FILE_MODE_READ,&viewer2);
42: PetscViewerBinaryGetDescriptor(viewer2,&fd2);
43: PetscBinaryRead(fd2,(char*)header,4,PETSC_INT);
45: /* error checking on files */
46: if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object");
47: MPI_Allreduce(header+2,&N,1,MPIU_INT,MPI_SUM,comm);
48: if (N != size*header[2]) SETERRQ(PETSC_COMM_SELF,1,"All files must have matrices with the same number of total columns");
50: /* number of rows in matrix is sum of rows in all files */
51: m = header[1]; N = header[2];
52: MPI_Allreduce(&m,&M,1,MPIU_INT,MPI_SUM,comm);
54: /* determine rows of matrices owned by each process */
55: PetscMalloc1(size+1,&rowners);
56: MPI_Allgather(&m,1,MPIU_INT,rowners+1,1,MPIU_INT,comm);
57: rowners[0] = 0;
58: for (i=2; i<=size; i++) {
59: rowners[i] += rowners[i-1];
60: }
61: rstart = rowners[rank];
62: rend = rowners[rank+1];
63: PetscFree(rowners);
65: /* determine column ownership if matrix is not square */
66: if (N != M) {
67: n = N/size + ((N % size) > rank);
68: MPI_Scan(&n,&cend,1,MPIU_INT,MPI_SUM,comm);
69: cstart = cend - n;
70: } else {
71: cstart = rstart;
72: cend = rend;
73: n = cend - cstart;
74: }
76: /* read in local row lengths */
77: PetscMalloc1(m,&ourlens);
78: PetscMalloc1(m,&offlens);
79: PetscBinaryRead(fd1,ourlens,m,PETSC_INT);
80: PetscBinaryRead(fd2,ourlens,m,PETSC_INT);
82: /* determine buffer space needed for column indices of any one row*/
83: maxnz = 0;
84: for (i=0; i<m; i++) {
85: maxnz = PetscMax(maxnz,ourlens[i]);
86: }
88: /* allocate enough memory to hold a single row of column indices */
89: PetscMalloc1(maxnz,&mycols);
91: /* loop over local rows, determining number of off diagonal entries */
92: PetscMemzero(offlens,m*sizeof(PetscInt));
93: for (i=0; i<m; i++) {
94: PetscBinaryRead(fd1,mycols,ourlens[i],PETSC_INT);
95: for (j=0; j<ourlens[i]; j++) {
96: if (mycols[j] < cstart || mycols[j] >= cend) offlens[i]++;
97: }
98: }
100: /* on diagonal entries are all that were not counted as off-diagonal */
101: for (i=0; i<m; i++) {
102: ourlens[i] -= offlens[i];
103: }
105: /* create our matrix */
106: MatCreate(comm,&A);
107: MatSetSizes(A,m,n,M,N);
108: MatSetType(A,MATMPIAIJ);
109: MatMPIAIJSetPreallocation(A,0,ourlens,0,offlens);
111: for (i=0; i<m; i++) {
112: ourlens[i] += offlens[i];
113: }
114: PetscFree(offlens);
116: /* allocate enough memory to hold a single row of matrix values */
117: PetscMalloc1(maxnz,&vals);
119: /* read in my part of the matrix numerical values and columns 1 row at a time and put in matrix */
120: jj = rstart;
121: for (i=0; i<m; i++) {
122: PetscBinaryRead(fd1,vals,ourlens[i],PETSC_SCALAR);
123: PetscBinaryRead(fd2,mycols,ourlens[i],PETSC_INT);
124: MatSetValues(A,1,&jj,ourlens[i],mycols,vals,INSERT_VALUES);
125: jj++;
126: }
127: PetscFree(ourlens);
128: PetscFree(vals);
129: PetscFree(mycols);
131: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
132: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
133: *newmat = A;
134: PetscViewerDestroy(&viewer1);
135: PetscViewerDestroy(&viewer2);
136: return(0);
137: }
139: int main(int argc,char **args)
140: {
142: Mat A;
143: char name[1024];
144: PetscBool flg;
146: PetscInitialize(&argc,&args,0,help);
147: PetscOptionsGetString(NULL,NULL,"-f",name,1024,&flg);
148: if (!flg) SETERRQ(PETSC_COMM_SELF,1,"Must pass in filename with -f option");
149: Mat_Parallel_Load(PETSC_COMM_WORLD,name,&A);
150: MatDestroy(&A);
151: PetscFinalize();
152: return 0;
153: }