GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
ls-mat4.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1996-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25 
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29 
30 #include <iomanip>
31 #include <istream>
32 #include <ostream>
33 #include <string>
34 
35 #include "byte-swap.h"
36 #include "dMatrix.h"
37 #include "dSparse.h"
38 #include "data-conv.h"
39 #include "file-ops.h"
40 #include "glob-match.h"
41 #include "lo-mappers.h"
42 #include "mach-info.h"
43 #include "oct-env.h"
44 #include "oct-locbuf.h"
45 #include "oct-time.h"
46 #include "quit.h"
47 
48 #include "ls-mat4.h"
49 #include "Cell.h"
50 #include "defun.h"
51 #include "error.h"
52 #include "errwarn.h"
53 #include "load-save.h"
54 #include "oct-map.h"
55 #include "ov-cell.h"
56 #include "ovl.h"
57 #include "pager.h"
58 #include "sysdep.h"
59 #include "utils.h"
60 #include "variables.h"
61 #include "version.h"
62 
63 
64 // Read LEN elements of data from IS in the format specified by
65 // PRECISION, placing the result in DATA. If SWAP is TRUE, swap
66 // the bytes of each element before copying to DATA. FLT_FMT
67 // specifies the format of the data if we are reading floating point
68 // numbers.
69 
70 static void
71 read_mat_binary_data (std::istream& is, double *data, int precision,
72  int len, bool swap,
74 {
75  switch (precision)
76  {
77  case 0:
78  read_doubles (is, data, LS_DOUBLE, len, swap, flt_fmt);
79  break;
80 
81  case 1:
82  read_doubles (is, data, LS_FLOAT, len, swap, flt_fmt);
83  break;
84 
85  case 2:
86  read_doubles (is, data, LS_INT, len, swap, flt_fmt);
87  break;
88 
89  case 3:
90  read_doubles (is, data, LS_SHORT, len, swap, flt_fmt);
91  break;
92 
93  case 4:
94  read_doubles (is, data, LS_U_SHORT, len, swap, flt_fmt);
95  break;
96 
97  case 5:
98  read_doubles (is, data, LS_U_CHAR, len, swap, flt_fmt);
99  break;
100 
101  default:
102  break;
103  }
104 }
105 
106 int
107 read_mat_file_header (std::istream& is, bool& swap, int32_t& mopt,
108  int32_t& nr, int32_t& nc,
109  int32_t& imag, int32_t& len,
110  int quiet)
111 {
112  swap = false;
113 
114  // We expect to fail here, at the beginning of a record, so not
115  // being able to read another mopt value should not result in an
116  // error.
117 
118  is.read (reinterpret_cast<char *> (&mopt), 4);
119  if (! is)
120  return 1;
121 
122  if (! is.read (reinterpret_cast<char *> (&nr), 4))
123  return -1;
124 
125  if (! is.read (reinterpret_cast<char *> (&nc), 4))
126  return -1;
127 
128  if (! is.read (reinterpret_cast<char *> (&imag), 4))
129  return -1;
130 
131  if (! is.read (reinterpret_cast<char *> (&len), 4))
132  return -1;
133 
134 // If mopt is nonzero and the byte order is swapped, mopt will be
135 // bigger than we expect, so we swap bytes.
136 //
137 // If mopt is zero, it means the file was written on a little endian machine,
138 // and we only need to swap if we are running on a big endian machine.
139 //
140 // Gag me.
141 
142  if (octave::mach_info::words_big_endian () && mopt == 0)
143  swap = true;
144 
145  // mopt is signed, therefore byte swap may result in negative value.
146 
147  if (mopt > 9999 || mopt < 0)
148  swap = true;
149 
150  if (swap)
151  {
152  swap_bytes<4> (&mopt);
153  swap_bytes<4> (&nr);
154  swap_bytes<4> (&nc);
155  swap_bytes<4> (&imag);
156  swap_bytes<4> (&len);
157  }
158 
159  if (mopt > 9999 || mopt < 0 || imag > 1 || imag < 0)
160  {
161  if (! quiet)
162  error ("load: can't read binary file");
163 
164  return -1;
165  }
166 
167  return 0;
168 }
169 
170 // We don't just use a cast here, because we need to be able to detect
171 // possible errors.
172 
175 {
177 
178  switch (mach)
179  {
180  case 0:
182  break;
183 
184  case 1:
186  break;
187 
188  case 2:
189  case 3:
190  case 4:
191  default:
193  break;
194  }
195 
196  return flt_fmt;
197 }
198 
199 int
201 {
202  int retval = -1;
203 
204  switch (flt_fmt)
205  {
207  retval = 0;
208  break;
209 
211  retval = 1;
212  break;
213 
214  default:
215  break;
216  }
217 
218  return retval;
219 }
220 
221 // Extract one value (scalar, matrix, string, etc.) from stream IS and
222 // place it in TC, returning the name of the variable.
223 //
224 // The data is expected to be in Matlab version 4 .mat format, though
225 // not all the features of that format are supported.
226 //
227 // FILENAME is used for error messages.
228 //
229 // This format provides no way to tag the data as global.
230 
231 std::string
232 read_mat_binary_data (std::istream& is, const std::string& filename,
233  octave_value& tc)
234 {
235  std::string retval;
236 
237  bool swap = false;
238  int32_t mopt, nr, nc, imag, len;
239 
240  int err = read_mat_file_header (is, swap, mopt, nr, nc, imag, len);
241  if (err)
242  {
243  if (err < 0)
244  error ("load: trouble reading binary file '%s'", filename.c_str ());
245 
246  return retval;
247  }
248 
249  int type = 0;
250  int prec = 0;
251  int order = 0;
252  int mach = 0;
253 
254  type = mopt % 10; // Full, sparse, etc.
255  mopt /= 10; // Eliminate first digit.
256  prec = mopt % 10; // double, float, int, etc.
257  mopt /= 10; // Eliminate second digit.
258  order = mopt % 10; // Row or column major ordering.
259  mopt /= 10; // Eliminate third digit.
260  mach = mopt % 10; // IEEE, VAX, etc.
261 
263  flt_fmt = mopt_digit_to_float_format (mach);
264 
265  if (flt_fmt == octave::mach_info::flt_fmt_unknown)
266  error ("load: unrecognized binary format!");
267 
268  if (imag && type == 1)
269  error ("load: encountered complex matrix with string flag set!");
270 
271  int dlen = 0;
272 
273  // LEN includes the terminating character, and the file is also
274  // supposed to include it, but apparently not all files do. Either
275  // way, I think this should work.
276 
277  {
278  OCTAVE_LOCAL_BUFFER (char, name, len+1);
279  name[len] = '\0';
280  if (! is.read (name, len))
281  error ("load: trouble reading binary file '%s'", filename.c_str ());
282  retval = name;
283 
284  dlen = nr * nc;
285  if (dlen < 0)
286  error ("load: trouble reading binary file '%s'", filename.c_str ());
287 
288  if (order)
289  {
290  octave_idx_type tmp = nr;
291  nr = nc;
292  nc = tmp;
293  }
294 
295  if (type == 2)
296  {
297  if (nc == 4)
298  {
299  octave_idx_type nr_new, nc_new;
300  Array<Complex> data (dim_vector (1, nr - 1));
301  Array<octave_idx_type> c (dim_vector (1, nr - 1));
302  Array<octave_idx_type> r (dim_vector (1, nr - 1));
303  OCTAVE_LOCAL_BUFFER (double, dtmp, nr);
304  OCTAVE_LOCAL_BUFFER (double, ctmp, nr);
305 
306  read_mat_binary_data (is, dtmp, prec, nr, swap, flt_fmt);
307  for (octave_idx_type i = 0; i < nr - 1; i++)
308  r.xelem (i) = dtmp[i] - 1;
309  nr_new = dtmp[nr - 1];
310  read_mat_binary_data (is, dtmp, prec, nr, swap, flt_fmt);
311  for (octave_idx_type i = 0; i < nr - 1; i++)
312  c.xelem (i) = dtmp[i] - 1;
313  nc_new = dtmp[nr - 1];
314  read_mat_binary_data (is, dtmp, prec, nr - 1, swap, flt_fmt);
315  read_mat_binary_data (is, ctmp, prec, 1, swap, flt_fmt);
316  read_mat_binary_data (is, ctmp, prec, nr - 1, swap, flt_fmt);
317 
318  for (octave_idx_type i = 0; i < nr - 1; i++)
319  data.xelem (i) = Complex (dtmp[i], ctmp[i]);
320  read_mat_binary_data (is, ctmp, prec, 1, swap, flt_fmt);
321 
322  SparseComplexMatrix smc = SparseComplexMatrix (data, r, c,
323  nr_new, nc_new);
324 
325  tc = (order ? smc.transpose () : smc);
326  }
327  else
328  {
329  octave_idx_type nr_new, nc_new;
330  Array<double> data (dim_vector (1, nr - 1));
331  Array<octave_idx_type> c (dim_vector (1, nr - 1));
332  Array<octave_idx_type> r (dim_vector (1, nr - 1));
333  OCTAVE_LOCAL_BUFFER (double, dtmp, nr);
334 
335  read_mat_binary_data (is, dtmp, prec, nr, swap, flt_fmt);
336  for (octave_idx_type i = 0; i < nr - 1; i++)
337  r.xelem (i) = dtmp[i] - 1;
338  nr_new = dtmp[nr - 1];
339  read_mat_binary_data (is, dtmp, prec, nr, swap, flt_fmt);
340  for (octave_idx_type i = 0; i < nr - 1; i++)
341  c.xelem (i) = dtmp[i] - 1;
342  nc_new = dtmp[nr - 1];
343  read_mat_binary_data (is, data.fortran_vec (), prec, nr - 1,
344  swap, flt_fmt);
345  read_mat_binary_data (is, dtmp, prec, 1, swap, flt_fmt);
346 
347  SparseMatrix sm = SparseMatrix (data, r, c, nr_new, nc_new);
348 
349  tc = (order ? sm.transpose () : sm);
350  }
351  }
352  else
353  {
354  Matrix re (nr, nc);
355 
356  read_mat_binary_data (is, re.fortran_vec (), prec, dlen, swap, flt_fmt);
357 
358  if (! is)
359  error ("load: reading matrix data for '%s'", name);
360 
361  if (imag)
362  {
363  Matrix im (nr, nc);
364 
365  read_mat_binary_data (is, im.fortran_vec (), prec, dlen, swap,
366  flt_fmt);
367 
368  if (! is)
369  error ("load: reading imaginary matrix data for '%s'", name);
370 
371  ComplexMatrix ctmp (nr, nc);
372 
373  for (octave_idx_type j = 0; j < nc; j++)
374  for (octave_idx_type i = 0; i < nr; i++)
375  ctmp (i,j) = Complex (re(i,j), im(i,j));
376 
377  tc = (order ? ctmp.transpose () : ctmp);
378  }
379  else
380  tc = (order ? re.transpose () : re);
381 
382  if (type == 1)
383  tc = tc.convert_to_str (false, true, '\'');
384  }
385 
386  return retval;
387  }
388 }
389 
390 // Save the data from TC along with the corresponding NAME on stream OS
391 // in the MatLab version 4 binary format.
392 
393 bool
394 save_mat_binary_data (std::ostream& os, const octave_value& tc,
395  const std::string& name)
396 {
397  int32_t mopt = 0;
398 
399  mopt += tc.issparse () ? 2 : tc.is_string () ? 1 : 0;
400 
403 
404  mopt += 1000 * float_format_to_mopt_digit (flt_fmt);
405 
406  os.write (reinterpret_cast<char *> (&mopt), 4);
407 
409  int32_t nr = tc.rows ();
410 
411  int32_t nc = tc.columns ();
412 
413  if (tc.issparse ())
414  {
415  len = tc.nnz ();
416  uint32_t nnz = len + 1;
417  os.write (reinterpret_cast<char *> (&nnz), 4);
418 
419  uint32_t iscmplx = (tc.iscomplex () ? 4 : 3);
420  os.write (reinterpret_cast<char *> (&iscmplx), 4);
421 
422  uint32_t tmp = 0;
423  os.write (reinterpret_cast<char *> (&tmp), 4);
424  }
425  else
426  {
427  os.write (reinterpret_cast<char *> (&nr), 4);
428  os.write (reinterpret_cast<char *> (&nc), 4);
429 
430  int32_t imag = (tc.iscomplex () ? 1 : 0);
431  os.write (reinterpret_cast<char *> (&imag), 4);
432 
433  len = nr * nc;
434  }
435 
436  // LEN includes the terminating character, and the file is also
437  // supposed to include it.
438 
439  int32_t name_len = name.length () + 1;
440 
441  os.write (reinterpret_cast<char *> (&name_len), 4);
442  os << name << '\0';
443 
444  if (tc.is_string ())
445  {
446  charMatrix chm = tc.char_matrix_value ();
447 
448  octave_idx_type nrow = chm.rows ();
449  octave_idx_type ncol = chm.cols ();
450 
451  OCTAVE_LOCAL_BUFFER (double, buf, ncol*nrow);
452 
453  for (octave_idx_type i = 0; i < nrow; i++)
454  {
455  std::string tstr = chm.row_as_string (i);
456  const char *s = tstr.data ();
457 
458  for (octave_idx_type j = 0; j < ncol; j++)
459  buf[j*nrow+i] = static_cast<double> (*s++ & 0x00FF);
460  }
461  std::streamsize n_bytes = static_cast<std::streamsize> (nrow) *
462  static_cast<std::streamsize> (ncol) *
463  sizeof (double);
464  os.write (reinterpret_cast<char *> (buf), n_bytes);
465  }
466  else if (tc.is_range ())
467  {
468  Range r = tc.range_value ();
469  double base = r.base ();
470  double inc = r.inc ();
471  octave_idx_type nel = r.numel ();
472  for (octave_idx_type i = 0; i < nel; i++)
473  {
474  double x = base + i * inc;
475  os.write (reinterpret_cast<char *> (&x), 8);
476  }
477  }
478  else if (tc.is_real_scalar ())
479  {
480  double tmp = tc.double_value ();
481  os.write (reinterpret_cast<char *> (&tmp), 8);
482  }
483  else if (tc.issparse ())
484  {
485  double ds;
486  OCTAVE_LOCAL_BUFFER (double, dtmp, len);
487  if (tc.is_complex_matrix ())
488  {
490 
491  for (octave_idx_type i = 0; i < len; i++)
492  dtmp[i] = m.ridx (i) + 1;
493  std::streamsize n_bytes = 8 * static_cast<std::streamsize> (len);
494  os.write (reinterpret_cast<const char *> (dtmp), n_bytes);
495  ds = nr;
496  os.write (reinterpret_cast<const char *> (&ds), 8);
497 
498  octave_idx_type ii = 0;
499  for (octave_idx_type j = 0; j < nc; j++)
500  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
501  dtmp[ii++] = j + 1;
502  os.write (reinterpret_cast<const char *> (dtmp), n_bytes);
503  ds = nc;
504  os.write (reinterpret_cast<const char *> (&ds), 8);
505 
506  for (octave_idx_type i = 0; i < len; i++)
507  dtmp[i] = std::real (m.data (i));
508  os.write (reinterpret_cast<const char *> (dtmp), n_bytes);
509  ds = 0.;
510  os.write (reinterpret_cast<const char *> (&ds), 8);
511 
512  for (octave_idx_type i = 0; i < len; i++)
513  dtmp[i] = std::imag (m.data (i));
514  os.write (reinterpret_cast<const char *> (dtmp), n_bytes);
515  os.write (reinterpret_cast<const char *> (&ds), 8);
516  }
517  else
518  {
520 
521  for (octave_idx_type i = 0; i < len; i++)
522  dtmp[i] = m.ridx (i) + 1;
523  std::streamsize n_bytes = 8 * static_cast<std::streamsize> (len);
524  os.write (reinterpret_cast<const char *> (dtmp), n_bytes);
525  ds = nr;
526  os.write (reinterpret_cast<const char *> (&ds), 8);
527 
528  octave_idx_type ii = 0;
529  for (octave_idx_type j = 0; j < nc; j++)
530  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
531  dtmp[ii++] = j + 1;
532  os.write (reinterpret_cast<const char *> (dtmp), n_bytes);
533  ds = nc;
534  os.write (reinterpret_cast<const char *> (&ds), 8);
535 
536  os.write (reinterpret_cast<const char *> (m.data ()), n_bytes);
537  ds = 0.;
538  os.write (reinterpret_cast<const char *> (&ds), 8);
539  }
540  }
541  else if (tc.is_real_matrix ())
542  {
543  Matrix m = tc.matrix_value ();
544  std::streamsize n_bytes = 8 * static_cast<std::streamsize> (len);
545  os.write (reinterpret_cast<const char *> (m.data ()), n_bytes);
546  }
547  else if (tc.is_complex_scalar ())
548  {
549  Complex tmp = tc.complex_value ();
550  os.write (reinterpret_cast<char *> (&tmp), 16);
551  }
552  else if (tc.is_complex_matrix ())
553  {
554  ComplexMatrix m_cmplx = tc.complex_matrix_value ();
555  Matrix m = ::real (m_cmplx);
556  std::streamsize n_bytes = 8 * static_cast<std::streamsize> (len);
557  os.write (reinterpret_cast<const char *> (m.data ()), n_bytes);
558  m = ::imag (m_cmplx);
559  os.write (reinterpret_cast<const char *> (m.data ()), n_bytes);
560  }
561  else
562  // FIXME: Should this just error out rather than warn?
563  warn_wrong_type_arg ("save", tc);
564 
565  return ! os.fail ();
566 }
void swap_bytes< 4 >(void *ptr)
Definition: byte-swap.h:63
N Dimensional Array with copy-on-write semantics.
Definition: Array.h:128
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:469
octave_idx_type cols(void) const
Definition: Array.h:423
octave_idx_type rows(void) const
Definition: Array.h:415
const T * fortran_vec(void) const
Size of the specified dimension.
Definition: Array.h:583
ComplexMatrix transpose(void) const
Definition: CMatrix.h:171
Definition: dMatrix.h:42
Matrix transpose(void) const
Definition: dMatrix.h:135
Definition: Range.h:40
SparseComplexMatrix transpose(void) const
Definition: CSparse.h:140
SparseMatrix transpose(void) const
Definition: dSparse.h:128
std::string row_as_string(octave_idx_type, bool strip_ws=false) const
Definition: chMatrix.cc:82
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:95
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:853
bool is_complex_matrix(void) const
Definition: ov.h:575
bool issparse(void) const
Definition: ov.h:706
charMatrix char_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:847
octave_idx_type nnz(void) const
Definition: ov.h:524
ComplexMatrix complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:824
Complex complex_value(bool frc_str_conv=false) const
Definition: ov.h:818
octave_idx_type rows(void) const
Definition: ov.h:504
Range range_value(void) const
Definition: ov.h:938
bool is_string(void) const
Definition: ov.h:593
octave_idx_type columns(void) const
Definition: ov.h:506
bool is_range(void) const
Definition: ov.h:602
bool is_real_matrix(void) const
Definition: ov.h:569
octave_value convert_to_str(bool pad=false, bool force=false, char type='\'') const
Definition: ov.h:1216
bool is_real_scalar(void) const
Definition: ov.h:566
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:806
bool iscomplex(void) const
Definition: ov.h:694
double double_value(bool frc_str_conv=false) const
Definition: ov.h:794
bool is_complex_scalar(void) const
Definition: ov.h:572
SparseComplexMatrix sparse_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:857
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:137
ColumnVector imag(const ComplexColumnVector &a)
Definition: dColVector.cc:143
void read_doubles(std::istream &is, double *data, save_type type, octave_idx_type len, bool swap, octave::mach_info::float_format fmt)
Definition: data-conv.cc:776
@ LS_U_CHAR
Definition: data-conv.h:88
@ LS_DOUBLE
Definition: data-conv.h:95
@ LS_U_SHORT
Definition: data-conv.h:89
@ LS_FLOAT
Definition: data-conv.h:94
@ LS_SHORT
Definition: data-conv.h:92
@ LS_INT
Definition: data-conv.h:93
void error(const char *fmt,...)
Definition: error.cc:968
void warn_wrong_type_arg(const char *name, const octave_value &tc)
Definition: errwarn.cc:372
QString name
F77_RET_T const F77_DBLE * x
octave::mach_info::float_format mopt_digit_to_float_format(int mach)
Definition: ls-mat4.cc:174
bool save_mat_binary_data(std::ostream &os, const octave_value &tc, const std::string &name)
Definition: ls-mat4.cc:394
int read_mat_file_header(std::istream &is, bool &swap, int32_t &mopt, int32_t &nr, int32_t &nc, int32_t &imag, int32_t &len, int quiet)
Definition: ls-mat4.cc:107
static void read_mat_binary_data(std::istream &is, double *data, int precision, int len, bool swap, octave::mach_info::float_format flt_fmt)
Definition: ls-mat4.cc:71
int float_format_to_mopt_digit(octave::mach_info::float_format flt_fmt)
Definition: ls-mat4.cc:200
T octave_idx_type m
Definition: mx-inlines.cc:773
T * r
Definition: mx-inlines.cc:773
float_format native_float_format(void)
Definition: mach-info.cc:65
@ flt_fmt_ieee_little_endian
Definition: mach-info.h:43
bool words_big_endian(void)
Definition: mach-info.cc:72
std::complex< double > Complex
Definition: oct-cmplx.h:33
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:44
octave_value::octave_value(const Array< char > &chm, char type) return retval
Definition: ov.cc:811
F77_RET_T len
Definition: xerbla.cc:61