xref: /petsc/share/petsc/matlab/PetscBinaryWrite.m (revision 89a64d68d2f8dc8c4ae19e87b21a1e0b1231e51b) 
1function PetscBinaryWrite(inarg,varargin)
2%
3%  Writes in PETSc binary file sparse matrices and vectors.
4%  If the array is multidimensional and dense it is saved
5%  as a one dimensional PETSc Vec. If you want to save the multidimensional
6%  array as a matrix that MatLoad() will read you must first convert it to
7%  a sparse matrix: for example PetscBinaryWrite('myfile',sparse(A));
8%
9%
10%   PetscBinaryWrite(inarg,args to write,['indices','int32' or 'int64'],['precision','float64' or 'float32'],['complex',true,false])
11%   inarg may be:
12%      filename
13%      socket number (0 for PETSc default)
14%      the object returned from PetscOpenSocket or PetscOpenFile
15%
16if ischar(inarg)
17  fd = PetscOpenFile(inarg,'w');
18else if isnumeric(inarg)
19  if inarg == 0
20    fd = PetscOpenSocket;
21  else
22    fd = PetscOpenSocket(inarg);
23  end
24else
25  fd = inarg;
26end
27end
28
29indices = 'int32';
30precision = 'float64';
31ispetsccomplex = false;
32tnargin = nargin;
33for l=1:nargin-2
34  if ischar(varargin{l}) && strcmpi(varargin{l},'indices')
35    tnargin = min(l,tnargin-1);
36    indices = varargin{l+1};
37  end
38  if ischar(varargin{l}) && strcmpi(varargin{l},'precision')
39    tnargin = min(l,tnargin-1);
40    precision = varargin{l+1};
41  end
42  if ischar(varargin{l}) && strcmpi(varargin{l},'complex')
43    tnargin = min(l,tnargin-1);
44    ispetsccomplex = varargin{l+1};
45  end
46end
47
48for l=1:nargin-1
49  A = varargin{l};
50  if issparse(A) || min(size(A)) > 1
51    % save sparse matrix in special MATLAB format
52    if ~issparse(A)
53        A = sparse(A);
54    end
55    [m,n] = size(A);
56
57    if min(size(A)) == 1     %a one-rank matrix will be compressed to a
58                             %scalar instead of a vectory by sum
59      n_nz = full(A' ~= 0);
60    else
61      n_nz = full(sum(A' ~= 0));
62    end
63    nz   = sum(n_nz);
64    write(fd,[1211216,m,n,nz],indices);
65
66    write(fd,n_nz,indices);   %nonzeros per row
67    [i,j,s] = find(A');
68    write(fd,i-1,indices);
69    if ~isreal(s) || ispetsccomplex
70      s = conj(s);
71      ll = length(s);
72      sr = real(s);
73      si = imag(s);
74      s(1:2:2*ll) = sr;
75      s(2:2:2*ll) = si;
76    end
77    write(fd,s,precision);
78  else
79    [m,n] = size(A);
80    if isinteger(A)
81      classid = 1211218;
82    else
83      classid = 1211214;
84    end
85    write(fd,[classid,m*n],indices);
86    if ~isinteger(A) && (~isreal(A) || ispetsccomplex)
87      ll = length(A);
88      sr = real(A);
89      si = imag(A);
90      A(1:2:2*ll) = sr;
91      A(2:2:2*ll) = si;
92    end
93    if isinteger(A)
94       write(fd,A-1,indices); % indices starts at 1
95    else
96       write(fd,A,precision);
97    end
98  end
99end
100if ischar(inarg) || isinteger(inarg)
101    close(fd)
102end
103