1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1998-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 <istream>
31 #include <limits>
32 #include <ostream>
33 #include <vector>
34
35 #include "lo-specfun.h"
36 #include "lo-mappers.h"
37 #include "oct-locbuf.h"
38
39 #include "mxarray.h"
40 #include "ov-base.h"
41 #include "ov-scalar.h"
42 #include "errwarn.h"
43
44 #include "oct-hdf5.h"
45 #include "ls-hdf5.h"
46
47 #include "ov-re-sparse.h"
48
49 #include "ov-base-sparse.h"
50 #include "ov-base-sparse.cc"
51
52 #include "ov-bool-sparse.h"
53
54
55 template class OCTINTERP_API octave_base_sparse<SparseMatrix>;
56
57 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_sparse_matrix, "sparse matrix",
58 "double");
59
60 idx_vector
index_vector(bool) const61 octave_sparse_matrix::index_vector (bool /* require_integers */) const
62 {
63 if (matrix.numel () == matrix.nnz ())
64 return idx_vector (array_value ());
65 else
66 {
67 std::string nm = '<' + type_name () + '>';
68 octave::err_invalid_index (nm.c_str ());
69 }
70 }
71
72 octave_base_value *
try_narrowing_conversion(void)73 octave_sparse_matrix::try_narrowing_conversion (void)
74 {
75 octave_base_value *retval = nullptr;
76
77 if (Vsparse_auto_mutate)
78 {
79 // Don't use numel, since it can overflow for very large matrices
80 // Note that for the second test, this means it becomes approximative
81 // since it involves a cast to double to avoid issues of overflow
82 if (matrix.rows () == 1 && matrix.cols () == 1)
83 {
84 // Const copy of the matrix, so the right version of () operator used
85 const SparseMatrix tmp (matrix);
86
87 retval = new octave_scalar (tmp (0));
88 }
89 else if (matrix.cols () > 0 && matrix.rows () > 0
90 && (double (matrix.byte_size ()) > double (matrix.rows ())
91 * double (matrix.cols ()) * sizeof (double)))
92 retval = new octave_matrix (matrix.matrix_value ());
93 }
94
95 return retval;
96 }
97
98 double
double_value(bool) const99 octave_sparse_matrix::double_value (bool) const
100 {
101 if (isempty ())
102 err_invalid_conversion ("real sparse matrix", "real scalar");
103
104 if (numel () > 1)
105 warn_implicit_conversion ("Octave:array-to-scalar",
106 "real sparse matrix", "real scalar");
107
108 return matrix(0, 0);
109 }
110
111 Complex
complex_value(bool) const112 octave_sparse_matrix::complex_value (bool) const
113 {
114 // FIXME: maybe this should be a function, valid_as_scalar()
115 if (rows () == 0 || columns () == 0)
116 err_invalid_conversion ("real sparse matrix", "complex scalar");
117
118 if (numel () > 1)
119 warn_implicit_conversion ("Octave:array-to-scalar",
120 "real sparse matrix", "complex scalar");
121
122 return Complex (matrix(0, 0), 0);
123 }
124
125 Matrix
matrix_value(bool) const126 octave_sparse_matrix::matrix_value (bool) const
127 {
128 return matrix.matrix_value ();
129 }
130
131 boolNDArray
bool_array_value(bool warn) const132 octave_sparse_matrix::bool_array_value (bool warn) const
133 {
134 NDArray m = matrix.matrix_value ();
135
136 if (m.any_element_is_nan ())
137 octave::err_nan_to_logical_conversion ();
138 if (warn && m.any_element_not_one_or_zero ())
139 warn_logical_conversion ();
140
141 return boolNDArray (m);
142 }
143
144 charNDArray
char_array_value(bool) const145 octave_sparse_matrix::char_array_value (bool) const
146 {
147 charNDArray retval (dims (), 0);
148 octave_idx_type nc = matrix.cols ();
149 octave_idx_type nr = matrix.rows ();
150
151 for (octave_idx_type j = 0; j < nc; j++)
152 for (octave_idx_type i = matrix.cidx (j); i < matrix.cidx (j+1); i++)
153 retval(matrix.ridx (i) + nr * j) = static_cast<char> (matrix.data (i));
154
155 return retval;
156 }
157
158 ComplexMatrix
complex_matrix_value(bool) const159 octave_sparse_matrix::complex_matrix_value (bool) const
160 {
161 return ComplexMatrix (matrix.matrix_value ());
162 }
163
164 ComplexNDArray
complex_array_value(bool) const165 octave_sparse_matrix::complex_array_value (bool) const
166 {
167 return ComplexNDArray (ComplexMatrix (matrix.matrix_value ()));
168 }
169
170 NDArray
array_value(bool) const171 octave_sparse_matrix::array_value (bool) const
172 {
173 return NDArray (matrix.matrix_value ());
174 }
175
176 SparseBoolMatrix
sparse_bool_matrix_value(bool warn) const177 octave_sparse_matrix::sparse_bool_matrix_value (bool warn) const
178 {
179 if (matrix.any_element_is_nan ())
180 octave::err_nan_to_logical_conversion ();
181 if (warn && matrix.any_element_not_one_or_zero ())
182 warn_logical_conversion ();
183
184 return mx_el_ne (matrix, 0.0);
185 }
186
187 octave_value
convert_to_str_internal(bool,bool,char type) const188 octave_sparse_matrix::convert_to_str_internal (bool, bool, char type) const
189 {
190 octave_value retval;
191 dim_vector dv = dims ();
192 octave_idx_type nel = dv.numel ();
193
194 if (nel == 0)
195 {
196 char s = '\0';
197 retval = octave_value (&s, type);
198 }
199 else
200 {
201 octave_idx_type nr = matrix.rows ();
202 octave_idx_type nc = matrix.cols ();
203 charNDArray chm (dv, static_cast<char> (0));
204
205 bool warned = false;
206
207 for (octave_idx_type j = 0; j < nc; j++)
208 for (octave_idx_type i = matrix.cidx (j);
209 i < matrix.cidx (j+1); i++)
210 {
211 octave_quit ();
212
213 double d = matrix.data (i);
214
215 if (octave::math::isnan (d))
216 octave::err_nan_to_character_conversion ();
217
218 int ival = octave::math::nint (d);
219
220 if (ival < 0 || ival > std::numeric_limits<unsigned char>::max ())
221 {
222 // FIXME: is there something better we could do?
223
224 ival = 0;
225
226 if (! warned)
227 {
228 ::warning ("range error for conversion to character value");
229 warned = true;
230 }
231 }
232
233 chm(matrix.ridx (i) + j * nr) = static_cast<char> (ival);
234 }
235
236 retval = octave_value (chm, type);
237 }
238
239 return retval;
240 }
241
242 octave_value
as_double(void) const243 octave_sparse_matrix::as_double (void) const
244 {
245 return this->matrix;
246 }
247
248 bool
save_binary(std::ostream & os,bool save_as_floats)249 octave_sparse_matrix::save_binary (std::ostream& os, bool save_as_floats)
250 {
251 dim_vector dv = this->dims ();
252 if (dv.ndims () < 1)
253 return false;
254
255 // Ensure that additional memory is deallocated
256 matrix.maybe_compress ();
257
258 int nr = dv(0);
259 int nc = dv(1);
260 int nz = nnz ();
261
262 int32_t itmp;
263 // Use negative value for ndims to be consistent with other formats
264 itmp = -2;
265 os.write (reinterpret_cast<char *> (&itmp), 4);
266
267 itmp = nr;
268 os.write (reinterpret_cast<char *> (&itmp), 4);
269
270 itmp = nc;
271 os.write (reinterpret_cast<char *> (&itmp), 4);
272
273 itmp = nz;
274 os.write (reinterpret_cast<char *> (&itmp), 4);
275
276 save_type st = LS_DOUBLE;
277 if (save_as_floats)
278 {
279 if (matrix.too_large_for_float ())
280 {
281 warning ("save: some values too large to save as floats --");
282 warning ("save: saving as doubles instead");
283 }
284 else
285 st = LS_FLOAT;
286 }
287 else if (matrix.nnz () > 8192) // FIXME: make this configurable.
288 {
289 double max_val, min_val;
290 if (matrix.all_integers (max_val, min_val))
291 st = get_save_type (max_val, min_val);
292 }
293
294 // add one to the printed indices to go from
295 // zero-based to one-based arrays
296 for (int i = 0; i < nc+1; i++)
297 {
298 octave_quit ();
299 itmp = matrix.cidx (i);
300 os.write (reinterpret_cast<char *> (&itmp), 4);
301 }
302
303 for (int i = 0; i < nz; i++)
304 {
305 octave_quit ();
306 itmp = matrix.ridx (i);
307 os.write (reinterpret_cast<char *> (&itmp), 4);
308 }
309
310 write_doubles (os, matrix.data (), st, nz);
311
312 return true;
313 }
314
315 bool
load_binary(std::istream & is,bool swap,octave::mach_info::float_format fmt)316 octave_sparse_matrix::load_binary (std::istream& is, bool swap,
317 octave::mach_info::float_format fmt)
318 {
319 int32_t nz, nc, nr, tmp;
320 char ctmp;
321
322 if (! is.read (reinterpret_cast<char *> (&tmp), 4))
323 return false;
324
325 if (swap)
326 swap_bytes<4> (&tmp);
327
328 if (tmp != -2)
329 error ("load: only 2-D sparse matrices are supported");
330
331 if (! is.read (reinterpret_cast<char *> (&nr), 4))
332 return false;
333 if (! is.read (reinterpret_cast<char *> (&nc), 4))
334 return false;
335 if (! is.read (reinterpret_cast<char *> (&nz), 4))
336 return false;
337
338 if (swap)
339 {
340 swap_bytes<4> (&nr);
341 swap_bytes<4> (&nc);
342 swap_bytes<4> (&nz);
343 }
344
345 SparseMatrix m (static_cast<octave_idx_type> (nr),
346 static_cast<octave_idx_type> (nc),
347 static_cast<octave_idx_type> (nz));
348
349 for (int i = 0; i < nc+1; i++)
350 {
351 octave_quit ();
352 if (! is.read (reinterpret_cast<char *> (&tmp), 4))
353 return false;
354 if (swap)
355 swap_bytes<4> (&tmp);
356 m.xcidx (i) = tmp;
357 }
358
359 for (int i = 0; i < nz; i++)
360 {
361 octave_quit ();
362 if (! is.read (reinterpret_cast<char *> (&tmp), 4))
363 return false;
364 if (swap)
365 swap_bytes<4> (&tmp);
366 m.xridx (i) = tmp;
367 }
368
369 if (! is.read (reinterpret_cast<char *> (&ctmp), 1))
370 return false;
371
372 read_doubles (is, m.xdata (), static_cast<save_type> (ctmp), nz, swap, fmt);
373
374 if (! is)
375 return false;
376
377 if (! m.indices_ok ())
378 return false;
379
380 matrix = m;
381
382 return true;
383 }
384
385 bool
save_hdf5(octave_hdf5_id loc_id,const char * name,bool save_as_floats)386 octave_sparse_matrix::save_hdf5 (octave_hdf5_id loc_id, const char *name,
387 bool save_as_floats)
388 {
389 bool retval = false;
390
391 #if defined (HAVE_HDF5)
392
393 dim_vector dv = dims ();
394 int empty = save_hdf5_empty (loc_id, name, dv);
395 if (empty)
396 return (empty > 0);
397
398 // Ensure that additional memory is deallocated
399 matrix.maybe_compress ();
400
401 #if defined (HAVE_HDF5_18)
402 hid_t group_hid = H5Gcreate (loc_id, name, octave_H5P_DEFAULT,
403 octave_H5P_DEFAULT, octave_H5P_DEFAULT);
404 #else
405 hid_t group_hid = H5Gcreate (loc_id, name, 0);
406 #endif
407 if (group_hid < 0)
408 return false;
409
410 hid_t space_hid, data_hid;
411 space_hid = data_hid = -1;
412 SparseMatrix m = sparse_matrix_value ();
413 octave_idx_type tmp;
414 hsize_t hdims[2];
415
416 space_hid = H5Screate_simple (0, hdims, nullptr);
417 if (space_hid < 0)
418 {
419 H5Gclose (group_hid);
420 return false;
421 }
422 #if defined (HAVE_HDF5_18)
423 data_hid = H5Dcreate (group_hid, "nr", H5T_NATIVE_IDX, space_hid,
424 octave_H5P_DEFAULT, octave_H5P_DEFAULT, octave_H5P_DEFAULT);
425 #else
426 data_hid = H5Dcreate (group_hid, "nr", H5T_NATIVE_IDX, space_hid,
427 octave_H5P_DEFAULT);
428 #endif
429 if (data_hid < 0)
430 {
431 H5Sclose (space_hid);
432 H5Gclose (group_hid);
433 return false;
434 }
435
436 tmp = m.rows ();
437 retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, octave_H5S_ALL, octave_H5S_ALL,
438 octave_H5P_DEFAULT, &tmp) >= 0;
439 H5Dclose (data_hid);
440 if (! retval)
441 {
442 H5Sclose (space_hid);
443 H5Gclose (group_hid);
444 return false;
445 }
446 #if defined (HAVE_HDF5_18)
447 data_hid = H5Dcreate (group_hid, "nc", H5T_NATIVE_IDX, space_hid,
448 octave_H5P_DEFAULT, octave_H5P_DEFAULT, octave_H5P_DEFAULT);
449 #else
450 data_hid = H5Dcreate (group_hid, "nc", H5T_NATIVE_IDX, space_hid,
451 octave_H5P_DEFAULT);
452 #endif
453 if (data_hid < 0)
454 {
455 H5Sclose (space_hid);
456 H5Gclose (group_hid);
457 return false;
458 }
459
460 tmp = m.cols ();
461 retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, octave_H5S_ALL, octave_H5S_ALL,
462 octave_H5P_DEFAULT, &tmp) >= 0;
463 H5Dclose (data_hid);
464 if (! retval)
465 {
466 H5Sclose (space_hid);
467 H5Gclose (group_hid);
468 return false;
469 }
470
471 #if defined (HAVE_HDF5_18)
472 data_hid = H5Dcreate (group_hid, "nz", H5T_NATIVE_IDX, space_hid,
473 octave_H5P_DEFAULT, octave_H5P_DEFAULT, octave_H5P_DEFAULT);
474 #else
475 data_hid = H5Dcreate (group_hid, "nz", H5T_NATIVE_IDX, space_hid,
476 octave_H5P_DEFAULT);
477 #endif
478 if (data_hid < 0)
479 {
480 H5Sclose (space_hid);
481 H5Gclose (group_hid);
482 return false;
483 }
484
485 tmp = m.nnz ();
486 retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, octave_H5S_ALL, octave_H5S_ALL,
487 octave_H5P_DEFAULT, &tmp) >= 0;
488 H5Dclose (data_hid);
489 if (! retval)
490 {
491 H5Sclose (space_hid);
492 H5Gclose (group_hid);
493 return false;
494 }
495
496 H5Sclose (space_hid);
497
498 hdims[0] = m.cols () + 1;
499 hdims[1] = 1;
500
501 space_hid = H5Screate_simple (2, hdims, nullptr);
502
503 if (space_hid < 0)
504 {
505 H5Gclose (group_hid);
506 return false;
507 }
508
509 #if defined (HAVE_HDF5_18)
510 data_hid = H5Dcreate (group_hid, "cidx", H5T_NATIVE_IDX, space_hid,
511 octave_H5P_DEFAULT, octave_H5P_DEFAULT, octave_H5P_DEFAULT);
512 #else
513 data_hid = H5Dcreate (group_hid, "cidx", H5T_NATIVE_IDX, space_hid,
514 octave_H5P_DEFAULT);
515 #endif
516 if (data_hid < 0)
517 {
518 H5Sclose (space_hid);
519 H5Gclose (group_hid);
520 return false;
521 }
522
523 octave_idx_type *itmp = m.xcidx ();
524 retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, octave_H5S_ALL, octave_H5S_ALL,
525 octave_H5P_DEFAULT, itmp) >= 0;
526 H5Dclose (data_hid);
527 if (! retval)
528 {
529 H5Sclose (space_hid);
530 H5Gclose (group_hid);
531 return false;
532 }
533
534 H5Sclose (space_hid);
535
536 hdims[0] = m.nnz ();
537 hdims[1] = 1;
538
539 space_hid = H5Screate_simple (2, hdims, nullptr);
540
541 if (space_hid < 0)
542 {
543 H5Gclose (group_hid);
544 return false;
545 }
546 #if defined (HAVE_HDF5_18)
547 data_hid = H5Dcreate (group_hid, "ridx", H5T_NATIVE_IDX, space_hid,
548 octave_H5P_DEFAULT, octave_H5P_DEFAULT, octave_H5P_DEFAULT);
549 #else
550 data_hid = H5Dcreate (group_hid, "ridx", H5T_NATIVE_IDX, space_hid,
551 octave_H5P_DEFAULT);
552 #endif
553 if (data_hid < 0)
554 {
555 H5Sclose (space_hid);
556 H5Gclose (group_hid);
557 return false;
558 }
559
560 itmp = m.xridx ();
561 retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, octave_H5S_ALL, octave_H5S_ALL,
562 octave_H5P_DEFAULT, itmp) >= 0;
563 H5Dclose (data_hid);
564 if (! retval)
565 {
566 H5Sclose (space_hid);
567 H5Gclose (group_hid);
568 return false;
569 }
570
571 hid_t save_type_hid = H5T_NATIVE_DOUBLE;
572
573 if (save_as_floats)
574 {
575 if (m.too_large_for_float ())
576 {
577 warning ("save: some values too large to save as floats --");
578 warning ("save: saving as doubles instead");
579 }
580 else
581 save_type_hid = H5T_NATIVE_FLOAT;
582 }
583 #if defined (HAVE_HDF5_INT2FLOAT_CONVERSIONS)
584 // hdf5 currently doesn't support float/integer conversions
585 else
586 {
587 double max_val, min_val;
588
589 if (m.all_integers (max_val, min_val))
590 save_type_hid
591 = save_type_to_hdf5 (get_save_type (max_val, min_val));
592 }
593 #endif
594
595 #if defined (HAVE_HDF5_18)
596 data_hid = H5Dcreate (group_hid, "data", save_type_hid, space_hid,
597 octave_H5P_DEFAULT, octave_H5P_DEFAULT, octave_H5P_DEFAULT);
598 #else
599 data_hid = H5Dcreate (group_hid, "data", save_type_hid, space_hid,
600 octave_H5P_DEFAULT);
601 #endif
602 if (data_hid < 0)
603 {
604 H5Sclose (space_hid);
605 H5Gclose (group_hid);
606 return false;
607 }
608
609 double *dtmp = m.xdata ();
610 retval = H5Dwrite (data_hid, H5T_NATIVE_DOUBLE, octave_H5S_ALL, octave_H5S_ALL,
611 octave_H5P_DEFAULT, dtmp) >= 0;
612 H5Dclose (data_hid);
613 H5Sclose (space_hid);
614 H5Gclose (group_hid);
615
616 #else
617 octave_unused_parameter (loc_id);
618 octave_unused_parameter (name);
619 octave_unused_parameter (save_as_floats);
620
621 warn_save ("hdf5");
622 #endif
623
624 return retval;
625 }
626
627 bool
load_hdf5(octave_hdf5_id loc_id,const char * name)628 octave_sparse_matrix::load_hdf5 (octave_hdf5_id loc_id, const char *name)
629 {
630 bool retval = false;
631
632 #if defined (HAVE_HDF5)
633
634 octave_idx_type nr, nc, nz;
635 hid_t group_hid, data_hid, space_hid;
636 hsize_t rank;
637
638 dim_vector dv;
639 int empty = load_hdf5_empty (loc_id, name, dv);
640 if (empty > 0)
641 matrix.resize (dv);
642 if (empty)
643 return (empty > 0);
644
645 #if defined (HAVE_HDF5_18)
646 group_hid = H5Gopen (loc_id, name, octave_H5P_DEFAULT);
647 #else
648 group_hid = H5Gopen (loc_id, name);
649 #endif
650 if (group_hid < 0) return false;
651
652 #if defined (HAVE_HDF5_18)
653 data_hid = H5Dopen (group_hid, "nr", octave_H5P_DEFAULT);
654 #else
655 data_hid = H5Dopen (group_hid, "nr");
656 #endif
657 space_hid = H5Dget_space (data_hid);
658 rank = H5Sget_simple_extent_ndims (space_hid);
659
660 if (rank != 0)
661 {
662 H5Dclose (data_hid);
663 H5Gclose (group_hid);
664 return false;
665 }
666
667 if (H5Dread (data_hid, H5T_NATIVE_IDX, octave_H5S_ALL, octave_H5S_ALL,
668 octave_H5P_DEFAULT, &nr) < 0)
669 {
670 H5Dclose (data_hid);
671 H5Gclose (group_hid);
672 return false;
673 }
674
675 H5Dclose (data_hid);
676
677 #if defined (HAVE_HDF5_18)
678 data_hid = H5Dopen (group_hid, "nc", octave_H5P_DEFAULT);
679 #else
680 data_hid = H5Dopen (group_hid, "nc");
681 #endif
682 space_hid = H5Dget_space (data_hid);
683 rank = H5Sget_simple_extent_ndims (space_hid);
684
685 if (rank != 0)
686 {
687 H5Dclose (data_hid);
688 H5Gclose (group_hid);
689 return false;
690 }
691
692 if (H5Dread (data_hid, H5T_NATIVE_IDX, octave_H5S_ALL, octave_H5S_ALL,
693 octave_H5P_DEFAULT, &nc) < 0)
694 {
695 H5Dclose (data_hid);
696 H5Gclose (group_hid);
697 return false;
698 }
699
700 H5Dclose (data_hid);
701
702 #if defined (HAVE_HDF5_18)
703 data_hid = H5Dopen (group_hid, "nz", octave_H5P_DEFAULT);
704 #else
705 data_hid = H5Dopen (group_hid, "nz");
706 #endif
707 space_hid = H5Dget_space (data_hid);
708 rank = H5Sget_simple_extent_ndims (space_hid);
709
710 if (rank != 0)
711 {
712 H5Dclose (data_hid);
713 H5Gclose (group_hid);
714 return false;
715 }
716
717 if (H5Dread (data_hid, H5T_NATIVE_IDX, octave_H5S_ALL, octave_H5S_ALL,
718 octave_H5P_DEFAULT, &nz) < 0)
719 {
720 H5Dclose (data_hid);
721 H5Gclose (group_hid);
722 return false;
723 }
724
725 H5Dclose (data_hid);
726
727 SparseMatrix m (static_cast<octave_idx_type> (nr),
728 static_cast<octave_idx_type> (nc),
729 static_cast<octave_idx_type> (nz));
730
731 #if defined (HAVE_HDF5_18)
732 data_hid = H5Dopen (group_hid, "cidx", octave_H5P_DEFAULT);
733 #else
734 data_hid = H5Dopen (group_hid, "cidx");
735 #endif
736 space_hid = H5Dget_space (data_hid);
737 rank = H5Sget_simple_extent_ndims (space_hid);
738
739 if (rank != 2)
740 {
741 H5Sclose (space_hid);
742 H5Dclose (data_hid);
743 H5Gclose (group_hid);
744 return false;
745 }
746
747 OCTAVE_LOCAL_BUFFER (hsize_t, hdims, rank);
748 OCTAVE_LOCAL_BUFFER (hsize_t, maxdims, rank);
749
750 H5Sget_simple_extent_dims (space_hid, hdims, maxdims);
751
752 if (static_cast<int> (hdims[0]) != nc + 1
753 || static_cast<int> (hdims[1]) != 1)
754 {
755 H5Sclose (space_hid);
756 H5Dclose (data_hid);
757 H5Gclose (group_hid);
758 return false;
759 }
760
761 octave_idx_type *itmp = m.xcidx ();
762 if (H5Dread (data_hid, H5T_NATIVE_IDX, octave_H5S_ALL, octave_H5S_ALL,
763 octave_H5P_DEFAULT, itmp) < 0)
764 {
765 H5Sclose (space_hid);
766 H5Dclose (data_hid);
767 H5Gclose (group_hid);
768 return false;
769 }
770
771 H5Sclose (space_hid);
772 H5Dclose (data_hid);
773
774 #if defined (HAVE_HDF5_18)
775 data_hid = H5Dopen (group_hid, "ridx", octave_H5P_DEFAULT);
776 #else
777 data_hid = H5Dopen (group_hid, "ridx");
778 #endif
779 space_hid = H5Dget_space (data_hid);
780 rank = H5Sget_simple_extent_ndims (space_hid);
781
782 if (rank != 2)
783 {
784 H5Sclose (space_hid);
785 H5Dclose (data_hid);
786 H5Gclose (group_hid);
787 return false;
788 }
789
790 H5Sget_simple_extent_dims (space_hid, hdims, maxdims);
791
792 if (static_cast<int> (hdims[0]) != nz || static_cast<int> (hdims[1]) != 1)
793 {
794 H5Sclose (space_hid);
795 H5Dclose (data_hid);
796 H5Gclose (group_hid);
797 return false;
798 }
799
800 itmp = m.xridx ();
801 if (H5Dread (data_hid, H5T_NATIVE_IDX, octave_H5S_ALL, octave_H5S_ALL,
802 octave_H5P_DEFAULT, itmp) < 0)
803 {
804 H5Sclose (space_hid);
805 H5Dclose (data_hid);
806 H5Gclose (group_hid);
807 return false;
808 }
809
810 H5Sclose (space_hid);
811 H5Dclose (data_hid);
812
813 #if defined (HAVE_HDF5_18)
814 data_hid = H5Dopen (group_hid, "data", octave_H5P_DEFAULT);
815 #else
816 data_hid = H5Dopen (group_hid, "data");
817 #endif
818 space_hid = H5Dget_space (data_hid);
819 rank = H5Sget_simple_extent_ndims (space_hid);
820
821 if (rank != 2)
822 {
823 H5Sclose (space_hid);
824 H5Dclose (data_hid);
825 H5Gclose (group_hid);
826 return false;
827 }
828
829 H5Sget_simple_extent_dims (space_hid, hdims, maxdims);
830
831 if (static_cast<int> (hdims[0]) != nz || static_cast<int> (hdims[1]) != 1)
832 {
833 H5Sclose (space_hid);
834 H5Dclose (data_hid);
835 H5Gclose (group_hid);
836 return false;
837 }
838
839 double *dtmp = m.xdata ();
840
841 if (H5Dread (data_hid, H5T_NATIVE_DOUBLE, octave_H5S_ALL, octave_H5S_ALL,
842 octave_H5P_DEFAULT, dtmp) >= 0
843 && m.indices_ok ())
844 {
845 retval = true;
846 matrix = m;
847 }
848
849 H5Sclose (space_hid);
850 H5Dclose (data_hid);
851 H5Gclose (group_hid);
852
853 #else
854 octave_unused_parameter (loc_id);
855 octave_unused_parameter (name);
856
857 warn_load ("hdf5");
858 #endif
859
860 return retval;
861 }
862
863 mxArray *
as_mxArray(void) const864 octave_sparse_matrix::as_mxArray (void) const
865 {
866 mwSize nz = nzmax ();
867 mwSize nr = rows ();
868 mwSize nc = columns ();
869 mxArray *retval = new mxArray (mxDOUBLE_CLASS, nr, nc, nz, mxREAL);
870 double *pr = static_cast<double *> (retval->get_data ());
871 mwIndex *ir = retval->get_ir ();
872 mwIndex *jc = retval->get_jc ();
873
874 for (mwIndex i = 0; i < nz; i++)
875 {
876 pr[i] = matrix.data (i);
877 ir[i] = matrix.ridx (i);
878 }
879
880 for (mwIndex i = 0; i < nc + 1; i++)
881 jc[i] = matrix.cidx (i);
882
883 return retval;
884 }
885
886 octave_value
map(unary_mapper_t umap) const887 octave_sparse_matrix::map (unary_mapper_t umap) const
888 {
889 switch (umap)
890 {
891 case umap_imag:
892 return SparseMatrix (matrix.rows (), matrix.cols (), 0.0);
893
894 case umap_real:
895 case umap_conj:
896 return matrix;
897
898 // Mappers handled specially.
899 #define ARRAY_METHOD_MAPPER(UMAP, FCN) \
900 case umap_ ## UMAP: \
901 return octave_value (matrix.FCN ())
902
903 ARRAY_METHOD_MAPPER (abs, abs);
904
905 #define ARRAY_MAPPER(UMAP, TYPE, FCN) \
906 case umap_ ## UMAP: \
907 return octave_value (matrix.map<TYPE> (FCN))
908
909 ARRAY_MAPPER (acos, Complex, octave::math::rc_acos);
910 ARRAY_MAPPER (acosh, Complex, octave::math::rc_acosh);
911 ARRAY_MAPPER (angle, double, std::arg);
912 ARRAY_MAPPER (arg, double,std::arg);
913 ARRAY_MAPPER (asin, Complex, octave::math::rc_asin);
914 ARRAY_MAPPER (asinh, double, octave::math::asinh);
915 ARRAY_MAPPER (atan, double, ::atan);
916 ARRAY_MAPPER (atanh, Complex, octave::math::rc_atanh);
917 ARRAY_MAPPER (erf, double, octave::math::erf);
918 ARRAY_MAPPER (erfinv, double, octave::math::erfinv);
919 ARRAY_MAPPER (erfcinv, double, octave::math::erfcinv);
920 ARRAY_MAPPER (erfc, double, octave::math::erfc);
921 ARRAY_MAPPER (erfcx, double, octave::math::erfcx);
922 ARRAY_MAPPER (erfi, double, octave::math::erfi);
923 ARRAY_MAPPER (dawson, double, octave::math::dawson);
924 ARRAY_MAPPER (gamma, double, octave::math::gamma);
925 ARRAY_MAPPER (lgamma, Complex, octave::math::rc_lgamma);
926 ARRAY_MAPPER (cbrt, double, octave::math::cbrt);
927 ARRAY_MAPPER (ceil, double, ::ceil);
928 ARRAY_MAPPER (cos, double, ::cos);
929 ARRAY_MAPPER (cosh, double, ::cosh);
930 ARRAY_MAPPER (exp, double, ::exp);
931 ARRAY_MAPPER (expm1, double, octave::math::expm1);
932 ARRAY_MAPPER (fix, double, octave::math::fix);
933 ARRAY_MAPPER (floor, double, ::floor);
934 ARRAY_MAPPER (log, Complex, octave::math::rc_log);
935 ARRAY_MAPPER (log2, Complex, octave::math::rc_log2);
936 ARRAY_MAPPER (log10, Complex, octave::math::rc_log10);
937 ARRAY_MAPPER (log1p, Complex, octave::math::rc_log1p);
938 ARRAY_MAPPER (round, double, octave::math::round);
939 ARRAY_MAPPER (roundb, double, octave::math::roundb);
940 ARRAY_MAPPER (signum, double, octave::math::signum);
941 ARRAY_MAPPER (sin, double, ::sin);
942 ARRAY_MAPPER (sinh, double, ::sinh);
943 ARRAY_MAPPER (sqrt, Complex, octave::math::rc_sqrt);
944 ARRAY_MAPPER (tan, double, ::tan);
945 ARRAY_MAPPER (tanh, double, ::tanh);
946 ARRAY_MAPPER (isnan, bool, octave::math::isnan);
947 ARRAY_MAPPER (isna, bool, octave::math::isna);
948 ARRAY_MAPPER (isinf, bool, octave::math::isinf);
949 ARRAY_MAPPER (isfinite, bool, octave::math::isfinite);
950
951 default: // Attempt to go via dense matrix.
952 return octave_base_sparse<SparseMatrix>::map (umap);
953 }
954 }
955