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