1 #ifndef lint 2 static char vcid[] = "$Id: vecio.c,v 1.2 1995/08/17 23:42:48 curfman Exp curfman $"; 3 #endif 4 5 /* 6 This file contains simple binary read/write routines for matrices. 7 */ 8 9 #include "petsc.h" 10 #include <unistd.h> 11 #include "vec/vecimpl.h" 12 #include "sysio.h" 13 #include "matimpl.h" 14 #include "row.h" 15 16 /* -------------------------------------------------------------------- */ 17 /* This version reads from MATROW and writes to MATAIJ/MATROW 18 implementation. Eventually should not use generic MatSetValues, 19 but rather directly read data into appropriate location. Also, 20 should be able to read/write to/from any implementation. */ 21 22 /* @ 23 MatLoadBinary - Loads a matrix that has been stored in binary format 24 with MatViewBinary(). 25 26 Input Parameters: 27 . comm - MPI communicator 28 . fd - file descriptor (not FILE pointer). Use open() for this. 29 . outtype - type of output matrix 30 . ind - optional index set of matrix rows to be locally owned 31 (or 0 for loading the entire matrix on each processor) 32 . ind2 - optional index set with new matrix ordering (size = global 33 number of rows) 34 35 Output Parameters: 36 . newmat - new matrix 37 38 Notes: 39 In parallel, each processor can load a subset of rows (or the 40 entire matrix). This routine is especially useful when a large 41 matrix is stored on disk and only part of it is desired on each 42 processor. For example, a parallel solver may access only some of 43 the rows from each processor. The algorithm used here reads 44 relatively small blocks of data rather than reading the entire 45 matrix and then subsetting it. 46 47 Currently, the _entire_ matrix must be loaded. This should 48 probably change. 49 50 .seealso: MatViewBinary(), VecLoadBinary() 51 */ 52 int MatLoadBinary(MPI_Comm comm,int fd,MatType outtype,IS ind,IS ind2, 53 Mat *newmat) 54 { 55 Mat mat; 56 int rows, i, nz, nnztot, *rlen, ierr, lsize, gsize, *rptr, j, dstore; 57 int *cwork, rstart, rend, readst, *pind, *pind2, iinput, iglobal; 58 Scalar *awork; 59 MatType type; 60 long startloc, mark; 61 62 /* Get the location of the beginning of the matrix data, in case the 63 file contains multiple elements */ 64 startloc = lseek(fd,0L,SEEK_CUR); 65 MPIU_printf(comm,"startloc=%d\n",startloc); 66 type = MATROW; 67 if (outtype != MATROW && outtype != MATAIJ && outtype != MATMPIROW && 68 outtype != MATMPIAIJ) SETERRQ(1, 69 "MatLoadBinary: Only MATROW, MATAIJ, MATMPIROW, & MATAMPIAIJ supported."); 70 71 /* Read matrix header. Should this really be the full header? */ 72 ierr = SYRead(fd,(char *)&type,sizeof(int),SYINT); CHKERRQ(ierr); 73 if (type != MATROW) 74 SETERRQ(1,"MatLoadBinary: Only MATROW input currently supported"); 75 ierr = SYRead(fd,(char *)&rows,sizeof(int),0); CHKERRQ(ierr); 76 ierr = SYRead(fd,(char *)&nnztot,sizeof(int),0); CHKERRQ(ierr); 77 MPIU_printf(comm,"Input matrix: rows=%d, nnztot=%d\n",rows,nnztot); 78 79 /* Check sizes, form index set if necessary */ 80 if (!ind) 81 {ierr = ISCreateStrideSequential(comm,rows,0,1,&ind); CHKERRQ(ierr);} 82 ierr = ISGetLocalSize(ind,&lsize); CHKERRQ(ierr); 83 MPI_Allreduce(&lsize,&gsize,1,MPI_INT,MPI_SUM,comm); 84 if (gsize != rows) 85 SETERRQ(1,"MatLoadBinary: Incompatible parallel matrix size."); 86 ierr = ISGetIndices(ind,&pind); CHKERRQ(ierr); 87 ierr = ISGetIndices(ind,&pind2); CHKERRQ(ierr); 88 89 /* Allocate work space */ 90 rlen = (int *) PETSCMALLOC( rows * sizeof(int)); CHKPTRQ(rlen); 91 rptr = (int *) PETSCMALLOC( (rows+1) * sizeof(int)); CHKPTRQ(rptr); 92 cwork = (int *) PETSCMALLOC( rows*sizeof(int)); CHKPTRQ(cwork); 93 awork = (Scalar *) PETSCMALLOC( rows*sizeof(Scalar)); CHKPTRQ(awork); 94 95 /* Read row length info and form matrix memory allocation size */ 96 ierr = SYRead(fd,(char *)rlen,rows*sizeof(int),SYINT); CHKERRQ(ierr); 97 ierr = SYReadBuffer(-1,(long)0,0,(char*)0,SYINT); CHKERRQ(ierr); 98 99 /* This should be fixed */ 100 dstore = 5; 101 for ( i=0; i<lsize; i++ ) rptr[i] = PETSCMAX(rlen[pind[i]] - dstore,0); 102 103 /* Form new matrix */ 104 if (outtype == MATROW) 105 ierr = MatCreateSequentialRow(comm,rows,rows,0,rlen,&mat); 106 else if (outtype == MATAIJ) 107 ierr = MatCreateSequentialAIJ(comm,rows,rows,0,rlen,&mat); 108 else if (outtype == MATMPIROW) 109 ierr = MatCreateMPIRow(comm,lsize,PETSC_DECIDE,gsize,gsize,dstore, 110 0,0,rptr,&mat); 111 else if (outtype == MATMPIAIJ) 112 ierr = MatCreateMPIAIJ(comm,lsize,PETSC_DECIDE,gsize,gsize,dstore, 113 0,0,rptr,&mat); 114 CHKERRQ(ierr); 115 116 /* Form row pointers */ 117 rptr[0] = 0; 118 for (i=0; i<rows; i++) rptr[i+1] = rptr[i] + rlen[i]; 119 120 MatGetOwnershipRange(mat,&rstart,&rend); 121 mark = startloc + (rows+2)*sizeof(int); 122 for ( i=0; i<lsize; i++ ) { 123 iglobal = i + rstart; 124 iinput = pind[i]; 125 nz = rlen[iinput]; 126 readst = mark + rptr[iinput]*(sizeof(int)+sizeof(Scalar)); 127 ierr = SYReadBuffer(fd,readst,nz*sizeof(int),(char*)cwork,SYINT); 128 CHKERRQ(ierr); 129 ierr = SYReadBuffer(fd,readst+nz*sizeof(int),nz*sizeof(Scalar), 130 (char *)awork,SYSCALAR); CHKERRQ(ierr); 131 for (j=0; j<nz; j++) cwork[j] = pind2[cwork[j]]; 132 ierr = MatSetValues(mat,1,&iglobal,nz,cwork,awork,INSERTVALUES); 133 CHKERRQ(ierr); 134 } 135 PETSCFREE(rlen); PETSCFREE(rptr); PETSCFREE(cwork); PETSCFREE(awork); 136 ierr = MatAssemblyBegin(mat,FINAL_ASSEMBLY); CHKERRQ(ierr); 137 ierr = MatAssemblyEnd(mat,FINAL_ASSEMBLY); CHKERRQ(ierr); 138 ierr = ISRestoreIndices(ind,&pind); CHKERRQ(ierr); 139 ierr = ISRestoreIndices(ind,&pind2); CHKERRQ(ierr); 140 141 *newmat = mat; 142 return 0; 143 } 144 /* ------------------------------------------------------------- */ 145 146 int MatViewBinary(Mat mat,int fd) 147 { 148 Mat_Row *mrow = (Mat_Row *) mat->data; 149 MatiVec *vs; 150 int i, rows = mrow->m, nz, nnztot, ierr, *lwork; 151 152 if (mat->type != MATROW) SETERRQ(1,"Only MATROW currently supported."); 153 nnztot = 0; 154 for (i=0; i<rows; i++) nnztot += mrow->rs[i]->nz; 155 156 /* Write header */ 157 ierr = SYWrite(fd,(char *)&rows,sizeof(int),0,0); CHKERRQ(ierr); 158 ierr = SYWrite(fd,(char *)&nnztot,sizeof(int),0,0); CHKERRQ(ierr); 159 160 /* Write row length info */ 161 lwork = (int *)PETSCMALLOC( rows * sizeof(int)); CHKPTRQ(lwork); 162 for (i=0;i<rows;i++) { 163 lwork[i] = mrow->rs[i]->nz; 164 } 165 ierr = SYWrite(fd,(char *)lwork,rows*sizeof(int),0,0); CHKERRQ(ierr); 166 PETSCFREE(lwork); 167 168 /* Write matrix data by rows */ 169 for (i=0;i<rows;i++) { 170 vs = mrow->rs[i]; 171 nz = vs->nz; 172 ierr = SYWrite(fd,(char *)vs->i,nz*sizeof(int),0,0); CHKERRQ(ierr); 173 ierr = SYWrite(fd,(char *)vs->v,nz*sizeof(Scalar),0,0); CHKERRQ(ierr); 174 } 175 return 0; 176 } 177