GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ov-cx-sparse.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2004-2013 David Bateman
4 Copyright (C) 1998-2004 Andy Adler
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include <iostream>
29 #include <limits>
30 #include <vector>
31 
32 #include "lo-specfun.h"
33 #include "lo-mappers.h"
34 #include "oct-locbuf.h"
35 
36 #include "mxarray.h"
37 #include "ov-base.h"
38 #include "ov-scalar.h"
39 #include "ov-complex.h"
40 #include "gripes.h"
41 
42 #include "ov-re-sparse.h"
43 #include "ov-cx-sparse.h"
44 
45 #include "ov-base-sparse.h"
46 #include "ov-base-sparse.cc"
47 
48 #include "ov-bool-sparse.h"
49 
51 
53 
55  "sparse complex matrix", "double");
56 
59 {
60  octave_base_value *retval = 0;
61 
63  {
64  int nr = matrix.rows ();
65  int nc = matrix.cols ();
66 
67  // Don't use numel, since it can overflow for very large matrices
68  // Note that for the tests on matrix size, they become approximative
69  // since they involves a cast to double to avoid issues of overflow
70  if (matrix.rows () == 1 && matrix.cols () == 1)
71  {
72  // Const copy of the matrix, so the right version of () operator used
73  const SparseComplexMatrix tmp (matrix);
74 
75  Complex c = tmp (0, 0);
76 
77  if (std::imag (c) == 0.0)
78  retval = new octave_scalar (std::real (c));
79  else
80  retval = new octave_complex (c);
81  }
82  else if (nr == 0 || nc == 0)
83  retval = new octave_matrix (Matrix (nr, nc));
84  else if (matrix.all_elements_are_real ())
85  if (matrix.cols () > 0 && matrix.rows () > 0
86  && (double (matrix.byte_size ()) > double (matrix.rows ())
87  * double (matrix.cols ()) * sizeof (double)))
88  retval = new octave_matrix (::real (matrix.matrix_value ()));
89  else
90  retval = new octave_sparse_matrix (::real (matrix));
91  else if (matrix.cols () > 0 && matrix.rows () > 0
92  && (double (matrix.byte_size ()) > double (matrix.rows ())
93  * double (matrix.cols ()) * sizeof (Complex)))
94  retval = new octave_complex_matrix (matrix.matrix_value ());
95  }
96  else
97  {
99  retval = new octave_sparse_matrix (::real (matrix));
100  }
101 
102  return retval;
103 }
104 
105 double
106 octave_sparse_complex_matrix::double_value (bool force_conversion) const
107 {
108  double retval = lo_ieee_nan_value ();
109 
110  if (! force_conversion)
111  gripe_implicit_conversion ("Octave:imag-to-real",
112  "complex sparse matrix", "real scalar");
113 
114  // FIXME: maybe this should be a function, valid_as_scalar()
115  if (numel () > 0)
116  {
117  if (numel () > 1)
118  gripe_implicit_conversion ("Octave:array-to-scalar",
119  "complex sparse matrix", "real scalar");
120 
121  retval = std::real (matrix (0, 0));
122  }
123  else
124  gripe_invalid_conversion ("complex sparse matrix", "real scalar");
125 
126  return retval;
127 }
128 
129 Matrix
130 octave_sparse_complex_matrix::matrix_value (bool force_conversion) const
131 {
132  Matrix retval;
133 
134  if (! force_conversion)
135  gripe_implicit_conversion ("Octave:imag-to-real",
136  "complex sparse matrix", "real matrix");
137 
138  retval = ::real (matrix.matrix_value ());
139 
140  return retval;
141 }
142 
143 Complex
145 {
146  double tmp = lo_ieee_nan_value ();
147 
148  Complex retval (tmp, tmp);
149 
150  // FIXME: maybe this should be a function, valid_as_scalar()
151  if (numel () > 0)
152  {
153  if (numel () > 1)
154  gripe_implicit_conversion ("Octave:array-to-scalar",
155  "complex sparse matrix", "real scalar");
156 
157  retval = matrix (0, 0);
158  }
159  else
160  gripe_invalid_conversion ("complex sparse matrix", "real scalar");
161 
162  return retval;
163 }
164 
167 {
168  return matrix.matrix_value ();
169 }
170 
173 {
174  return ComplexNDArray (matrix.matrix_value ());
175 }
176 
179 {
180  charNDArray retval;
181 
182  if (! frc_str_conv)
183  gripe_implicit_conversion ("Octave:num-to-str",
184  "sparse complex matrix", "string");
185  else
186  {
187  retval = charNDArray (dims (), 0);
188  octave_idx_type nc = matrix.cols ();
189  octave_idx_type nr = matrix.rows ();
190 
191  for (octave_idx_type j = 0; j < nc; j++)
192  for (octave_idx_type i = matrix.cidx (j); i < matrix.cidx (j+1); i++)
193  retval(matrix.ridx (i) + nr * j) =
194  static_cast<char>(std::real (matrix.data (i)));
195  }
196 
197  return retval;
198 }
199 
202 {
203  SparseMatrix retval;
204 
205  if (! force_conversion)
206  gripe_implicit_conversion ("Octave:imag-to-real",
207  "complex sparse matrix",
208  "real sparse matrix");
209 
210  retval = ::real (matrix);
211 
212  return retval;
213 }
214 
217 {
218  if (matrix.any_element_is_nan ())
220  else if (warn && (! matrix.all_elements_are_real ()
221  || real (matrix).any_element_not_one_or_zero ()))
223 
224  return mx_el_ne (matrix, Complex (0.0));
225 }
226 
227 bool
229  bool&save_as_floats)
230 {
231  dim_vector d = this->dims ();
232  if (d.length () < 1)
233  return false;
234 
235  // Ensure that additional memory is deallocated
237 
238  int nr = d(0);
239  int nc = d(1);
240  int nz = nnz ();
241 
242  int32_t itmp;
243  // Use negative value for ndims to be consistent with other formats
244  itmp= -2;
245  os.write (reinterpret_cast<char *> (&itmp), 4);
246 
247  itmp= nr;
248  os.write (reinterpret_cast<char *> (&itmp), 4);
249 
250  itmp= nc;
251  os.write (reinterpret_cast<char *> (&itmp), 4);
252 
253  itmp= nz;
254  os.write (reinterpret_cast<char *> (&itmp), 4);
255 
256  save_type st = LS_DOUBLE;
257  if (save_as_floats)
258  {
260  {
261  warning ("save: some values too large to save as floats --");
262  warning ("save: saving as doubles instead");
263  }
264  else
265  st = LS_FLOAT;
266  }
267  else if (matrix.nnz () > 8192) // FIXME: make this configurable.
268  {
269  double max_val, min_val;
270  if (matrix.all_integers (max_val, min_val))
271  st = get_save_type (max_val, min_val);
272  }
273 
274  // add one to the printed indices to go from
275  // zero-based to one-based arrays
276  for (int i = 0; i < nc+1; i++)
277  {
278  octave_quit ();
279  itmp = matrix.cidx (i);
280  os.write (reinterpret_cast<char *> (&itmp), 4);
281  }
282 
283  for (int i = 0; i < nz; i++)
284  {
285  octave_quit ();
286  itmp = matrix.ridx (i);
287  os.write (reinterpret_cast<char *> (&itmp), 4);
288  }
289 
290  write_doubles (os, reinterpret_cast<const double *> (matrix.data ()), st,
291  2 * nz);
292 
293  return true;
294 }
295 
296 bool
297 octave_sparse_complex_matrix::load_binary (std::istream& is, bool swap,
299 {
300  int32_t nz, nc, nr, tmp;
301  char ctmp;
302 
303  if (! is.read (reinterpret_cast<char *> (&tmp), 4))
304  return false;
305 
306  if (swap)
307  swap_bytes<4> (&tmp);
308 
309  if (tmp != -2)
310  {
311  error ("load: only 2-D sparse matrices are supported");
312  return false;
313  }
314 
315  if (! is.read (reinterpret_cast<char *> (&nr), 4))
316  return false;
317  if (! is.read (reinterpret_cast<char *> (&nc), 4))
318  return false;
319  if (! is.read (reinterpret_cast<char *> (&nz), 4))
320  return false;
321 
322  if (swap)
323  {
324  swap_bytes<4> (&nr);
325  swap_bytes<4> (&nc);
326  swap_bytes<4> (&nz);
327  }
328 
329  SparseComplexMatrix m (static_cast<octave_idx_type> (nr),
330  static_cast<octave_idx_type> (nc),
331  static_cast<octave_idx_type> (nz));
332 
333  for (int i = 0; i < nc+1; i++)
334  {
335  octave_quit ();
336  if (! is.read (reinterpret_cast<char *> (&tmp), 4))
337  return false;
338  if (swap)
339  swap_bytes<4> (&tmp);
340  m.cidx (i) = tmp;
341  }
342 
343  for (int i = 0; i < nz; i++)
344  {
345  octave_quit ();
346  if (! is.read (reinterpret_cast<char *> (&tmp), 4))
347  return false;
348  if (swap)
349  swap_bytes<4> (&tmp);
350  m.ridx (i) = tmp;
351  }
352 
353  if (! is.read (reinterpret_cast<char *> (&ctmp), 1))
354  return false;
355 
356  read_doubles (is, reinterpret_cast<double *> (m.data ()),
357  static_cast<save_type> (ctmp), 2 * nz, swap, fmt);
358 
359  if (error_state || ! is)
360  return false;
361 
362  if (! m.indices_ok ())
363  return false;
364 
365  matrix = m;
366 
367  return true;
368 }
369 
370 #if defined (HAVE_HDF5)
371 
372 bool
373 octave_sparse_complex_matrix::save_hdf5 (hid_t loc_id, const char *name,
374  bool save_as_floats)
375 {
376  dim_vector dv = dims ();
377  int empty = save_hdf5_empty (loc_id, name, dv);
378  if (empty)
379  return (empty > 0);
380 
381  // Ensure that additional memory is deallocated
383 
384 #if HAVE_HDF5_18
385  hid_t group_hid = H5Gcreate (loc_id, name, H5P_DEFAULT, H5P_DEFAULT,
386  H5P_DEFAULT);
387 #else
388  hid_t group_hid = H5Gcreate (loc_id, name, 0);
389 #endif
390  if (group_hid < 0)
391  return false;
392 
393  hid_t space_hid = -1, data_hid = -1;
394  bool retval = true;
396  octave_idx_type tmp;
397  hsize_t hdims[2];
398 
399  space_hid = H5Screate_simple (0, hdims, 0);
400  if (space_hid < 0)
401  {
402  H5Gclose (group_hid);
403  return false;
404  }
405 
406 #if HAVE_HDF5_18
407  data_hid = H5Dcreate (group_hid, "nr", H5T_NATIVE_IDX, space_hid,
408  H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
409 #else
410  data_hid = H5Dcreate (group_hid, "nr", H5T_NATIVE_IDX, space_hid,
411  H5P_DEFAULT);
412 #endif
413  if (data_hid < 0)
414  {
415  H5Sclose (space_hid);
416  H5Gclose (group_hid);
417  return false;
418  }
419 
420  tmp = m.rows ();
421  retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL,
422  H5P_DEFAULT, &tmp) >= 0;
423  H5Dclose (data_hid);
424  if (!retval)
425  {
426  H5Sclose (space_hid);
427  H5Gclose (group_hid);
428  return false;
429  }
430 
431 #if HAVE_HDF5_18
432  data_hid = H5Dcreate (group_hid, "nc", H5T_NATIVE_IDX, space_hid,
433  H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
434 #else
435  data_hid = H5Dcreate (group_hid, "nc", H5T_NATIVE_IDX, space_hid,
436  H5P_DEFAULT);
437 #endif
438  if (data_hid < 0)
439  {
440  H5Sclose (space_hid);
441  H5Gclose (group_hid);
442  return false;
443  }
444 
445  tmp = m.cols ();
446  retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL,
447  H5P_DEFAULT, &tmp) >= 0;
448  H5Dclose (data_hid);
449  if (!retval)
450  {
451  H5Sclose (space_hid);
452  H5Gclose (group_hid);
453  return false;
454  }
455 
456 #if HAVE_HDF5_18
457  data_hid = H5Dcreate (group_hid, "nz", H5T_NATIVE_IDX, space_hid,
458  H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
459 #else
460  data_hid = H5Dcreate (group_hid, "nz", H5T_NATIVE_IDX, space_hid,
461  H5P_DEFAULT);
462 #endif
463  if (data_hid < 0)
464  {
465  H5Sclose (space_hid);
466  H5Gclose (group_hid);
467  return false;
468  }
469 
470  tmp = m.nnz ();
471  retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL,
472  H5P_DEFAULT, &tmp) >= 0;
473  H5Dclose (data_hid);
474  if (!retval)
475  {
476  H5Sclose (space_hid);
477  H5Gclose (group_hid);
478  return false;
479  }
480 
481  H5Sclose (space_hid);
482 
483  hdims[0] = m.cols () + 1;
484  hdims[1] = 1;
485 
486  space_hid = H5Screate_simple (2, hdims, 0);
487 
488  if (space_hid < 0)
489  {
490  H5Gclose (group_hid);
491  return false;
492  }
493 
494 #if HAVE_HDF5_18
495  data_hid = H5Dcreate (group_hid, "cidx", H5T_NATIVE_IDX, space_hid,
496  H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
497 #else
498  data_hid = H5Dcreate (group_hid, "cidx", H5T_NATIVE_IDX, space_hid,
499  H5P_DEFAULT);
500 #endif
501  if (data_hid < 0)
502  {
503  H5Sclose (space_hid);
504  H5Gclose (group_hid);
505  return false;
506  }
507 
508  octave_idx_type * itmp = m.xcidx ();
509  retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL,
510  H5P_DEFAULT, itmp) >= 0;
511  H5Dclose (data_hid);
512  if (!retval)
513  {
514  H5Sclose (space_hid);
515  H5Gclose (group_hid);
516  return false;
517  }
518 
519  H5Sclose (space_hid);
520 
521  hdims[0] = m.nnz ();
522  hdims[1] = 1;
523 
524  space_hid = H5Screate_simple (2, hdims, 0);
525 
526  if (space_hid < 0)
527  {
528  H5Gclose (group_hid);
529  return false;
530  }
531 
532 #if HAVE_HDF5_18
533  data_hid = H5Dcreate (group_hid, "ridx", H5T_NATIVE_IDX, space_hid,
534  H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
535 #else
536  data_hid = H5Dcreate (group_hid, "ridx", H5T_NATIVE_IDX, space_hid,
537  H5P_DEFAULT);
538 #endif
539  if (data_hid < 0)
540  {
541  H5Sclose (space_hid);
542  H5Gclose (group_hid);
543  return false;
544  }
545 
546  itmp = m.xridx ();
547  retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL, H5P_DEFAULT,
548  itmp) >= 0;
549  H5Dclose (data_hid);
550  if (!retval)
551  {
552  H5Sclose (space_hid);
553  H5Gclose (group_hid);
554  return false;
555  }
556 
557  hid_t save_type_hid = H5T_NATIVE_DOUBLE;
558 
559  if (save_as_floats)
560  {
561  if (m.too_large_for_float ())
562  {
563  warning ("save: some values too large to save as floats --");
564  warning ("save: saving as doubles instead");
565  }
566  else
567  save_type_hid = H5T_NATIVE_FLOAT;
568  }
569 #if HAVE_HDF5_INT2FLOAT_CONVERSIONS
570  // hdf5 currently doesn't support float/integer conversions
571  else
572  {
573  double max_val, min_val;
574 
575  if (m.all_integers (max_val, min_val))
576  save_type_hid
577  = save_type_to_hdf5 (get_save_type (max_val, min_val));
578  }
579 #endif /* HAVE_HDF5_INT2FLOAT_CONVERSIONS */
580 
581  hid_t type_hid = hdf5_make_complex_type (save_type_hid);
582  if (type_hid < 0)
583  {
584  H5Sclose (space_hid);
585  H5Gclose (group_hid);
586  return false;
587  }
588 #if HAVE_HDF5_18
589  data_hid = H5Dcreate (group_hid, "data", type_hid, space_hid,
590  H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
591 #else
592  data_hid = H5Dcreate (group_hid, "data", type_hid, space_hid, H5P_DEFAULT);
593 #endif
594  if (data_hid < 0)
595  {
596  H5Sclose (space_hid);
597  H5Tclose (type_hid);
598  H5Gclose (group_hid);
599  return false;
600  }
601 
602  hid_t complex_type_hid = hdf5_make_complex_type (H5T_NATIVE_DOUBLE);
603  retval = false;
604  if (complex_type_hid >= 0)
605  {
606  Complex * ctmp = m.xdata ();
607 
608  retval = H5Dwrite (data_hid, complex_type_hid, H5S_ALL, H5S_ALL,
609  H5P_DEFAULT, ctmp) >= 0;
610  }
611 
612  H5Dclose (data_hid);
613  H5Sclose (space_hid);
614  H5Tclose (type_hid);
615  H5Gclose (group_hid);
616 
617  return retval;
618 }
619 
620 bool
621 octave_sparse_complex_matrix::load_hdf5 (hid_t loc_id, const char *name)
622 {
623  octave_idx_type nr, nc, nz;
624  hid_t group_hid, data_hid, space_hid;
625  hsize_t rank;
626 
627  dim_vector dv;
628  int empty = load_hdf5_empty (loc_id, name, dv);
629  if (empty > 0)
630  matrix.resize (dv);
631  if (empty)
632  return (empty > 0);
633 
634 #if HAVE_HDF5_18
635  group_hid = H5Gopen (loc_id, name, H5P_DEFAULT);
636 #else
637  group_hid = H5Gopen (loc_id, name);
638 #endif
639  if (group_hid < 0 ) return false;
640 
641 #if HAVE_HDF5_18
642  data_hid = H5Dopen (group_hid, "nr", H5P_DEFAULT);
643 #else
644  data_hid = H5Dopen (group_hid, "nr");
645 #endif
646  space_hid = H5Dget_space (data_hid);
647  rank = H5Sget_simple_extent_ndims (space_hid);
648 
649  if (rank != 0)
650  {
651  H5Dclose (data_hid);
652  H5Gclose (group_hid);
653  return false;
654  }
655 
656  if (H5Dread (data_hid, H5T_NATIVE_IDX, H5S_ALL,
657  H5S_ALL, H5P_DEFAULT, &nr) < 0)
658  {
659  H5Dclose (data_hid);
660  H5Gclose (group_hid);
661  return false;
662  }
663 
664  H5Dclose (data_hid);
665 
666 #if HAVE_HDF5_18
667  data_hid = H5Dopen (group_hid, "nc", H5P_DEFAULT);
668 #else
669  data_hid = H5Dopen (group_hid, "nc");
670 #endif
671  space_hid = H5Dget_space (data_hid);
672  rank = H5Sget_simple_extent_ndims (space_hid);
673 
674  if (rank != 0)
675  {
676  H5Dclose (data_hid);
677  H5Gclose (group_hid);
678  return false;
679  }
680 
681  if (H5Dread (data_hid, H5T_NATIVE_IDX, H5S_ALL,
682  H5S_ALL, H5P_DEFAULT, &nc) < 0)
683  {
684  H5Dclose (data_hid);
685  H5Gclose (group_hid);
686  return false;
687  }
688 
689  H5Dclose (data_hid);
690 
691 #if HAVE_HDF5_18
692  data_hid = H5Dopen (group_hid, "nz", H5P_DEFAULT);
693 #else
694  data_hid = H5Dopen (group_hid, "nz");
695 #endif
696  space_hid = H5Dget_space (data_hid);
697  rank = H5Sget_simple_extent_ndims (space_hid);
698 
699  if (rank != 0)
700  {
701  H5Dclose (data_hid);
702  H5Gclose (group_hid);
703  return false;
704  }
705 
706  if (H5Dread (data_hid, H5T_NATIVE_IDX, H5S_ALL,
707  H5S_ALL, H5P_DEFAULT, &nz) < 0)
708  {
709  H5Dclose (data_hid);
710  H5Gclose (group_hid);
711  return false;
712  }
713 
714  H5Dclose (data_hid);
715 
716  SparseComplexMatrix m (static_cast<octave_idx_type> (nr),
717  static_cast<octave_idx_type> (nc),
718  static_cast<octave_idx_type> (nz));
719 
720 #if HAVE_HDF5_18
721  data_hid = H5Dopen (group_hid, "cidx", H5P_DEFAULT);
722 #else
723  data_hid = H5Dopen (group_hid, "cidx");
724 #endif
725  space_hid = H5Dget_space (data_hid);
726  rank = H5Sget_simple_extent_ndims (space_hid);
727 
728  if (rank != 2)
729  {
730  H5Sclose (space_hid);
731  H5Dclose (data_hid);
732  H5Gclose (group_hid);
733  return false;
734  }
735 
736  OCTAVE_LOCAL_BUFFER (hsize_t, hdims, rank);
737  OCTAVE_LOCAL_BUFFER (hsize_t, maxdims, rank);
738 
739  H5Sget_simple_extent_dims (space_hid, hdims, maxdims);
740 
741  if (static_cast<int> (hdims[0]) != nc + 1
742  || static_cast<int> (hdims[1]) != 1)
743  {
744  H5Sclose (space_hid);
745  H5Dclose (data_hid);
746  H5Gclose (group_hid);
747  return false;
748  }
749 
750  octave_idx_type *itmp = m.xcidx ();
751  if (H5Dread (data_hid, H5T_NATIVE_IDX, H5S_ALL,
752  H5S_ALL, H5P_DEFAULT, itmp) < 0)
753  {
754  H5Sclose (space_hid);
755  H5Dclose (data_hid);
756  H5Gclose (group_hid);
757  return false;
758  }
759 
760  H5Sclose (space_hid);
761  H5Dclose (data_hid);
762 
763 #if HAVE_HDF5_18
764  data_hid = H5Dopen (group_hid, "ridx", H5P_DEFAULT);
765 #else
766  data_hid = H5Dopen (group_hid, "ridx");
767 #endif
768  space_hid = H5Dget_space (data_hid);
769  rank = H5Sget_simple_extent_ndims (space_hid);
770 
771  if (rank != 2)
772  {
773  H5Sclose (space_hid);
774  H5Dclose (data_hid);
775  H5Gclose (group_hid);
776  return false;
777  }
778 
779  H5Sget_simple_extent_dims (space_hid, hdims, maxdims);
780 
781  if (static_cast<int> (hdims[0]) != nz
782  || static_cast<int> (hdims[1]) != 1)
783  {
784  H5Sclose (space_hid);
785  H5Dclose (data_hid);
786  H5Gclose (group_hid);
787  return false;
788  }
789 
790  itmp = m.xridx ();
791  if (H5Dread (data_hid, H5T_NATIVE_IDX, H5S_ALL,
792  H5S_ALL, H5P_DEFAULT, itmp) < 0)
793  {
794  H5Sclose (space_hid);
795  H5Dclose (data_hid);
796  H5Gclose (group_hid);
797  return false;
798  }
799 
800  H5Sclose (space_hid);
801  H5Dclose (data_hid);
802 
803 #if HAVE_HDF5_18
804  data_hid = H5Dopen (group_hid, "data", H5P_DEFAULT);
805 #else
806  data_hid = H5Dopen (group_hid, "data");
807 #endif
808  hid_t type_hid = H5Dget_type (data_hid);
809 
810  hid_t complex_type = hdf5_make_complex_type (H5T_NATIVE_DOUBLE);
811 
812  if (! hdf5_types_compatible (type_hid, complex_type))
813  {
814  H5Tclose (complex_type);
815  H5Dclose (data_hid);
816  H5Gclose (group_hid);
817  return false;
818  }
819 
820  space_hid = H5Dget_space (data_hid);
821  rank = H5Sget_simple_extent_ndims (space_hid);
822 
823  if (rank != 2)
824  {
825  H5Sclose (space_hid);
826  H5Dclose (data_hid);
827  H5Gclose (group_hid);
828  return false;
829  }
830 
831  H5Sget_simple_extent_dims (space_hid, hdims, maxdims);
832 
833  if (static_cast<int> (hdims[0]) != nz
834  || static_cast<int> (hdims[1]) != 1)
835  {
836  H5Sclose (space_hid);
837  H5Dclose (data_hid);
838  H5Gclose (group_hid);
839  return false;
840  }
841 
842  Complex *ctmp = m.xdata ();
843  bool retval = false;
844  if (H5Dread (data_hid, complex_type, H5S_ALL, H5S_ALL,
845  H5P_DEFAULT, ctmp) >= 0
846  && m.indices_ok ())
847  {
848  retval = true;
849  matrix = m;
850  }
851 
852  H5Tclose (complex_type);
853  H5Sclose (space_hid);
854  H5Dclose (data_hid);
855  H5Gclose (group_hid);
856 
857  return retval;
858 }
859 
860 #endif
861 
862 mxArray *
864 {
865  mwSize nz = nzmax ();
866  mxArray *retval = new mxArray (mxDOUBLE_CLASS, rows (), columns (),
867  nz, mxCOMPLEX);
868  double *pr = static_cast<double *> (retval->get_data ());
869  double *pi = static_cast<double *> (retval->get_imag_data ());
870  mwIndex *ir = retval->get_ir ();
871  mwIndex *jc = retval->get_jc ();
872 
873  for (mwIndex i = 0; i < nz; i++)
874  {
875  Complex val = matrix.data (i);
876  pr[i] = std::real (val);
877  pi[i] = std::imag (val);
878  ir[i] = matrix.ridx (i);
879  }
880 
881  for (mwIndex i = 0; i < columns () + 1; i++)
882  jc[i] = matrix.cidx (i);
883 
884  return retval;
885 }
886 
889 {
890  switch (umap)
891  {
892  // Mappers handled specially.
893  case umap_real:
895  case umap_imag:
897 
898 #define ARRAY_METHOD_MAPPER(UMAP, FCN) \
899  case umap_ ## UMAP: \
900  return octave_value (matrix.FCN ())
901 
903 
904 #define ARRAY_MAPPER(UMAP, TYPE, FCN) \
905  case umap_ ## UMAP: \
906  return octave_value (matrix.map<TYPE> (FCN))
907 
910  ARRAY_MAPPER (angle, double, std::arg);
911  ARRAY_MAPPER (arg, double, std::arg);
916  ARRAY_MAPPER (erf, Complex, ::erf);
922  ARRAY_MAPPER (conj, Complex, std::conj<double>);
923  ARRAY_MAPPER (cos, Complex, std::cos);
924  ARRAY_MAPPER (cosh, Complex, std::cosh);
925  ARRAY_MAPPER (exp, Complex, std::exp);
927  ARRAY_MAPPER (fix, Complex, ::fix);
929  ARRAY_MAPPER (log, Complex, std::log);
930  ARRAY_MAPPER (log2, Complex, xlog2);
931  ARRAY_MAPPER (log10, Complex, std::log10);
933  ARRAY_MAPPER (round, Complex, xround);
934  ARRAY_MAPPER (roundb, Complex, xroundb);
936  ARRAY_MAPPER (sin, Complex, std::sin);
937  ARRAY_MAPPER (sinh, Complex, std::sinh);
938  ARRAY_MAPPER (sqrt, Complex, std::sqrt);
939  ARRAY_MAPPER (tan, Complex, std::tan);
940  ARRAY_MAPPER (tanh, Complex, std::tanh);
941  ARRAY_MAPPER (isnan, bool, xisnan);
942  ARRAY_MAPPER (isna, bool, octave_is_NA);
943  ARRAY_MAPPER (isinf, bool, xisinf);
944  ARRAY_MAPPER (finite, bool, xfinite);
945 
946  default: // Attempt to go via dense matrix.
948  }
949 }