xref: /petsc/share/petsc/matlab/PetscBinaryRead.m (revision 09b68a49ed2854d1e4985cc2aa6af33c7c4e69b3)
1function [varargout] = PetscBinaryRead(inarg,varargin)
2%
3%   [varargout] = PetscBinaryRead(inarg,['complex',false or true],['indices','int32' or 'int64'],['cell',cnt],['precision','float64' or 'float32'])
4%
5%  Reads in PETSc binary file matrices or vectors
6%  emits as MATLAB sparse matrice or vectors.
7%
8%  [] indices optional arguments
9%  There are no [] in the arguments
10%
11%  Examples: A = PetscBinaryRead('myfile'); read from file
12%            b = PetscBinaryRead(1024);   read from socket
13%            c = PetscBinaryRead();       read from default socket PETSc uses
14%
15%  Argument may be file name (string), socket number (integer)
16%  or any MATLAB class that provides the read() and close() methods
17%  [We provide PetscOpenFile() and PetscOpenSocket() for binary files and sockets]
18%  For example: fd = PetscOpenFile('filename');
19%                a = PetscBinaryRead(fd);
20%                b = PetscBinaryRead(fd);
21%
22%  'complex', true indicates the numbers in the file are complex, that is PETSc was built with --with-scalar-type=complex
23%  'indices','int64' indicates the PETSc program was built with --with-64-bit-indices
24%  'cell',cnt  means return a MATLAB cell array containing the first cnt objects in the file, use 10,000 to read in all objects
25%  'precision','float32' indicates the PETSc program was built with --with-precision=single
26%
27%  Examples:  A = PetscBinaryRead('myfile','cell',10000);  read all objects in file
28%             A = PetscBinaryRead(1024,'cell',2);  read two objects from socket
29%
30if nargin == 0
31  fd = PetscOpenSocket();
32else if ischar(inarg)
33  fd = PetscOpenFile(inarg);
34else if isnumeric(inarg)
35  fd = PetscOpenSocket(inarg);
36else % assume it is a PetscOpenFile or PetscOpenSocket object and handles read()
37  fd = inarg;
38end
39end
40end
41
42indices = 'int32';
43precision = 'float64';
44arecell = 0;
45arecomplex = false;
46
47tnargin = nargin;
48for l=1:nargin-2
49  if ischar(varargin{l}) && strcmpi(varargin{l},'indices')
50    tnargin = min(l,tnargin-1);
51    indices = varargin{l+1};
52  end
53  if ischar(varargin{l}) && strcmpi(varargin{l},'precision')
54    tnargin = min(l,tnargin-1);
55    precision = varargin{l+1};
56  end
57  if ischar(varargin{l}) && strcmpi(varargin{l},'cell')
58    tnargin = min(l,tnargin-1);
59    arecell = varargin{l+1};
60  end
61  if ischar(varargin{l}) && strcmpi(varargin{l},'complex')
62    tnargin = min(l,tnargin-1);
63    arecomplex = varargin{l+1};
64  end
65end
66
67if strcmp(precision,'float128')
68  precision = 'float64';
69  system(['./convert -f ' inarg]);
70  fd = PetscOpenFile([inarg '_double']);
71end
72
73if arecell
74  narg = arecell;
75  rsult = cell(1);
76else
77  narg = nargout;
78end
79
80for l=1:narg
81  header = double(read(fd,1,indices));
82  if isempty(header)
83    if arecell
84      varargout(1) = {result};
85      return
86    else
87      disp('File/Socket does not have that many items')
88    end
89    return
90  end
91  if header == 1211216 % PETSc Mat Object
92
93    header = double(read(fd,3,indices));
94    m      = header(1);
95    n      = header(2);
96    nz     = header(3);
97    if (nz == -1)
98      if arecomplex
99        s     = read(fd,2*m*n,precision);
100        iReal = 1:2:n*m*2-1;
101        iImag = iReal +1 ;
102        A     = complex(reshape(s(iReal),n,m)',reshape(s(iImag),n,m)') ;
103      else
104        s   = read(fd,m*n,precision);
105        A   = reshape(s,n,m)';
106      end
107    else
108      nnz = double(read(fd,m,indices));  %nonzeros per row
109      sum_nz = sum(nnz);
110      if(sum_nz ~=nz)
111        str = sprintf('No-Nonzeros sum-rowlengths do not match %d %d',nz,sum_nz);
112        error(str);
113      end
114      j   = double(read(fd,nz,indices)) + 1;
115      if arecomplex
116        s   = read(fd,2*nz,precision);
117      else
118        s   = read(fd,nz,precision);
119      end
120      i   = ones(nz,1);
121      cnt = 1;
122      for k=1:m
123        next = cnt+nnz(k)-1;
124        i(cnt:next,1) = (double(k))*ones(nnz(k),1);
125        cnt = next+1;
126      end
127      if arecomplex
128        A = sparse(i,j,complex(s(1:2:2*nz),s(2:2:2*nz)),m,n,nz);
129      else
130        A = sparse(i,j,s,m,n,nz);
131      end
132    end
133    if arecell
134      result{l} = A;
135    else
136      varargout(l) = {A};
137    end
138  elseif  header == 1211214 % PETSc Vec Object
139    m = double(read(fd,1,indices));
140    if arecomplex
141      v = read(fd,2*m,precision);
142      v = complex(v(1:2:2*m),v(2:2:2*m));
143    else
144      v = read(fd,m,precision);
145    end
146    if arecell
147      result{l} = v;
148    else
149      varargout(l) = {v};
150    end
151
152  elseif  header == 1211213 % single real number
153    v = read(fd,1,precision);
154
155    if arecell
156      result{l} = v;
157    else
158      varargout(l) = {v};
159    end
160
161  elseif  header == 1211218 % PETSc IS Object
162    m = double(read(fd,1,indices));
163    v = read(fd,m,'int') + 1; % Indexing in MATLAB starts at 1, 0 in PETSc
164    if arecell
165      result{l} = v;
166    else
167      varargout(l) = {v};
168    end
169
170  elseif header == 1211219 % PETSc Bag Object
171    b = PetscBagRead(fd);
172    if arecell
173      result{l} = b;
174    else
175      varargout(l) = {b};
176    end
177
178  elseif header == 1211221 % PETSc DM Object
179    m  = double(read(fd,7,indices));
180    me = double(read(fd,5,indices));
181    b = [' dm ' int2str(m(3)) ' by ' int2str(m(4)) ' by ' int2str(m(5))];
182    if arecell
183      result{l} = b;
184    else
185      varargout(l) = {b};
186    end
187
188  else
189    disp(['Found unrecognized header ' int2str(header) ' in file. If your file contains complex numbers'])
190    disp(' then call PetscBinaryRead() with "complex",true as two additional arguments')
191    return
192  end
193
194end
195
196if arecell
197  varargout(1) = {result};
198end
199
200% close the reader if we opened it
201
202if nargin > 0
203  if (ischar(inarg) || isinteger(inarg)) close(fd); end;
204end
205