1: function [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: %
30: if nargin == 0
31: fd = PetscOpenSocket();
32: else if ischar(inarg) 33: fd = PetscOpenFile(inarg);
34: else if isnumeric(inarg) 35: fd = PetscOpenSocket(inarg);
36: else % assume it is a PetscOpenFile or PetscOpenSocket object and handles read()
37: fd = inarg;
38: end
39: end
40: end
42: indices = 'int32';
43: precision = 'float64';
44: arecell = 0;
45: arecomplex = false;
47: tnargin = nargin;
48: for 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
65: end
67: if strcmp(precision,'float128')
68: precision = 'float64';
69: system(['./convert -f ' inarg]);
70: fd = PetscOpenFile([inarg '_double']);
71: end
72: 73: if arecell
74: narg = arecell;
75: rsult = cell(1);
76: else
77: narg = nargout;
78: end
80: for 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
152: elseif header == 1211213 % single real number
153: v = read(fd,1,precision);
155: if arecell
156: result{l} = v;
157: else
158: varargout(l) = {v};
159: end
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
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
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
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
194: end
196: if arecell
197: varargout(1) = {result};
198: end
200: % close the reader if we opened it
202: if nargin > 0
203: if (ischar(inarg) || isinteger(inarg)) close(fd); end;
204: end