Actual source code: ex5.c
petsc-3.10.5 2019-03-28
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: }
154: /*TEST
156: test:
157: nsize: 2
158: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
159: args: -f ${DATAFILESPATH}/matrices/medium -viewer_binary_skip_info
160: TODO: Need to develop comparison test
162: TEST*/