Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 21667) +++ /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 21668) @@ -683,4 +683,5 @@ case 7: case 8: + case 10: /*We are not interested in this record, too memory intensive. Skip it: */ /*skip: */ @@ -825,4 +826,5 @@ case 7: break; //do nothing. not interested in this type of data, which is memory intensive. case 8: break; //do nothing. not interested in this type of data, which is memory intensive. + case 10: break; //do nothing. not interested in this type of data, which is memory intensive. case 9: ISSM_MPI_Bcast(&numstrings,1,ISSM_MPI_INT,0,IssmComm::GetComm()); @@ -1132,5 +1134,5 @@ /*Set file pointer to beginning of the data: */ fid=this->SetFilePointerToData(&code,NULL,data_name); - if(code!=5 && code!=6 && code!=7)_error_("expecting a IssmDouble, integer or boolean matrix for \""<(M*N); - - /*Read matrix on node 0, then broadcast: */ - if(my_rank==0){ - if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix "); - } - ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); - - *pmatrix=xNew(M*N); - for(int i=0;i(matrix); - } - else - *pmatrix=NULL; + if(code==10){ + /*Special case for Compressed mat*/ + IssmPDouble offset,range; + if(my_rank==0) if(fread(&offset,sizeof(IssmPDouble),1,fid)!=1) _error_("could not read offset"); + ISSM_MPI_Bcast(&offset,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); + + if(my_rank==0) if(fread(&range,sizeof(IssmPDouble),1,fid)!=1) _error_("could not read range"); + ISSM_MPI_Bcast(&range,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); + + *pmatrix=xNew(M*N); + + /*Read matrix*/ + uint8_t* rawmatrix=xNew((M-1)*N); + if(my_rank==0) if(fread(rawmatrix,(M-1)*N*sizeof(char),1,fid)!=1) _error_("could not read matrix "); + ISSM_MPI_Bcast(rawmatrix,(M-1)*N,ISSM_MPI_CHAR,0,IssmComm::GetComm()); + + for(int i=0;i<(M-1)*N;++i) (*pmatrix)[i]=offset+range*reCast(rawmatrix[i])/255.; + xDelete(rawmatrix); + + /*read time now*/ + IssmPDouble* timematrix=xNew(N); + if(my_rank==0) if(fread(timematrix,N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read time in compressed matrix"); + ISSM_MPI_Bcast(timematrix,N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); + + for(int i=0;i(timematrix); + + } + else{ + /*Read matrix on node 0, then broadcast: */ + matrix=xNew(M*N); + if(my_rank==0) if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix "); + ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); + + *pmatrix=xNew(M*N); + for(int i=0;i(matrix); + } + } + else{ + *pmatrix=NULL; + } /*Assign output pointers: */ if(pM) *pM=M; @@ -1542,5 +1572,5 @@ /*Defaulting only supported for double arrays*/ - if(code!=7) _error_(vector_name<<" is not a double array"); + if(code!=7 && code!=10) _error_(vector_name<<" is not a double array"); this->FetchData(&doublearray,&M,&N,vector_name); @@ -2466,5 +2496,5 @@ /*if record_code points to a vector, get its type (nodal or elementary): */ - if(5<=record_code && record_code<=7){ + if((5<=record_code && record_code<=7) || record_code==10){ if(fread(&vector_type,sizeof(int),1,fid)!=1) _error_("Could not read vector_type"); } Index: /issm/trunk-jpl/src/m/classes/SMBforcing.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBforcing.m (revision 21667) +++ /issm/trunk-jpl/src/m/classes/SMBforcing.m (revision 21668) @@ -69,4 +69,5 @@ WriteData(fid,prefix,'name','md.smb.model','data',1,'format','Integer'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','mass_balance','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + %WriteData(fid,prefix,'object',self,'class','smb','fieldname','mass_balance','format','CompressedMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); %process requested outputs Index: /issm/trunk-jpl/src/m/solve/WriteData.m =================================================================== --- /issm/trunk-jpl/src/m/solve/WriteData.m (revision 21667) +++ /issm/trunk-jpl/src/m/solve/WriteData.m (revision 21668) @@ -168,4 +168,41 @@ fwrite(fid,data','double'); %get to the "c" convention, hence the transpose end + % }}} +elseif strcmpi(format,'CompressedMat'), % {{{ + + %Get size + s=size(data); + + if s(1)==1 | s(2)==1, + %No need to use Compressed format + error('Not needed (should call WriteData with DoubleMat)'); + end + + %first write length of record + recordlength=4+4+8+8+1*(s(1)-1)*s(2)+8*s(2)+4+4; %2 integers (32 bits) + the matrix + code + matrix type + if recordlength>2^31; error(['field ' name ' cannot be marshalled because it is larger than 2^31 bytes!']); end + fwrite(fid,recordlength,'int'); + + %write data code and matrix type: + fwrite(fid,FormatToCode(format),'int'); + fwrite(fid,mattype,'int'); + + %write matrix size + fwrite(fid,s(1),'int'); + fwrite(fid,s(2),'int'); + + %Write offset and range + A = data(1:end-1,:); + offset = min(A(:)); + range = max(A(:)) - offset; + fwrite(fid,offset,'double'); + fwrite(fid,range,'double'); + + %Convert data to uint8 and write it + A=uint8((A-offset)/range*255); + fwrite(fid,A','uint8'); %get to the "c" convention, hence the transpose + + %Write last row as double (time) + fwrite(fid,data(end,:),'double'); % }}} elseif strcmpi(format,'MatArray'), % {{{ @@ -258,4 +295,6 @@ elseif strcmpi(format,'StringArray'), code=9; + elseif strcmpi(format,'CompressedMat'), + code=10; else error('FormatToCode error message: data type not supported yet!');