Mercurial > hg > octave-lyh

comparison src/octave-value/ov-cx-sparse.cc @ 15057:46b19589b593

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Partition src/ directory with more code in subdirs. Create new octave-value dir for ov* code. Create new parse-tre dir for pt* code. Move OPERATORS and TEMPLATE-INST directories to lowercase names * octave-value/module.mk: Hook code in octave-value dir into build system. * octave-value/ov-base-diag.cc, octave-value/ov-base-diag.h, octave-value/ov-base-int.cc, octave-value/ov-base-int.h, octave-value/ov-base-mat.cc, octave-value/ov-base-mat.h, octave-value/ov-base-scalar.cc, octave-value/ov-base-scalar.h, octave-value/ov-base-sparse.cc, octave-value/ov-base-sparse.h, octave-value/ov-base.cc, octave-value/ov-base.h, octave-value/ov-bool-mat.cc, octave-value/ov-bool-mat.h, octave-value/ov-bool-sparse.cc, octave-value/ov-bool-sparse.h, octave-value/ov-bool.cc, octave-value/ov-bool.h, octave-value/ov-builtin.cc, octave-value/ov-builtin.h, octave-value/ov-cell.cc, octave-value/ov-cell.h, octave-value/ov-ch-mat.cc, octave-value/ov-ch-mat.h, octave-value/ov-class.cc, octave-value/ov-class.h, octave-value/ov-colon.cc, octave-value/ov-colon.h, octave-value/ov-complex.cc, octave-value/ov-complex.h, octave-value/ov-cs-list.cc, octave-value/ov-cs-list.h, octave-value/ov-cx-diag.cc, octave-value/ov-cx-diag.h, octave-value/ov-cx-mat.cc, octave-value/ov-cx-mat.h, octave-value/ov-cx-sparse.cc, octave-value/ov-cx-sparse.h, octave-value/ov-dld-fcn.cc, octave-value/ov-dld-fcn.h, octave-value/ov-fcn-handle.cc, octave-value/ov-fcn-handle.h, octave-value/ov-fcn-inline.cc, octave-value/ov-fcn-inline.h, octave-value/ov-fcn.cc, octave-value/ov-fcn.h, octave-value/ov-float.cc, octave-value/ov-float.h, octave-value/ov-flt-complex.cc, octave-value/ov-flt-complex.h, octave-value/ov-flt-cx-diag.cc, octave-value/ov-flt-cx-diag.h, octave-value/ov-flt-cx-mat.cc, octave-value/ov-flt-cx-mat.h, octave-value/ov-flt-re-diag.cc, octave-value/ov-flt-re-diag.h, octave-value/ov-flt-re-mat.cc, octave-value/ov-flt-re-mat.h, octave-value/ov-int-traits.h, octave-value/ov-int16.cc, octave-value/ov-int16.h, octave-value/ov-int32.cc, octave-value/ov-int32.h, octave-value/ov-int64.cc, octave-value/ov-int64.h, octave-value/ov-int8.cc, octave-value/ov-int8.h, octave-value/ov-intx.h, octave-value/ov-lazy-idx.cc, octave-value/ov-lazy-idx.h, octave-value/ov-mex-fcn.cc, octave-value/ov-mex-fcn.h, octave-value/ov-null-mat.cc, octave-value/ov-null-mat.h, octave-value/ov-oncleanup.cc, octave-value/ov-oncleanup.h, octave-value/ov-perm.cc, octave-value/ov-perm.h, octave-value/ov-range.cc, octave-value/ov-range.h, octave-value/ov-re-diag.cc, octave-value/ov-re-diag.h, octave-value/ov-re-mat.cc, octave-value/ov-re-mat.h, octave-value/ov-re-sparse.cc, octave-value/ov-re-sparse.h, octave-value/ov-scalar.cc, octave-value/ov-scalar.h, octave-value/ov-str-mat.cc, octave-value/ov-str-mat.h, octave-value/ov-struct.cc, octave-value/ov-struct.h, octave-value/ov-type-conv.h, octave-value/ov-typeinfo.cc, octave-value/ov-typeinfo.h, octave-value/ov-uint16.cc, octave-value/ov-uint16.h, octave-value/ov-uint32.cc, octave-value/ov-uint32.h, octave-value/ov-uint64.cc, octave-value/ov-uint64.h, octave-value/ov-uint8.cc, octave-value/ov-uint8.h, octave-value/ov-usr-fcn.cc, octave-value/ov-usr-fcn.h, octave-value/ov.cc, octave-value/ov.h: Moved from src/ dir to octave-value dir. * operators/module.mk, operators/op-b-b.cc, operators/op-b-bm.cc, operators/op-b-sbm.cc, operators/op-bm-b.cc, operators/op-bm-bm.cc, operators/op-bm-sbm.cc, operators/op-cdm-cdm.cc, operators/op-cdm-cm.cc, operators/op-cdm-cs.cc, operators/op-cdm-dm.cc, operators/op-cdm-m.cc, operators/op-cdm-s.cc, operators/op-cell.cc, operators/op-chm.cc, operators/op-class.cc, operators/op-cm-cdm.cc, operators/op-cm-cm.cc, operators/op-cm-cs.cc, operators/op-cm-dm.cc, operators/op-cm-m.cc, operators/op-cm-pm.cc, operators/op-cm-s.cc, operators/op-cm-scm.cc, operators/op-cm-sm.cc, operators/op-cs-cm.cc, operators/op-cs-cs.cc, operators/op-cs-m.cc, operators/op-cs-s.cc, operators/op-cs-scm.cc, operators/op-cs-sm.cc, operators/op-dm-cdm.cc, operators/op-dm-cm.cc, operators/op-dm-cs.cc, operators/op-dm-dm.cc, operators/op-dm-m.cc, operators/op-dm-s.cc, operators/op-dm-scm.cc, operators/op-dm-sm.cc, operators/op-dm-template.cc, operators/op-dms-template.cc, operators/op-double-conv.cc, operators/op-fcdm-fcdm.cc, operators/op-fcdm-fcm.cc, operators/op-fcdm-fcs.cc, operators/op-fcdm-fdm.cc, operators/op-fcdm-fm.cc, operators/op-fcdm-fs.cc, operators/op-fcm-fcdm.cc, operators/op-fcm-fcm.cc, operators/op-fcm-fcs.cc, operators/op-fcm-fdm.cc, operators/op-fcm-fm.cc, operators/op-fcm-fs.cc, operators/op-fcm-pm.cc, operators/op-fcn.cc, operators/op-fcs-fcm.cc, operators/op-fcs-fcs.cc, operators/op-fcs-fm.cc, operators/op-fcs-fs.cc, operators/op-fdm-fcdm.cc, operators/op-fdm-fcm.cc, operators/op-fdm-fcs.cc, operators/op-fdm-fdm.cc, operators/op-fdm-fm.cc, operators/op-fdm-fs.cc, operators/op-float-conv.cc, operators/op-fm-fcdm.cc, operators/op-fm-fcm.cc, operators/op-fm-fcs.cc, operators/op-fm-fdm.cc, operators/op-fm-fm.cc, operators/op-fm-fs.cc, operators/op-fm-pm.cc, operators/op-fs-fcm.cc, operators/op-fs-fcs.cc, operators/op-fs-fm.cc, operators/op-fs-fs.cc, operators/op-i16-i16.cc, operators/op-i32-i32.cc, operators/op-i64-i64.cc, operators/op-i8-i8.cc, operators/op-int-concat.cc, operators/op-int-conv.cc, operators/op-int.h, operators/op-m-cdm.cc, operators/op-m-cm.cc, operators/op-m-cs.cc, operators/op-m-dm.cc, operators/op-m-m.cc, operators/op-m-pm.cc, operators/op-m-s.cc, operators/op-m-scm.cc, operators/op-m-sm.cc, operators/op-pm-cm.cc, operators/op-pm-fcm.cc, operators/op-pm-fm.cc, operators/op-pm-m.cc, operators/op-pm-pm.cc, operators/op-pm-scm.cc, operators/op-pm-sm.cc, operators/op-pm-template.cc, operators/op-range.cc, operators/op-s-cm.cc, operators/op-s-cs.cc, operators/op-s-m.cc, operators/op-s-s.cc, operators/op-s-scm.cc, operators/op-s-sm.cc, operators/op-sbm-b.cc, operators/op-sbm-bm.cc, operators/op-sbm-sbm.cc, operators/op-scm-cm.cc, operators/op-scm-cs.cc, operators/op-scm-m.cc, operators/op-scm-s.cc, operators/op-scm-scm.cc, operators/op-scm-sm.cc, operators/op-sm-cm.cc, operators/op-sm-cs.cc, operators/op-sm-m.cc, operators/op-sm-s.cc, operators/op-sm-scm.cc, operators/op-sm-sm.cc, operators/op-str-m.cc, operators/op-str-s.cc, operators/op-str-str.cc, operators/op-struct.cc, operators/op-ui16-ui16.cc, operators/op-ui32-ui32.cc, operators/op-ui64-ui64.cc, operators/op-ui8-ui8.cc: Moved from OPERATORS/ dir to operators/ directory. * mkops: Correctly print comment in generated file ops.cc that it is made by mkops. Change sed expression for OPERATORS/ to operators/. * parse-tree/module.mk: Hook code in parse-tree dir into build system. * parse-tree/pt-all.h, parse-tree/pt-arg-list.cc, parse-tree/pt-arg-list.h, parse-tree/pt-assign.cc, parse-tree/pt-assign.h, parse-tree/pt-binop.cc, parse-tree/pt-binop.h, parse-tree/pt-bp.cc, parse-tree/pt-bp.h, parse-tree/pt-cbinop.cc, parse-tree/pt-cbinop.h, parse-tree/pt-cell.cc, parse-tree/pt-cell.h, parse-tree/pt-check.cc, parse-tree/pt-check.h, parse-tree/pt-cmd.cc, parse-tree/pt-cmd.h, parse-tree/pt-colon.cc, parse-tree/pt-colon.h, parse-tree/pt-const.cc, parse-tree/pt-const.h, parse-tree/pt-decl.cc, parse-tree/pt-decl.h, parse-tree/pt-eval.cc, parse-tree/pt-eval.h, parse-tree/pt-except.cc, parse-tree/pt-except.h, parse-tree/pt-exp.cc, parse-tree/pt-exp.h, parse-tree/pt-fcn-handle.cc, parse-tree/pt-fcn-handle.h, parse-tree/pt-id.cc, parse-tree/pt-id.h, parse-tree/pt-idx.cc, parse-tree/pt-idx.h, parse-tree/pt-jump.cc, parse-tree/pt-jump.h, parse-tree/pt-loop.cc, parse-tree/pt-loop.h, parse-tree/pt-mat.cc, parse-tree/pt-mat.h, parse-tree/pt-misc.cc, parse-tree/pt-misc.h, parse-tree/pt-pr-code.cc, parse-tree/pt-pr-code.h, parse-tree/pt-select.cc, parse-tree/pt-select.h, parse-tree/pt-stmt.cc, parse-tree/pt-stmt.h, parse-tree/pt-unop.cc, parse-tree/pt-unop.h, parse-tree/pt-walk.h, parse-tree/pt.cc, parse-tree/pt.h: Moved from src/ dir to parse-tree dir. * template-inst/Array-jit.cc, template-inst/Array-os.cc, template-inst/Array-sym.cc, template-inst/Array-tc.cc, template-inst/module.mk: Moved from TEMPLATE-INST dir to template-inst/ directory. * src/Makefile.am: Add new directories to build system. * corefcn/module.mk: Use COREFCN_SRC with all capitals to indicate it is not an Automake special target.
author Rik <rik@octave.org>
date Mon, 30 Jul 2012 15:29:19 -0700
parents src/ov-cx-sparse.cc@f7afecdd87ef
children 62a35ae7d6a2
comparison
equal deleted inserted replaced
15056:bc32288f4a42 15057:46b19589b593
1 /*
2
3 Copyright (C) 2004-2012 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 <climits>
29
30 #include <iostream>
31 #include <vector>
32
33 #include "lo-specfun.h"
34 #include "lo-mappers.h"
35 #include "oct-locbuf.h"
36
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
50 template class OCTINTERP_API octave_base_sparse<SparseComplexMatrix>;
51
52 DEFINE_OCTAVE_ALLOCATOR (octave_sparse_complex_matrix);
53
54 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_sparse_complex_matrix, "sparse complex matrix", "double");
55
56 octave_base_value *
57 octave_sparse_complex_matrix::try_narrowing_conversion (void)
58 {
59 octave_base_value *retval = 0;
60
61 if (Vsparse_auto_mutate)
62 {
63 int nr = matrix.rows ();
64 int nc = matrix.cols ();
65
66 // Don't use numel, since it can overflow for very large matrices
67 // Note that for the tests on matrix size, they become approximative
68 // since they involves a cast to double to avoid issues of overflow
69 if (matrix.rows () == 1 && matrix.cols () == 1)
70 {
71 // Const copy of the matrix, so the right version of () operator used
72 const SparseComplexMatrix tmp (matrix);
73
74 Complex c = tmp (0, 0);
75
76 if (std::imag (c) == 0.0)
77 retval = new octave_scalar (std::real (c));
78 else
79 retval = new octave_complex (c);
80 }
81 else if (nr == 0 || nc == 0)
82 retval = new octave_matrix (Matrix (nr, nc));
83 else if (matrix.all_elements_are_real ())
84 if (matrix.cols () > 0 && matrix.rows () > 0
85 && (double (matrix.byte_size ()) > double (matrix.rows ())
86 * double (matrix.cols ()) * sizeof (double)))
87 retval = new octave_matrix (::real (matrix.matrix_value ()));
88 else
89 retval = new octave_sparse_matrix (::real (matrix));
90 else if (matrix.cols () > 0 && matrix.rows () > 0
91 && (double (matrix.byte_size ()) > double (matrix.rows ())
92 * double (matrix.cols ()) * sizeof (Complex)))
93 retval = new octave_complex_matrix (matrix.matrix_value ());
94 }
95 else
96 {
97 if (matrix.all_elements_are_real ())
98 retval = new octave_sparse_matrix (::real (matrix));
99 }
100
101 return retval;
102 }
103
104 double
105 octave_sparse_complex_matrix::double_value (bool force_conversion) const
106 {
107 double retval = lo_ieee_nan_value ();
108
109 if (! force_conversion)
110 gripe_implicit_conversion ("Octave:imag-to-real",
111 "complex sparse matrix", "real scalar");
112
113 // FIXME -- maybe this should be a function, valid_as_scalar()
114 if (numel () > 0)
115 {
116 if (numel () > 1)
117 gripe_implicit_conversion ("Octave:array-to-scalar",
118 "complex sparse matrix", "real scalar");
119
120 retval = std::real (matrix (0, 0));
121 }
122 else
123 gripe_invalid_conversion ("complex sparse matrix", "real scalar");
124
125 return retval;
126 }
127
128 Matrix
129 octave_sparse_complex_matrix::matrix_value (bool force_conversion) const
130 {
131 Matrix retval;
132
133 if (! force_conversion)
134 gripe_implicit_conversion ("Octave:imag-to-real",
135 "complex sparse matrix", "real matrix");
136
137 retval = ::real (matrix.matrix_value ());
138
139 return retval;
140 }
141
142 Complex
143 octave_sparse_complex_matrix::complex_value (bool) const
144 {
145 double tmp = lo_ieee_nan_value ();
146
147 Complex retval (tmp, tmp);
148
149 // FIXME -- maybe this should be a function, valid_as_scalar()
150 if (numel () > 0)
151 {
152 if (numel () > 1)
153 gripe_implicit_conversion ("Octave:array-to-scalar",
154 "complex sparse matrix", "real scalar");
155
156 retval = matrix (0, 0);
157 }
158 else
159 gripe_invalid_conversion ("complex sparse matrix", "real scalar");
160
161 return retval;
162 }
163
164 ComplexMatrix
165 octave_sparse_complex_matrix::complex_matrix_value (bool) const
166 {
167 return matrix.matrix_value ();
168 }
169
170 ComplexNDArray
171 octave_sparse_complex_matrix::complex_array_value (bool) const
172 {
173 return ComplexNDArray (matrix.matrix_value ());
174 }
175
176 charNDArray
177 octave_sparse_complex_matrix::char_array_value (bool frc_str_conv) const
178 {
179 charNDArray retval;
180
181 if (! frc_str_conv)
182 gripe_implicit_conversion ("Octave:num-to-str",
183 "sparse complex matrix", "string");
184 else
185 {
186 retval = charNDArray (dims (), 0);
187 octave_idx_type nc = matrix.cols ();
188 octave_idx_type nr = matrix.rows ();
189
190 for (octave_idx_type j = 0; j < nc; j++)
191 for (octave_idx_type i = matrix.cidx (j); i < matrix.cidx (j+1); i++)
192 retval(matrix.ridx (i) + nr * j) =
193 static_cast<char>(std::real (matrix.data (i)));
194 }
195
196 return retval;
197 }
198
199 SparseMatrix
200 octave_sparse_complex_matrix::sparse_matrix_value (bool force_conversion) const
201 {
202 SparseMatrix retval;
203
204 if (! force_conversion)
205 gripe_implicit_conversion ("Octave:imag-to-real",
206 "complex sparse matrix",
207 "real sparse matrix");
208
209 retval = ::real (matrix);
210
211 return retval;
212 }
213
214 SparseBoolMatrix
215 octave_sparse_complex_matrix::sparse_bool_matrix_value (bool warn) const
216 {
217 if (matrix.any_element_is_nan ())
218 gripe_nan_to_logical_conversion ();
219 else if (warn && (! matrix.all_elements_are_real ()
220 || real (matrix).any_element_not_one_or_zero ()))
221 gripe_logical_conversion ();
222
223 return mx_el_ne (matrix, Complex (0.0));
224 }
225
226 bool
227 octave_sparse_complex_matrix::save_binary (std::ostream& os,
228 bool&save_as_floats)
229 {
230 dim_vector d = this->dims ();
231 if (d.length () < 1)
232 return false;
233
234 // Ensure that additional memory is deallocated
235 matrix.maybe_compress ();
236
237 int nr = d(0);
238 int nc = d(1);
239 int nz = nnz ();
240
241 int32_t itmp;
242 // Use negative value for ndims to be consistent with other formats
243 itmp= -2;
244 os.write (reinterpret_cast<char *> (&itmp), 4);
245
246 itmp= nr;
247 os.write (reinterpret_cast<char *> (&itmp), 4);
248
249 itmp= nc;
250 os.write (reinterpret_cast<char *> (&itmp), 4);
251
252 itmp= nz;
253 os.write (reinterpret_cast<char *> (&itmp), 4);
254
255 save_type st = LS_DOUBLE;
256 if (save_as_floats)
257 {
258 if (matrix.too_large_for_float ())
259 {
260 warning ("save: some values too large to save as floats --");
261 warning ("save: saving as doubles instead");
262 }
263 else
264 st = LS_FLOAT;
265 }
266 else if (matrix.nnz () > 8192) // FIXME -- make this configurable.
267 {
268 double max_val, min_val;
269 if (matrix.all_integers (max_val, min_val))
270 st = get_save_type (max_val, min_val);
271 }
272
273 // add one to the printed indices to go from
274 // zero-based to one-based arrays
275 for (int i = 0; i < nc+1; i++)
276 {
277 octave_quit ();
278 itmp = matrix.cidx (i);
279 os.write (reinterpret_cast<char *> (&itmp), 4);
280 }
281
282 for (int i = 0; i < nz; i++)
283 {
284 octave_quit ();
285 itmp = matrix.ridx (i);
286 os.write (reinterpret_cast<char *> (&itmp), 4);
287 }
288
289 write_doubles (os, reinterpret_cast<const double *> (matrix.data ()), st, 2 * nz);
290
291 return true;
292 }
293
294 bool
295 octave_sparse_complex_matrix::load_binary (std::istream& is, bool swap,
296 oct_mach_info::float_format fmt)
297 {
298 int32_t nz, nc, nr, tmp;
299 char ctmp;
300
301 if (! is.read (reinterpret_cast<char *> (&tmp), 4))
302 return false;
303
304 if (swap)
305 swap_bytes<4> (&tmp);
306
307 if (tmp != -2) {
308 error ("load: only 2D sparse matrices are supported");
309 return false;
310 }
311
312 if (! is.read (reinterpret_cast<char *> (&nr), 4))
313 return false;
314 if (! is.read (reinterpret_cast<char *> (&nc), 4))
315 return false;
316 if (! is.read (reinterpret_cast<char *> (&nz), 4))
317 return false;
318
319 if (swap)
320 {
321 swap_bytes<4> (&nr);
322 swap_bytes<4> (&nc);
323 swap_bytes<4> (&nz);
324 }
325
326 SparseComplexMatrix m (static_cast<octave_idx_type> (nr),
327 static_cast<octave_idx_type> (nc),
328 static_cast<octave_idx_type> (nz));
329
330 for (int i = 0; i < nc+1; i++)
331 {
332 octave_quit ();
333 if (! is.read (reinterpret_cast<char *> (&tmp), 4))
334 return false;
335 if (swap)
336 swap_bytes<4> (&tmp);
337 m.cidx (i) = tmp;
338 }
339
340 for (int i = 0; i < nz; i++)
341 {
342 octave_quit ();
343 if (! is.read (reinterpret_cast<char *> (&tmp), 4))
344 return false;
345 if (swap)
346 swap_bytes<4> (&tmp);
347 m.ridx (i) = tmp;
348 }
349
350 if (! is.read (reinterpret_cast<char *> (&ctmp), 1))
351 return false;
352
353 read_doubles (is, reinterpret_cast<double *> (m.data ()),
354 static_cast<save_type> (ctmp), 2 * nz, swap, fmt);
355
356 if (error_state || ! is)
357 return false;
358
359 if (! m.indices_ok ())
360 return false;
361
362 matrix = m;
363
364 return true;
365 }
366
367 #if defined (HAVE_HDF5)
368
369 bool
370 octave_sparse_complex_matrix::save_hdf5 (hid_t loc_id, const char *name,
371 bool save_as_floats)
372 {
373 dim_vector dv = dims ();
374 int empty = save_hdf5_empty (loc_id, name, dv);
375 if (empty)
376 return (empty > 0);
377
378 // Ensure that additional memory is deallocated
379 matrix.maybe_compress ();
380
381 #if HAVE_HDF5_18
382 hid_t group_hid = H5Gcreate (loc_id, name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
383 #else
384 hid_t group_hid = H5Gcreate (loc_id, name, 0);
385 #endif
386 if (group_hid < 0)
387 return false;
388
389 hid_t space_hid = -1, data_hid = -1;
390 bool retval = true;
391 SparseComplexMatrix m = sparse_complex_matrix_value ();
392 octave_idx_type tmp;
393 hsize_t hdims[2];
394
395 space_hid = H5Screate_simple (0, hdims, 0);
396 if (space_hid < 0)
397 {
398 H5Gclose (group_hid);
399 return false;
400 }
401
402 #if HAVE_HDF5_18
403 data_hid = H5Dcreate (group_hid, "nr", H5T_NATIVE_IDX, space_hid,
404 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
405 #else
406 data_hid = H5Dcreate (group_hid, "nr", H5T_NATIVE_IDX, space_hid,
407 H5P_DEFAULT);
408 #endif
409 if (data_hid < 0)
410 {
411 H5Sclose (space_hid);
412 H5Gclose (group_hid);
413 return false;
414 }
415
416 tmp = m.rows ();
417 retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL,
418 H5P_DEFAULT, &tmp) >= 0;
419 H5Dclose (data_hid);
420 if (!retval)
421 {
422 H5Sclose (space_hid);
423 H5Gclose (group_hid);
424 return false;
425 }
426
427 #if HAVE_HDF5_18
428 data_hid = H5Dcreate (group_hid, "nc", H5T_NATIVE_IDX, space_hid,
429 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
430 #else
431 data_hid = H5Dcreate (group_hid, "nc", H5T_NATIVE_IDX, space_hid,
432 H5P_DEFAULT);
433 #endif
434 if (data_hid < 0)
435 {
436 H5Sclose (space_hid);
437 H5Gclose (group_hid);
438 return false;
439 }
440
441 tmp = m.cols ();
442 retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL,
443 H5P_DEFAULT, &tmp) >= 0;
444 H5Dclose (data_hid);
445 if (!retval)
446 {
447 H5Sclose (space_hid);
448 H5Gclose (group_hid);
449 return false;
450 }
451
452 #if HAVE_HDF5_18
453 data_hid = H5Dcreate (group_hid, "nz", H5T_NATIVE_IDX, space_hid,
454 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
455 #else
456 data_hid = H5Dcreate (group_hid, "nz", H5T_NATIVE_IDX, space_hid,
457 H5P_DEFAULT);
458 #endif
459 if (data_hid < 0)
460 {
461 H5Sclose (space_hid);
462 H5Gclose (group_hid);
463 return false;
464 }
465
466 tmp = m.nnz ();
467 retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL,
468 H5P_DEFAULT, &tmp) >= 0;
469 H5Dclose (data_hid);
470 if (!retval)
471 {
472 H5Sclose (space_hid);
473 H5Gclose (group_hid);
474 return false;
475 }
476
477 H5Sclose (space_hid);
478
479 hdims[0] = m.cols () + 1;
480 hdims[1] = 1;
481
482 space_hid = H5Screate_simple (2, hdims, 0);
483
484 if (space_hid < 0)
485 {
486 H5Gclose (group_hid);
487 return false;
488 }
489
490 #if HAVE_HDF5_18
491 data_hid = H5Dcreate (group_hid, "cidx", H5T_NATIVE_IDX, space_hid,
492 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
493 #else
494 data_hid = H5Dcreate (group_hid, "cidx", H5T_NATIVE_IDX, space_hid,
495 H5P_DEFAULT);
496 #endif
497 if (data_hid < 0)
498 {
499 H5Sclose (space_hid);
500 H5Gclose (group_hid);
501 return false;
502 }
503
504 octave_idx_type * itmp = m.xcidx ();
505 retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL,
506 H5P_DEFAULT, itmp) >= 0;
507 H5Dclose (data_hid);
508 if (!retval)
509 {
510 H5Sclose (space_hid);
511 H5Gclose (group_hid);
512 return false;
513 }
514
515 H5Sclose (space_hid);
516
517 hdims[0] = m.nnz ();
518 hdims[1] = 1;
519
520 space_hid = H5Screate_simple (2, hdims, 0);
521
522 if (space_hid < 0)
523 {
524 H5Gclose (group_hid);
525 return false;
526 }
527
528 #if HAVE_HDF5_18
529 data_hid = H5Dcreate (group_hid, "ridx", H5T_NATIVE_IDX, space_hid,
530 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
531 #else
532 data_hid = H5Dcreate (group_hid, "ridx", H5T_NATIVE_IDX, space_hid,
533 H5P_DEFAULT);
534 #endif
535 if (data_hid < 0)
536 {
537 H5Sclose (space_hid);
538 H5Gclose (group_hid);
539 return false;
540 }
541
542 itmp = m.xridx ();
543 retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL, H5P_DEFAULT, itmp) >= 0;
544 H5Dclose (data_hid);
545 if (!retval)
546 {
547 H5Sclose (space_hid);
548 H5Gclose (group_hid);
549 return false;
550 }
551
552 hid_t save_type_hid = H5T_NATIVE_DOUBLE;
553
554 if (save_as_floats)
555 {
556 if (m.too_large_for_float ())
557 {
558 warning ("save: some values too large to save as floats --");
559 warning ("save: saving as doubles instead");
560 }
561 else
562 save_type_hid = H5T_NATIVE_FLOAT;
563 }
564 #if HAVE_HDF5_INT2FLOAT_CONVERSIONS
565 // hdf5 currently doesn't support float/integer conversions
566 else
567 {
568 double max_val, min_val;
569
570 if (m.all_integers (max_val, min_val))
571 save_type_hid
572 = save_type_to_hdf5 (get_save_type (max_val, min_val));
573 }
574 #endif /* HAVE_HDF5_INT2FLOAT_CONVERSIONS */
575
576 hid_t type_hid = hdf5_make_complex_type (save_type_hid);
577 if (type_hid < 0)
578 {
579 H5Sclose (space_hid);
580 H5Gclose (group_hid);
581 return false;
582 }
583 #if HAVE_HDF5_18
584 data_hid = H5Dcreate (group_hid, "data", type_hid, space_hid,
585 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
586 #else
587 data_hid = H5Dcreate (group_hid, "data", type_hid, space_hid, H5P_DEFAULT);
588 #endif
589 if (data_hid < 0)
590 {
591 H5Sclose (space_hid);
592 H5Tclose (type_hid);
593 H5Gclose (group_hid);
594 return false;
595 }
596
597 hid_t complex_type_hid = hdf5_make_complex_type (H5T_NATIVE_DOUBLE);
598 retval = false;
599 if (complex_type_hid >= 0)
600 {
601 Complex * ctmp = m.xdata ();
602
603 retval = H5Dwrite (data_hid, complex_type_hid, H5S_ALL, H5S_ALL,
604 H5P_DEFAULT, ctmp) >= 0;
605 }
606
607 H5Dclose (data_hid);
608 H5Sclose (space_hid);
609 H5Tclose (type_hid);
610 H5Gclose (group_hid);
611
612 return retval;
613 }
614
615 bool
616 octave_sparse_complex_matrix::load_hdf5 (hid_t loc_id, const char *name)
617 {
618 octave_idx_type nr, nc, nz;
619 hid_t group_hid, data_hid, space_hid;
620 hsize_t rank;
621
622 dim_vector dv;
623 int empty = load_hdf5_empty (loc_id, name, dv);
624 if (empty > 0)
625 matrix.resize (dv);
626 if (empty)
627 return (empty > 0);
628
629 #if HAVE_HDF5_18
630 group_hid = H5Gopen (loc_id, name, H5P_DEFAULT);
631 #else
632 group_hid = H5Gopen (loc_id, name);
633 #endif
634 if (group_hid < 0 ) return false;
635
636 #if HAVE_HDF5_18
637 data_hid = H5Dopen (group_hid, "nr", H5P_DEFAULT);
638 #else
639 data_hid = H5Dopen (group_hid, "nr");
640 #endif
641 space_hid = H5Dget_space (data_hid);
642 rank = H5Sget_simple_extent_ndims (space_hid);
643
644 if (rank != 0)
645 {
646 H5Dclose (data_hid);
647 H5Gclose (group_hid);
648 return false;
649 }
650
651 if (H5Dread (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL, H5P_DEFAULT, &nr) < 0)
652 {
653 H5Dclose (data_hid);
654 H5Gclose (group_hid);
655 return false;
656 }
657
658 H5Dclose (data_hid);
659
660 #if HAVE_HDF5_18
661 data_hid = H5Dopen (group_hid, "nc", H5P_DEFAULT);
662 #else
663 data_hid = H5Dopen (group_hid, "nc");
664 #endif
665 space_hid = H5Dget_space (data_hid);
666 rank = H5Sget_simple_extent_ndims (space_hid);
667
668 if (rank != 0)
669 {
670 H5Dclose (data_hid);
671 H5Gclose (group_hid);
672 return false;
673 }
674
675 if (H5Dread (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL, H5P_DEFAULT, &nc) < 0)
676 {
677 H5Dclose (data_hid);
678 H5Gclose (group_hid);
679 return false;
680 }
681
682 H5Dclose (data_hid);
683
684 #if HAVE_HDF5_18
685 data_hid = H5Dopen (group_hid, "nz", H5P_DEFAULT);
686 #else
687 data_hid = H5Dopen (group_hid, "nz");
688 #endif
689 space_hid = H5Dget_space (data_hid);
690 rank = H5Sget_simple_extent_ndims (space_hid);
691
692 if (rank != 0)
693 {
694 H5Dclose (data_hid);
695 H5Gclose (group_hid);
696 return false;
697 }
698
699 if (H5Dread (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL, H5P_DEFAULT, &nz) < 0)
700 {
701 H5Dclose (data_hid);
702 H5Gclose (group_hid);
703 return false;
704 }
705
706 H5Dclose (data_hid);
707
708 SparseComplexMatrix m (static_cast<octave_idx_type> (nr),
709 static_cast<octave_idx_type> (nc),
710 static_cast<octave_idx_type> (nz));
711
712 #if HAVE_HDF5_18
713 data_hid = H5Dopen (group_hid, "cidx", H5P_DEFAULT);
714 #else
715 data_hid = H5Dopen (group_hid, "cidx");
716 #endif
717 space_hid = H5Dget_space (data_hid);
718 rank = H5Sget_simple_extent_ndims (space_hid);
719
720 if (rank != 2)
721 {
722 H5Sclose (space_hid);
723 H5Dclose (data_hid);
724 H5Gclose (group_hid);
725 return false;
726 }
727
728 OCTAVE_LOCAL_BUFFER (hsize_t, hdims, rank);
729 OCTAVE_LOCAL_BUFFER (hsize_t, maxdims, rank);
730
731 H5Sget_simple_extent_dims (space_hid, hdims, maxdims);
732
733 if (static_cast<int> (hdims[0]) != nc + 1
734 || static_cast<int> (hdims[1]) != 1)
735 {
736 H5Sclose (space_hid);
737 H5Dclose (data_hid);
738 H5Gclose (group_hid);
739 return false;
740 }
741
742 octave_idx_type *itmp = m.xcidx ();
743 if (H5Dread (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL, H5P_DEFAULT, itmp) < 0)
744 {
745 H5Sclose (space_hid);
746 H5Dclose (data_hid);
747 H5Gclose (group_hid);
748 return false;
749 }
750
751 H5Sclose (space_hid);
752 H5Dclose (data_hid);
753
754 #if HAVE_HDF5_18
755 data_hid = H5Dopen (group_hid, "ridx", H5P_DEFAULT);
756 #else
757 data_hid = H5Dopen (group_hid, "ridx");
758 #endif
759 space_hid = H5Dget_space (data_hid);
760 rank = H5Sget_simple_extent_ndims (space_hid);
761
762 if (rank != 2)
763 {
764 H5Sclose (space_hid);
765 H5Dclose (data_hid);
766 H5Gclose (group_hid);
767 return false;
768 }
769
770 H5Sget_simple_extent_dims (space_hid, hdims, maxdims);
771
772 if (static_cast<int> (hdims[0]) != nz
773 || static_cast<int> (hdims[1]) != 1)
774 {
775 H5Sclose (space_hid);
776 H5Dclose (data_hid);
777 H5Gclose (group_hid);
778 return false;
779 }
780
781 itmp = m.xridx ();
782 if (H5Dread (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL, H5P_DEFAULT, itmp) < 0)
783 {
784 H5Sclose (space_hid);
785 H5Dclose (data_hid);
786 H5Gclose (group_hid);
787 return false;
788 }
789
790 H5Sclose (space_hid);
791 H5Dclose (data_hid);
792
793 #if HAVE_HDF5_18
794 data_hid = H5Dopen (group_hid, "data", H5P_DEFAULT);
795 #else
796 data_hid = H5Dopen (group_hid, "data");
797 #endif
798 hid_t type_hid = H5Dget_type (data_hid);
799
800 hid_t complex_type = hdf5_make_complex_type (H5T_NATIVE_DOUBLE);
801
802 if (! hdf5_types_compatible (type_hid, complex_type))
803 {
804 H5Tclose (complex_type);
805 H5Dclose (data_hid);
806 H5Gclose (group_hid);
807 return false;
808 }
809
810 space_hid = H5Dget_space (data_hid);
811 rank = H5Sget_simple_extent_ndims (space_hid);
812
813 if (rank != 2)
814 {
815 H5Sclose (space_hid);
816 H5Dclose (data_hid);
817 H5Gclose (group_hid);
818 return false;
819 }
820
821 H5Sget_simple_extent_dims (space_hid, hdims, maxdims);
822
823 if (static_cast<int> (hdims[0]) != nz
824 || static_cast<int> (hdims[1]) != 1)
825 {
826 H5Sclose (space_hid);
827 H5Dclose (data_hid);
828 H5Gclose (group_hid);
829 return false;
830 }
831
832 Complex *ctmp = m.xdata ();
833 bool retval = false;
834 if (H5Dread (data_hid, complex_type, H5S_ALL, H5S_ALL,
835 H5P_DEFAULT, ctmp) >= 0
836 && m.indices_ok ())
837 {
838 retval = true;
839 matrix = m;
840 }
841
842 H5Tclose (complex_type);
843 H5Sclose (space_hid);
844 H5Dclose (data_hid);
845 H5Gclose (group_hid);
846
847 return retval;
848 }
849
850 #endif
851
852 mxArray *
853 octave_sparse_complex_matrix::as_mxArray (void) const
854 {
855 mwSize nz = nzmax ();
856 mxArray *retval = new mxArray (mxDOUBLE_CLASS, rows (), columns (),
857 nz, mxCOMPLEX);
858 double *pr = static_cast<double *> (retval->get_data ());
859 double *pi = static_cast<double *> (retval->get_imag_data ());
860 mwIndex *ir = retval->get_ir ();
861 mwIndex *jc = retval->get_jc ();
862
863 for (mwIndex i = 0; i < nz; i++)
864 {
865 Complex val = matrix.data (i);
866 pr[i] = std::real (val);
867 pi[i] = std::imag (val);
868 ir[i] = matrix.ridx (i);
869 }
870
871 for (mwIndex i = 0; i < columns () + 1; i++)
872 jc[i] = matrix.cidx (i);
873
874 return retval;
875 }
876
877 octave_value
878 octave_sparse_complex_matrix::map (unary_mapper_t umap) const
879 {
880 switch (umap)
881 {
882 // Mappers handled specially.
883 case umap_real:
884 return ::real (matrix);
885 case umap_imag:
886 return ::imag (matrix);
887
888 #define ARRAY_METHOD_MAPPER(UMAP, FCN) \
889 case umap_ ## UMAP: \
890 return octave_value (matrix.FCN ())
891
892 ARRAY_METHOD_MAPPER (abs, abs);
893
894 #define ARRAY_MAPPER(UMAP, TYPE, FCN) \
895 case umap_ ## UMAP: \
896 return octave_value (matrix.map<TYPE> (FCN))
897
898 ARRAY_MAPPER (acos, Complex, ::acos);
899 ARRAY_MAPPER (acosh, Complex, ::acosh);
900 ARRAY_MAPPER (angle, double, std::arg);
901 ARRAY_MAPPER (arg, double, std::arg);
902 ARRAY_MAPPER (asin, Complex, ::asin);
903 ARRAY_MAPPER (asinh, Complex, ::asinh);
904 ARRAY_MAPPER (atan, Complex, ::atan);
905 ARRAY_MAPPER (atanh, Complex, ::atanh);
906 ARRAY_MAPPER (ceil, Complex, ::ceil);
907 ARRAY_MAPPER (conj, Complex, std::conj<double>);
908 ARRAY_MAPPER (cos, Complex, std::cos);
909 ARRAY_MAPPER (cosh, Complex, std::cosh);
910 ARRAY_MAPPER (exp, Complex, std::exp);
911 ARRAY_MAPPER (expm1, Complex, ::expm1);
912 ARRAY_MAPPER (fix, Complex, ::fix);
913 ARRAY_MAPPER (floor, Complex, ::floor);
914 ARRAY_MAPPER (log, Complex, std::log);
915 ARRAY_MAPPER (log2, Complex, xlog2);
916 ARRAY_MAPPER (log10, Complex, std::log10);
917 ARRAY_MAPPER (log1p, Complex, ::log1p);
918 ARRAY_MAPPER (round, Complex, xround);
919 ARRAY_MAPPER (roundb, Complex, xroundb);
920 ARRAY_MAPPER (signum, Complex, ::signum);
921 ARRAY_MAPPER (sin, Complex, std::sin);
922 ARRAY_MAPPER (sinh, Complex, std::sinh);
923 ARRAY_MAPPER (sqrt, Complex, std::sqrt);
924 ARRAY_MAPPER (tan, Complex, std::tan);
925 ARRAY_MAPPER (tanh, Complex, std::tanh);
926 ARRAY_MAPPER (isnan, bool, xisnan);
927 ARRAY_MAPPER (isna, bool, octave_is_NA);
928 ARRAY_MAPPER (isinf, bool, xisinf);
929 ARRAY_MAPPER (finite, bool, xfinite);
930
931 default: // Attempt to go via dense matrix.
932 return octave_base_sparse<SparseComplexMatrix>::map (umap);
933 }
934 }