00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #if !defined (octave_Sparse_h)
00027 #define octave_Sparse_h 1
00028
00029 #include <cassert>
00030 #include <cstddef>
00031
00032 #include <iosfwd>
00033 #include <algorithm>
00034
00035 #include "Array.h"
00036 #include "dim-vector.h"
00037 #include "lo-error.h"
00038 #include "lo-utils.h"
00039
00040 #include "oct-sort.h"
00041 #include "oct-mem.h"
00042
00043 class idx_vector;
00044
00045
00046
00047
00048 template <class T>
00049 class
00050 Sparse
00051 {
00052 public:
00053
00054 typedef T element_type;
00055
00056 protected:
00057
00058
00059
00060
00061 class OCTAVE_API SparseRep
00062 {
00063 public:
00064
00065 T *d;
00066 octave_idx_type *r;
00067 octave_idx_type *c;
00068 octave_idx_type nzmx;
00069 octave_idx_type nrows;
00070 octave_idx_type ncols;
00071 octave_refcount<int> count;
00072
00073 SparseRep (void) : d (0), r (0), c (new octave_idx_type [1]), nzmx (0), nrows (0),
00074 ncols (0), count (1) { c[0] = 0; }
00075
00076 SparseRep (octave_idx_type n) : d (0), r (0), c (new octave_idx_type [n+1]), nzmx (0), nrows (n),
00077 ncols (n), count (1)
00078 {
00079 for (octave_idx_type i = 0; i < n + 1; i++)
00080 c[i] = 0;
00081 }
00082
00083 SparseRep (octave_idx_type nr, octave_idx_type nc) : d (0), r (0), c (new octave_idx_type [nc+1]), nzmx (0),
00084 nrows (nr), ncols (nc), count (1)
00085 {
00086 for (octave_idx_type i = 0; i < nc + 1; i++)
00087 c[i] = 0;
00088 }
00089
00090 SparseRep (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz) : d (new T [nz]),
00091 r (new octave_idx_type [nz]), c (new octave_idx_type [nc+1]), nzmx (nz), nrows (nr),
00092 ncols (nc), count (1)
00093 {
00094 for (octave_idx_type i = 0; i < nc + 1; i++)
00095 c[i] = 0;
00096 }
00097
00098 SparseRep (const SparseRep& a)
00099 : d (new T [a.nzmx]), r (new octave_idx_type [a.nzmx]), c (new octave_idx_type [a.ncols + 1]),
00100 nzmx (a.nzmx), nrows (a.nrows), ncols (a.ncols), count (1)
00101 {
00102 octave_idx_type nz = a.nnz ();
00103 copy_or_memcpy (nz, a.d, d);
00104 copy_or_memcpy (nz, a.r, r);
00105 copy_or_memcpy (ncols + 1, a.c, c);
00106 }
00107
00108 ~SparseRep (void) { delete [] d; delete [] r; delete [] c; }
00109
00110 octave_idx_type length (void) const { return nzmx; }
00111
00112 octave_idx_type nnz (void) const { return c [ncols]; }
00113
00114 T& elem (octave_idx_type _r, octave_idx_type _c);
00115
00116 T celem (octave_idx_type _r, octave_idx_type _c) const;
00117
00118 T& data (octave_idx_type i) { return d[i]; }
00119
00120 T cdata (octave_idx_type i) const { return d[i]; }
00121
00122 octave_idx_type& ridx (octave_idx_type i) { return r[i]; }
00123
00124 octave_idx_type cridx (octave_idx_type i) const { return r[i]; }
00125
00126 octave_idx_type& cidx (octave_idx_type i) { return c[i]; }
00127
00128 octave_idx_type ccidx (octave_idx_type i) const { return c[i]; }
00129
00130 void maybe_compress (bool remove_zeros);
00131
00132 void change_length (octave_idx_type nz);
00133
00134 bool indices_ok (void) const;
00135
00136 private:
00137
00138
00139
00140 SparseRep& operator = (const SparseRep& a);
00141 };
00142
00143
00144
00145 void make_unique (void)
00146 {
00147 if (rep->count > 1)
00148 {
00149 --rep->count;
00150 rep = new SparseRep (*rep);
00151 }
00152 }
00153
00154 public:
00155
00156
00157
00158
00159 typename Sparse<T>::SparseRep *rep;
00160
00161 dim_vector dimensions;
00162
00163 private:
00164
00165 typename Sparse<T>::SparseRep *nil_rep (void) const
00166 {
00167 static typename Sparse<T>::SparseRep *nr
00168 = new typename Sparse<T>::SparseRep ();
00169
00170 nr->count++;
00171
00172 return nr;
00173 }
00174
00175 public:
00176
00177 Sparse (void)
00178 : rep (nil_rep ()), dimensions (dim_vector(0,0)) { }
00179
00180 explicit Sparse (octave_idx_type n)
00181 : rep (new typename Sparse<T>::SparseRep (n)),
00182 dimensions (dim_vector (n, n)) { }
00183
00184 explicit Sparse (octave_idx_type nr, octave_idx_type nc)
00185 : rep (new typename Sparse<T>::SparseRep (nr, nc)),
00186 dimensions (dim_vector (nr, nc)) { }
00187
00188 explicit Sparse (octave_idx_type nr, octave_idx_type nc, T val);
00189
00190 Sparse (const dim_vector& dv, octave_idx_type nz)
00191 : rep (new typename Sparse<T>::SparseRep (dv(0), dv(1), nz)),
00192 dimensions (dv) { }
00193
00194 Sparse (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz)
00195 : rep (new typename Sparse<T>::SparseRep (nr, nc, nz)),
00196 dimensions (dim_vector (nr, nc)) { }
00197
00198
00199 template <class U>
00200 Sparse (const Sparse<U>& a)
00201 : rep (new typename Sparse<T>::SparseRep (a.rep->nrows, a.rep->ncols, a.rep->nzmx)),
00202 dimensions (a.dimensions)
00203 {
00204 octave_idx_type nz = a.nnz ();
00205 std::copy (a.rep->d, a.rep->d + nz, rep->d);
00206 copy_or_memcpy (nz, a.rep->r, rep->r);
00207 copy_or_memcpy (rep->ncols + 1, a.rep->c, rep->c);
00208 }
00209
00210
00211 Sparse (const Sparse<T>& a)
00212 : rep (a.rep), dimensions (a.dimensions)
00213 {
00214 rep->count++;
00215 }
00216
00217 public:
00218
00219 Sparse (const dim_vector& dv);
00220
00221 Sparse (const Sparse<T>& a, const dim_vector& dv);
00222
00223 Sparse (const Array<T>& a, const idx_vector& r, const idx_vector& c,
00224 octave_idx_type nr = -1, octave_idx_type nc = -1,
00225 bool sum_terms = true, octave_idx_type nzm = -1);
00226
00227
00228 Sparse (const Array<T>& a);
00229
00230 virtual ~Sparse (void);
00231
00232 Sparse<T>& operator = (const Sparse<T>& a);
00233
00234
00235
00236
00237 octave_idx_type nzmax (void) const { return rep->length (); }
00238 octave_idx_type capacity (void) const { return nzmax (); }
00239 octave_idx_type nnz (void) const { return rep->nnz (); }
00240
00241
00242
00243 octave_idx_type numel (void) const
00244 {
00245 return dimensions.safe_numel ();
00246 }
00247
00248 octave_idx_type nelem (void) const { return capacity (); }
00249 octave_idx_type length (void) const { return numel (); }
00250
00251 octave_idx_type dim1 (void) const { return dimensions(0); }
00252 octave_idx_type dim2 (void) const { return dimensions(1); }
00253
00254 octave_idx_type rows (void) const { return dim1 (); }
00255 octave_idx_type cols (void) const { return dim2 (); }
00256 octave_idx_type columns (void) const { return dim2 (); }
00257
00258 octave_idx_type get_row_index (octave_idx_type k) { return ridx (k); }
00259 octave_idx_type get_col_index (octave_idx_type k)
00260 {
00261 octave_idx_type ret = 0;
00262 while (cidx(ret+1) < k)
00263 ret++;
00264 return ret;
00265 }
00266
00267 size_t byte_size (void) const
00268 {
00269 return (static_cast<size_t>(cols () + 1) * sizeof (octave_idx_type)
00270 + static_cast<size_t> (capacity ()) * (sizeof (T) + sizeof (octave_idx_type)));
00271 }
00272
00273 dim_vector dims (void) const { return dimensions; }
00274
00275 Sparse<T> squeeze (void) const { return *this; }
00276
00277 octave_idx_type compute_index (const Array<octave_idx_type>& ra_idx) const;
00278
00279 T range_error (const char *fcn, octave_idx_type n) const;
00280 T& range_error (const char *fcn, octave_idx_type n);
00281
00282 T range_error (const char *fcn, octave_idx_type i, octave_idx_type j) const;
00283 T& range_error (const char *fcn, octave_idx_type i, octave_idx_type j);
00284
00285 T range_error (const char *fcn, const Array<octave_idx_type>& ra_idx) const;
00286 T& range_error (const char *fcn, const Array<octave_idx_type>& ra_idx);
00287
00288
00289
00290 T& xelem (octave_idx_type n)
00291 {
00292 octave_idx_type i = n % rows (), j = n / rows();
00293 return xelem (i, j);
00294 }
00295
00296 T xelem (octave_idx_type n) const
00297 {
00298 octave_idx_type i = n % rows (), j = n / rows();
00299 return xelem (i, j);
00300 }
00301
00302 T& xelem (octave_idx_type i, octave_idx_type j) { return rep->elem (i, j); }
00303 T xelem (octave_idx_type i, octave_idx_type j) const { return rep->celem (i, j); }
00304
00305 T& xelem (const Array<octave_idx_type>& ra_idx)
00306 { return xelem (compute_index (ra_idx)); }
00307
00308 T xelem (const Array<octave_idx_type>& ra_idx) const
00309 { return xelem (compute_index (ra_idx)); }
00310
00311
00312
00313
00314
00315 T& checkelem (octave_idx_type n)
00316 {
00317 if (n < 0 || n >= numel ())
00318 return range_error ("T& Sparse<T>::checkelem", n);
00319 else
00320 {
00321 make_unique ();
00322 return xelem (n);
00323 }
00324 }
00325
00326 T& checkelem (octave_idx_type i, octave_idx_type j)
00327 {
00328 if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ())
00329 return range_error ("T& Sparse<T>::checkelem", i, j);
00330 else
00331 {
00332 make_unique ();
00333 return xelem (i, j);
00334 }
00335 }
00336
00337 T& checkelem (const Array<octave_idx_type>& ra_idx)
00338 {
00339 octave_idx_type i = compute_index (ra_idx);
00340
00341 if (i < 0)
00342 return range_error ("T& Sparse<T>::checkelem", ra_idx);
00343 else
00344 return elem (i);
00345 }
00346
00347 T& elem (octave_idx_type n)
00348 {
00349 make_unique ();
00350 return xelem (n);
00351 }
00352
00353 T& elem (octave_idx_type i, octave_idx_type j)
00354 {
00355 make_unique ();
00356 return xelem (i, j);
00357 }
00358
00359 T& elem (const Array<octave_idx_type>& ra_idx)
00360 { return Sparse<T>::elem (compute_index (ra_idx)); }
00361
00362 #if defined (BOUNDS_CHECKING)
00363 T& operator () (octave_idx_type n) { return checkelem (n); }
00364 T& operator () (octave_idx_type i, octave_idx_type j) { return checkelem (i, j); }
00365 T& operator () (const Array<octave_idx_type>& ra_idx) { return checkelem (ra_idx); }
00366 #else
00367 T& operator () (octave_idx_type n) { return elem (n); }
00368 T& operator () (octave_idx_type i, octave_idx_type j) { return elem (i, j); }
00369 T& operator () (const Array<octave_idx_type>& ra_idx) { return elem (ra_idx); }
00370 #endif
00371
00372 T checkelem (octave_idx_type n) const
00373 {
00374 if (n < 0 || n >= numel ())
00375 return range_error ("T Sparse<T>::checkelem", n);
00376 else
00377 return xelem (n);
00378 }
00379
00380 T checkelem (octave_idx_type i, octave_idx_type j) const
00381 {
00382 if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ())
00383 return range_error ("T Sparse<T>::checkelem", i, j);
00384 else
00385 return xelem (i, j);
00386 }
00387
00388 T checkelem (const Array<octave_idx_type>& ra_idx) const
00389 {
00390 octave_idx_type i = compute_index (ra_idx);
00391
00392 if (i < 0)
00393 return range_error ("T Sparse<T>::checkelem", ra_idx);
00394 else
00395 return Sparse<T>::elem (i);
00396 }
00397
00398 T elem (octave_idx_type n) const { return xelem (n); }
00399
00400 T elem (octave_idx_type i, octave_idx_type j) const { return xelem (i, j); }
00401
00402 T elem (const Array<octave_idx_type>& ra_idx) const
00403 { return Sparse<T>::elem (compute_index (ra_idx)); }
00404
00405 #if defined (BOUNDS_CHECKING)
00406 T operator () (octave_idx_type n) const { return checkelem (n); }
00407 T operator () (octave_idx_type i, octave_idx_type j) const { return checkelem (i, j); }
00408 T operator () (const Array<octave_idx_type>& ra_idx) const { return checkelem (ra_idx); }
00409 #else
00410 T operator () (octave_idx_type n) const { return elem (n); }
00411 T operator () (octave_idx_type i, octave_idx_type j) const { return elem (i, j); }
00412 T operator () (const Array<octave_idx_type>& ra_idx) const { return elem (ra_idx); }
00413 #endif
00414
00415 Sparse<T> maybe_compress (bool remove_zeros = false)
00416 {
00417 if (remove_zeros)
00418 make_unique ();
00419
00420 rep->maybe_compress (remove_zeros);
00421 return (*this);
00422 }
00423
00424 Sparse<T> reshape (const dim_vector& new_dims) const;
00425
00426 Sparse<T> permute (const Array<octave_idx_type>& vec, bool inv = false) const;
00427
00428 Sparse<T> ipermute (const Array<octave_idx_type>& vec) const
00429 { return permute (vec, true); }
00430
00431 void resize1 (octave_idx_type n);
00432
00433 void resize (octave_idx_type r, octave_idx_type c);
00434
00435 void resize (const dim_vector& dv);
00436
00437 void change_capacity (octave_idx_type nz)
00438 {
00439 if (nz < nnz ())
00440 make_unique ();
00441 rep->change_length (nz);
00442 }
00443
00444 Sparse<T>& insert (const Sparse<T>& a, octave_idx_type r, octave_idx_type c);
00445 Sparse<T>& insert (const Sparse<T>& a, const Array<octave_idx_type>& idx);
00446
00447 bool is_square (void) const { return (dim1 () == dim2 ()); }
00448
00449 bool is_empty (void) const { return (rows () < 1 && cols () < 1); }
00450
00451 Sparse<T> transpose (void) const;
00452
00453 T* data (void) { make_unique (); return rep->d; }
00454 T& data (octave_idx_type i) { make_unique (); return rep->data (i); }
00455 T* xdata (void) { return rep->d; }
00456 T& xdata (octave_idx_type i) { return rep->data (i); }
00457
00458 T data (octave_idx_type i) const { return rep->data (i); }
00459
00460 T* data (void) const { return rep->d; }
00461
00462 octave_idx_type* ridx (void) { make_unique (); return rep->r; }
00463 octave_idx_type& ridx (octave_idx_type i) { make_unique (); return rep->ridx (i); }
00464 octave_idx_type* xridx (void) { return rep->r; }
00465 octave_idx_type& xridx (octave_idx_type i) { return rep->ridx (i); }
00466
00467 octave_idx_type ridx (octave_idx_type i) const { return rep->cridx (i); }
00468
00469 octave_idx_type* ridx (void) const { return rep->r; }
00470
00471 octave_idx_type* cidx (void) { make_unique (); return rep->c; }
00472 octave_idx_type& cidx (octave_idx_type i) { make_unique (); return rep->cidx (i); }
00473 octave_idx_type* xcidx (void) { return rep->c; }
00474 octave_idx_type& xcidx (octave_idx_type i) { return rep->cidx (i); }
00475
00476 octave_idx_type cidx (octave_idx_type i) const { return rep->ccidx (i); }
00477
00478 octave_idx_type* cidx (void) const { return rep->c; }
00479
00480 octave_idx_type ndims (void) const { return dimensions.length (); }
00481
00482 void delete_elements (const idx_vector& i);
00483
00484 void delete_elements (int dim, const idx_vector& i);
00485
00486 void delete_elements (const idx_vector& i, const idx_vector& j);
00487
00488 Sparse<T> index (const idx_vector& i, bool resize_ok = false) const;
00489
00490 Sparse<T> index (const idx_vector& i, const idx_vector& j, bool resize_ok = false) const;
00491
00492 void assign (const idx_vector& i, const Sparse<T>& rhs);
00493
00494 void assign (const idx_vector& i, const idx_vector& j, const Sparse<T>& rhs);
00495
00496 void print_info (std::ostream& os, const std::string& prefix) const;
00497
00498
00499
00500 void *mex_get_data (void) const { return const_cast<T *> (data ()); }
00501
00502 octave_idx_type *mex_get_ir (void) const { return const_cast<octave_idx_type *> (ridx ()); }
00503
00504 octave_idx_type *mex_get_jc (void) const { return const_cast<octave_idx_type *> (cidx ()); }
00505
00506 Sparse<T> sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const;
00507 Sparse<T> sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
00508 sortmode mode = ASCENDING) const;
00509
00510 Sparse<T> diag (octave_idx_type k = 0) const;
00511
00512
00513 static Sparse<T>
00514 cat (int dim, octave_idx_type n, const Sparse<T> *sparse_list);
00515
00516 Array<T> array_value (void) const;
00517
00518 template <class U, class F>
00519 Sparse<U>
00520 map (F fcn) const
00521 {
00522 Sparse<U> result;
00523 U f_zero = fcn (0.);
00524
00525 if (f_zero != 0.)
00526 {
00527 octave_idx_type nr = rows ();
00528 octave_idx_type nc = cols ();
00529
00530 result = Sparse<U> (nr, nc, f_zero);
00531
00532 for (octave_idx_type j = 0; j < nc; j++)
00533 for (octave_idx_type i = cidx(j); i < cidx (j+1); i++)
00534 {
00535 octave_quit ();
00536
00537 result.data (ridx (i) + j * nr) = fcn (data(i));
00538 }
00539
00540 result.maybe_compress (true);
00541 }
00542 else
00543 {
00544 octave_idx_type nz = nnz ();
00545 octave_idx_type nr = rows ();
00546 octave_idx_type nc = cols ();
00547
00548 result = Sparse<U> (nr, nc, nz);
00549 octave_idx_type ii = 0;
00550 result.cidx (ii) = 0;
00551
00552 for (octave_idx_type j = 0; j < nc; j++)
00553 {
00554 for (octave_idx_type i = cidx(j); i < cidx (j+1); i++)
00555 {
00556 U val = fcn (data (i));
00557 if (val != 0.0)
00558 {
00559 result.data (ii) = val;
00560 result.ridx (ii++) = ridx (i);
00561 }
00562 octave_quit ();
00563 }
00564 result.cidx (j+1) = ii;
00565 }
00566
00567 result.maybe_compress (false);
00568 }
00569
00570 return result;
00571 }
00572
00573
00574 template <class U>
00575 Sparse<U>
00576 map (U (&fcn) (T)) const
00577 { return map<U, U (&) (T)> (fcn); }
00578
00579 template <class U>
00580 Sparse<U>
00581 map (U (&fcn) (const T&)) const
00582 { return map<U, U (&) (const T&)> (fcn); }
00583
00584 bool indices_ok (void) const { return rep->indices_ok (); }
00585 };
00586
00587 template<typename T>
00588 std::istream&
00589 read_sparse_matrix (std::istream& is, Sparse<T>& a,
00590 T (*read_fcn) (std::istream&))
00591 {
00592 octave_idx_type nr = a.rows ();
00593 octave_idx_type nc = a.cols ();
00594 octave_idx_type nz = a.nzmax ();
00595
00596 if (nr > 0 && nc > 0)
00597 {
00598 octave_idx_type itmp;
00599 octave_idx_type jtmp;
00600 octave_idx_type iold = 0;
00601 octave_idx_type jold = 0;
00602 octave_idx_type ii = 0;
00603 T tmp;
00604
00605 a.cidx (0) = 0;
00606 for (octave_idx_type i = 0; i < nz; i++)
00607 {
00608 itmp = 0; jtmp = 0;
00609 is >> itmp;
00610 itmp--;
00611
00612 is >> jtmp;
00613 jtmp--;
00614
00615 if (itmp < 0 || itmp >= nr)
00616 {
00617 (*current_liboctave_error_handler)
00618 ("invalid sparse matrix: row index = %d out of range",
00619 itmp + 1);
00620 is.setstate (std::ios::failbit);
00621 goto done;
00622 }
00623
00624 if (jtmp < 0 || jtmp >= nc)
00625 {
00626 (*current_liboctave_error_handler)
00627 ("invalid sparse matrix: column index = %d out of range",
00628 jtmp + 1);
00629 is.setstate (std::ios::failbit);
00630 goto done;
00631 }
00632
00633 if (jtmp < jold)
00634 {
00635 (*current_liboctave_error_handler)
00636 ("invalid sparse matrix: column indices must appear in ascending order");
00637 is.setstate (std::ios::failbit);
00638 goto done;
00639 }
00640 else if (jtmp > jold)
00641 {
00642 for (octave_idx_type j = jold; j < jtmp; j++)
00643 a.cidx(j+1) = ii;
00644 }
00645 else if (itmp < iold)
00646 {
00647 (*current_liboctave_error_handler)
00648 ("invalid sparse matrix: row indices must appear in ascending order in each column");
00649 is.setstate (std::ios::failbit);
00650 goto done;
00651 }
00652
00653 iold = itmp;
00654 jold = jtmp;
00655
00656 tmp = read_fcn (is);
00657
00658 if (is)
00659 {
00660 a.data (ii) = tmp;
00661 a.ridx (ii++) = itmp;
00662 }
00663 else
00664 goto done;
00665 }
00666
00667 for (octave_idx_type j = jold; j < nc; j++)
00668 a.cidx(j+1) = ii;
00669 }
00670
00671 done:
00672
00673 return is;
00674 }
00675
00676 #endif