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