5164
|
1 // Template sparse classes |
|
2 /* |
|
3 |
|
4 Copyright (C) 2004 David Bateman |
|
5 Copyright (C) 1998-2004 Andy Adler |
|
6 |
|
7 Octave is free software; you can redistribute it and/or modify it |
|
8 under the terms of the GNU General Public License as published by the |
|
9 Free Software Foundation; either version 2, or (at your option) any |
|
10 later version. |
|
11 |
|
12 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
15 for more details. |
|
16 |
|
17 You should have received a copy of the GNU General Public License |
|
18 along with this program; see the file COPYING. If not, write to the Free |
|
19 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
|
20 |
|
21 */ |
|
22 |
|
23 #if !defined (octave_Sparse_h) |
|
24 #define octave_Sparse_h 1 |
|
25 |
|
26 #include <cassert> |
|
27 #include <cstddef> |
|
28 |
|
29 #include <iostream> |
|
30 |
|
31 #include "Array.h" |
|
32 #include "Array2.h" |
|
33 #include "dim-vector.h" |
|
34 #include "lo-utils.h" |
|
35 |
|
36 class idx_vector; |
|
37 |
|
38 // Two dimensional sparse class. Handles the reference counting for |
|
39 // all the derived classes. |
|
40 |
|
41 template <class T> |
|
42 class |
|
43 Sparse |
|
44 { |
|
45 protected: |
|
46 //-------------------------------------------------------------------- |
|
47 // The real representation of all Sparse arrays. |
|
48 //-------------------------------------------------------------------- |
|
49 |
|
50 class SparseRep |
|
51 { |
|
52 public: |
|
53 |
|
54 T *d; |
5275
|
55 octave_idx_type *r; |
|
56 octave_idx_type *c; |
|
57 octave_idx_type nnz; |
|
58 octave_idx_type nrows; |
|
59 octave_idx_type ncols; |
5164
|
60 int count; |
|
61 |
5275
|
62 SparseRep (void) : d (0), r (0), c (new octave_idx_type [1]), nnz (0), nrows (0), |
5164
|
63 ncols (0), count (1) { c[0] = 0; } |
|
64 |
5275
|
65 SparseRep (octave_idx_type n) : d (0), r (0), c (new octave_idx_type [n+1]), nnz (0), nrows (n), |
5164
|
66 ncols (n), count (1) |
|
67 { |
5275
|
68 for (octave_idx_type i = 0; i < n + 1; i++) |
5164
|
69 c[i] = 0; |
|
70 } |
|
71 |
5275
|
72 SparseRep (octave_idx_type nr, octave_idx_type nc) : d (0), r (0), c (new octave_idx_type [nc+1]), nnz (0), |
5164
|
73 nrows (nr), ncols (nc), count (1) |
|
74 { |
5275
|
75 for (octave_idx_type i = 0; i < nc + 1; i++) |
5164
|
76 c[i] = 0; |
|
77 } |
|
78 |
5275
|
79 SparseRep (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz) : d (new T [nz]), |
|
80 r (new octave_idx_type [nz]), c (new octave_idx_type [nc+1]), nnz (nz), nrows (nr), |
5164
|
81 ncols (nc), count (1) |
|
82 { |
5275
|
83 for (octave_idx_type i = 0; i < nc + 1; i++) |
5164
|
84 c[i] = 0; |
|
85 } |
|
86 |
|
87 SparseRep (const SparseRep& a) |
5275
|
88 : d (new T [a.nnz]), r (new octave_idx_type [a.nnz]), c (new octave_idx_type [a.ncols + 1]), |
5164
|
89 nnz (a.nnz), nrows (a.nrows), ncols (a.ncols), count (1) |
|
90 { |
5275
|
91 for (octave_idx_type i = 0; i < nnz; i++) |
5164
|
92 { |
|
93 d[i] = a.d[i]; |
|
94 r[i] = a.r[i]; |
|
95 } |
5275
|
96 for (octave_idx_type i = 0; i < ncols + 1; i++) |
5164
|
97 c[i] = a.c[i]; |
|
98 } |
|
99 |
|
100 ~SparseRep (void) { delete [] d; delete [] r; delete [] c; } |
|
101 |
5275
|
102 octave_idx_type length (void) const { return nnz; } |
5164
|
103 |
5275
|
104 octave_idx_type nonzero (void) const { return c [ncols]; } |
5164
|
105 |
5275
|
106 T& elem (octave_idx_type _r, octave_idx_type _c); |
5164
|
107 |
5275
|
108 T celem (octave_idx_type _r, octave_idx_type _c) const; |
5164
|
109 |
5275
|
110 T& data (octave_idx_type i) { return d[i]; } |
5164
|
111 |
5275
|
112 T cdata (octave_idx_type i) const { return d[i]; } |
5164
|
113 |
5275
|
114 octave_idx_type& ridx (octave_idx_type i) { return r[i]; } |
5164
|
115 |
5275
|
116 octave_idx_type cridx (octave_idx_type i) const { return r[i]; } |
5164
|
117 |
5275
|
118 octave_idx_type& cidx (octave_idx_type i) { return c[i]; } |
5164
|
119 |
5275
|
120 octave_idx_type ccidx (octave_idx_type i) const { return c[i]; } |
5164
|
121 |
|
122 void maybe_compress (bool remove_zeros); |
|
123 |
5275
|
124 void change_length (octave_idx_type nz); |
5164
|
125 |
|
126 private: |
|
127 |
|
128 // No assignment! |
|
129 |
|
130 SparseRep& operator = (const SparseRep& a); |
|
131 }; |
|
132 |
|
133 //-------------------------------------------------------------------- |
|
134 |
|
135 void make_unique (void) |
|
136 { |
|
137 if (rep->count > 1) |
|
138 { |
|
139 --rep->count; |
|
140 rep = new SparseRep (*rep); |
|
141 } |
|
142 } |
|
143 |
|
144 public: |
|
145 |
|
146 // !!! WARNING !!! -- these should be protected, not public. You |
|
147 // should not access these data members directly! |
|
148 |
|
149 typename Sparse<T>::SparseRep *rep; |
|
150 |
|
151 dim_vector dimensions; |
|
152 |
|
153 protected: |
|
154 idx_vector *idx; |
5275
|
155 octave_idx_type idx_count; |
5164
|
156 |
|
157 private: |
|
158 |
|
159 typename Sparse<T>::SparseRep *nil_rep (void) const |
|
160 { |
|
161 static typename Sparse<T>::SparseRep *nr |
|
162 = new typename Sparse<T>::SparseRep (); |
|
163 |
|
164 nr->count++; |
|
165 |
|
166 return nr; |
|
167 } |
|
168 |
|
169 public: |
|
170 |
|
171 Sparse (void) |
|
172 : rep (nil_rep ()), dimensions (dim_vector(0,0)), |
|
173 idx (0), idx_count (0) { } |
|
174 |
5275
|
175 explicit Sparse (octave_idx_type n) |
5164
|
176 : rep (new typename Sparse<T>::SparseRep (n)), |
|
177 dimensions (dim_vector (n, n)), idx (0), idx_count (0) { } |
|
178 |
5275
|
179 explicit Sparse (octave_idx_type nr, octave_idx_type nc) |
5164
|
180 : rep (new typename Sparse<T>::SparseRep (nr, nc)), |
|
181 dimensions (dim_vector (nr, nc)), idx (0), idx_count (0) { } |
|
182 |
5275
|
183 explicit Sparse (octave_idx_type nr, octave_idx_type nc, T val); |
5164
|
184 |
5275
|
185 Sparse (const dim_vector& dv, octave_idx_type nz) |
5164
|
186 : rep (new typename Sparse<T>::SparseRep (dv(0), dv(1), nz)), |
|
187 dimensions (dv), idx (0), idx_count (0) { } |
|
188 |
5275
|
189 Sparse (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz) |
5164
|
190 : rep (new typename Sparse<T>::SparseRep (nr, nc, nz)), |
|
191 dimensions (dim_vector (nr, nc)), idx (0), idx_count (0) { } |
|
192 |
|
193 // Type conversion case. |
|
194 template <class U> Sparse (const Sparse<U>& a); |
|
195 |
|
196 // No type conversion case. |
|
197 Sparse (const Sparse<T>& a) |
|
198 : rep (a.rep), dimensions (a.dimensions), idx (0), idx_count (0) |
|
199 { |
|
200 rep->count++; |
|
201 } |
|
202 |
|
203 public: |
|
204 |
|
205 Sparse (const dim_vector& dv); |
|
206 |
|
207 Sparse (const Sparse<T>& a, const dim_vector& dv); |
|
208 |
5275
|
209 Sparse (const Array<T>& a, const Array<octave_idx_type>& r, const Array<octave_idx_type>& c, |
|
210 octave_idx_type nr, octave_idx_type nc, bool sum_terms); |
5164
|
211 |
|
212 Sparse (const Array<T>& a, const Array<double>& r, const Array<double>& c, |
5275
|
213 octave_idx_type nr, octave_idx_type nc, bool sum_terms); |
5164
|
214 |
|
215 // Sparsify a normal matrix |
|
216 Sparse (const Array2<T>& a); |
|
217 Sparse (const Array<T>& a); |
|
218 |
|
219 virtual ~Sparse (void); |
|
220 |
|
221 Sparse<T>& operator = (const Sparse<T>& a) |
|
222 { |
|
223 if (this != &a) |
|
224 { |
|
225 if (--rep->count <= 0) |
|
226 delete rep; |
|
227 |
|
228 rep = a.rep; |
|
229 rep->count++; |
|
230 |
|
231 dimensions = a.dimensions; |
|
232 } |
|
233 |
|
234 idx_count = 0; |
|
235 idx = 0; |
|
236 |
|
237 return *this; |
|
238 } |
|
239 |
|
240 // Note that capacity and nnz are the amount of storage for non-zero |
|
241 // elements, while nonzero is the actual number of non-zero terms |
5275
|
242 octave_idx_type capacity (void) const { return rep->length (); } |
|
243 octave_idx_type nnz (void) const { return capacity (); } |
|
244 octave_idx_type nonzero (void) const { return rep->nonzero (); } |
5164
|
245 |
|
246 // Paranoid number of elements test for case of dims = (-1,-1) |
5275
|
247 octave_idx_type numel (void) const |
5164
|
248 { |
|
249 if (dim1() < 0 || dim2() < 0) |
|
250 return 0; |
|
251 else |
|
252 return dimensions.numel (); |
|
253 } |
|
254 |
5275
|
255 octave_idx_type nelem (void) const { return capacity (); } |
|
256 octave_idx_type length (void) const { return numel (); } |
5164
|
257 |
5275
|
258 octave_idx_type dim1 (void) const { return dimensions(0); } |
|
259 octave_idx_type dim2 (void) const { return dimensions(1); } |
5164
|
260 |
5275
|
261 octave_idx_type rows (void) const { return dim1 (); } |
|
262 octave_idx_type cols (void) const { return dim2 (); } |
|
263 octave_idx_type columns (void) const { return dim2 (); } |
5164
|
264 |
5275
|
265 octave_idx_type get_row_index (octave_idx_type k) { return ridx (k); } |
|
266 octave_idx_type get_col_index (octave_idx_type k) |
5164
|
267 { |
5275
|
268 octave_idx_type ret = 0; |
5164
|
269 while (cidx(ret+1) < k) |
|
270 ret++; |
|
271 return ret; |
|
272 } |
5275
|
273 size_t byte_size (void) const { return (cols () + 1) * sizeof (octave_idx_type) + |
|
274 capacity () * (sizeof (T) + sizeof (octave_idx_type)); } |
5164
|
275 |
|
276 dim_vector dims (void) const { return dimensions; } |
|
277 |
|
278 Sparse<T> squeeze (void) const { return *this; } |
|
279 |
5275
|
280 octave_idx_type compute_index (const Array<octave_idx_type>& ra_idx) const; |
5164
|
281 |
5275
|
282 T range_error (const char *fcn, octave_idx_type n) const; |
|
283 T& range_error (const char *fcn, octave_idx_type n); |
5164
|
284 |
5275
|
285 T range_error (const char *fcn, octave_idx_type i, octave_idx_type j) const; |
|
286 T& range_error (const char *fcn, octave_idx_type i, octave_idx_type j); |
5164
|
287 |
5275
|
288 T range_error (const char *fcn, const Array<octave_idx_type>& ra_idx) const; |
|
289 T& range_error (const char *fcn, const Array<octave_idx_type>& ra_idx); |
5164
|
290 |
|
291 // No checking, even for multiple references, ever. |
|
292 |
5275
|
293 T& xelem (octave_idx_type n) |
5164
|
294 { |
5275
|
295 octave_idx_type i = n % rows (), j = n / rows(); |
5164
|
296 return xelem (i, j); |
|
297 } |
|
298 |
5275
|
299 T xelem (octave_idx_type n) const |
5164
|
300 { |
5275
|
301 octave_idx_type i = n % rows (), j = n / rows(); |
5164
|
302 return xelem (i, j); |
|
303 } |
|
304 |
5275
|
305 T& xelem (octave_idx_type i, octave_idx_type j) { return rep->elem (i, j); } |
|
306 T xelem (octave_idx_type i, octave_idx_type j) const { return rep->celem (i, j); } |
5164
|
307 |
5275
|
308 T& xelem (const Array<octave_idx_type>& ra_idx) |
5164
|
309 { return xelem (compute_index (ra_idx)); } |
|
310 |
5275
|
311 T xelem (const Array<octave_idx_type>& ra_idx) const |
5164
|
312 { return xelem (compute_index (ra_idx)); } |
|
313 |
|
314 // XXX FIXME XXX -- would be nice to fix this so that we don't |
|
315 // unnecessarily force a copy, but that is not so easy, and I see no |
|
316 // clean way to do it. |
|
317 |
5275
|
318 T& checkelem (octave_idx_type n) |
5164
|
319 { |
|
320 if (n < 0 || n >= numel ()) |
|
321 return range_error ("T& Sparse<T>::checkelem", n); |
|
322 else |
|
323 { |
|
324 make_unique (); |
|
325 return xelem (n); |
|
326 } |
|
327 } |
|
328 |
5275
|
329 T& checkelem (octave_idx_type i, octave_idx_type j) |
5164
|
330 { |
|
331 if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ()) |
|
332 return range_error ("T& Sparse<T>::checkelem", i, j); |
|
333 else |
|
334 { |
|
335 make_unique (); |
|
336 return xelem (i, j); |
|
337 } |
|
338 } |
|
339 |
5275
|
340 T& checkelem (const Array<octave_idx_type>& ra_idx) |
5164
|
341 { |
5275
|
342 octave_idx_type i = compute_index (ra_idx); |
5164
|
343 |
|
344 if (i < 0) |
|
345 return range_error ("T& Sparse<T>::checkelem", ra_idx); |
|
346 else |
|
347 return elem (i); |
|
348 } |
|
349 |
5275
|
350 T& elem (octave_idx_type n) |
5164
|
351 { |
|
352 make_unique (); |
|
353 return xelem (n); |
|
354 } |
|
355 |
5275
|
356 T& elem (octave_idx_type i, octave_idx_type j) |
5164
|
357 { |
|
358 make_unique (); |
|
359 return xelem (i, j); |
|
360 } |
|
361 |
5275
|
362 T& elem (const Array<octave_idx_type>& ra_idx) |
5164
|
363 { return Sparse<T>::elem (compute_index (ra_idx)); } |
|
364 |
|
365 #if defined (BOUNDS_CHECKING) |
5275
|
366 T& operator () (octave_idx_type n) { return checkelem (n); } |
|
367 T& operator () (octave_idx_type i, octave_idx_type j) { return checkelem (i, j); } |
|
368 T& operator () (const Array<octave_idx_type>& ra_idx) { return checkelem (ra_idx); } |
5164
|
369 #else |
5275
|
370 T& operator () (octave_idx_type n) { return elem (n); } |
|
371 T& operator () (octave_idx_type i, octave_idx_type j) { return elem (i, j); } |
|
372 T& operator () (const Array<octave_idx_type>& ra_idx) { return elem (ra_idx); } |
5164
|
373 #endif |
|
374 |
5275
|
375 T checkelem (octave_idx_type n) const |
5164
|
376 { |
|
377 if (n < 0 || n >= numel ()) |
|
378 return range_error ("T Sparse<T>::checkelem", n); |
|
379 else |
|
380 return xelem (n); |
|
381 } |
|
382 |
5275
|
383 T checkelem (octave_idx_type i, octave_idx_type j) const |
5164
|
384 { |
|
385 if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ()) |
|
386 return range_error ("T Sparse<T>::checkelem", i, j); |
|
387 else |
|
388 return xelem (i, j); |
|
389 } |
|
390 |
5275
|
391 T checkelem (const Array<octave_idx_type>& ra_idx) const |
5164
|
392 { |
5275
|
393 octave_idx_type i = compute_index (ra_idx); |
5164
|
394 |
|
395 if (i < 0) |
|
396 return range_error ("T Sparse<T>::checkelem", ra_idx); |
|
397 else |
|
398 return Sparse<T>::elem (i); |
|
399 } |
|
400 |
5275
|
401 T elem (octave_idx_type n) const { return xelem (n); } |
5164
|
402 |
5275
|
403 T elem (octave_idx_type i, octave_idx_type j) const { return xelem (i, j); } |
5164
|
404 |
5275
|
405 T elem (const Array<octave_idx_type>& ra_idx) const |
5164
|
406 { return Sparse<T>::elem (compute_index (ra_idx)); } |
|
407 |
|
408 #if defined (BOUNDS_CHECKING) |
5275
|
409 T operator () (octave_idx_type n) const { return checkelem (n); } |
|
410 T operator () (octave_idx_type i, octave_idx_type j) const { return checkelem (i, j); } |
|
411 T operator () (const Array<octave_idx_type>& ra_idx) const { return checkelem (ra_idx); } |
5164
|
412 #else |
5275
|
413 T operator () (octave_idx_type n) const { return elem (n); } |
|
414 T operator () (octave_idx_type i, octave_idx_type j) const { return elem (i, j); } |
|
415 T operator () (const Array<octave_idx_type>& ra_idx) const { return elem (ra_idx); } |
5164
|
416 #endif |
|
417 |
|
418 Sparse<T> maybe_compress (bool remove_zeros = false) |
|
419 { rep->maybe_compress (remove_zeros); return (*this); } |
|
420 |
|
421 Sparse<T> reshape (const dim_vector& new_dims) const; |
|
422 |
|
423 // !!! WARNING !!! -- the following resize_no_fill functions are |
|
424 // public because template friends don't work properly with versions |
|
425 // of gcc earlier than 3.3. You should use these functions only in |
|
426 // classes that are derived from Sparse<T>. |
|
427 |
|
428 // protected: |
|
429 |
5275
|
430 void resize_no_fill (octave_idx_type r, octave_idx_type c); |
5164
|
431 |
|
432 void resize_no_fill (const dim_vector& dv); |
|
433 |
|
434 public: |
5275
|
435 Sparse<T> permute (const Array<octave_idx_type>& vec, bool inv = false) const; |
5164
|
436 |
5275
|
437 Sparse<T> ipermute (const Array<octave_idx_type>& vec) const |
5164
|
438 { return permute (vec, true); } |
|
439 |
5275
|
440 void resize (octave_idx_type r, octave_idx_type c) { resize_no_fill (r, c); } |
5164
|
441 |
|
442 void resize (const dim_vector& dv) { resize_no_fill (dv); } |
|
443 |
5275
|
444 void change_capacity (octave_idx_type nz) { rep->change_length (nz); } |
5164
|
445 |
5275
|
446 Sparse<T>& insert (const Sparse<T>& a, octave_idx_type r, octave_idx_type c); |
|
447 Sparse<T>& insert (const Sparse<T>& a, const Array<octave_idx_type>& idx); |
5164
|
448 |
|
449 bool is_square (void) const { return (dim1 () == dim2 ()); } |
|
450 |
|
451 bool is_empty (void) const { return (rows () < 1 && cols () < 1); } |
|
452 |
|
453 Sparse<T> transpose (void) const; |
|
454 |
|
455 T* data (void) { make_unique (); return rep->d; } |
5275
|
456 T& data (octave_idx_type i) { make_unique (); return rep->data (i); } |
5164
|
457 T* xdata (void) { return rep->d; } |
5275
|
458 T& xdata (octave_idx_type i) { return rep->data (i); } |
5164
|
459 |
5275
|
460 T data (octave_idx_type i) const { return rep->data (i); } |
5164
|
461 T* data (void) const { return rep->d; } |
|
462 |
5275
|
463 octave_idx_type* ridx (void) { make_unique (); return rep->r; } |
|
464 octave_idx_type& ridx (octave_idx_type i) { make_unique (); return rep->ridx (i); } |
|
465 octave_idx_type* xridx (void) { return rep->r; } |
|
466 octave_idx_type& xridx (octave_idx_type i) { return rep->ridx (i); } |
5164
|
467 |
5275
|
468 octave_idx_type ridx (octave_idx_type i) const { return rep->cridx (i); } |
|
469 octave_idx_type* ridx (void) const { return rep->r; } |
5164
|
470 |
5275
|
471 octave_idx_type* cidx (void) { make_unique (); return rep->c; } |
|
472 octave_idx_type& cidx (octave_idx_type i) { make_unique (); return rep->cidx (i); } |
|
473 octave_idx_type* xcidx (void) { return rep->c; } |
|
474 octave_idx_type& xcidx (octave_idx_type i) { return rep->cidx (i); } |
5164
|
475 |
5275
|
476 octave_idx_type cidx (octave_idx_type i) const { return rep->ccidx (i); } |
|
477 octave_idx_type* cidx (void) const { return rep->c; } |
5164
|
478 |
5275
|
479 octave_idx_type ndims (void) const { return dimensions.length (); } |
5164
|
480 |
|
481 void clear_index (void); |
|
482 |
|
483 void set_index (const idx_vector& i); |
|
484 |
5275
|
485 octave_idx_type index_count (void) const { return idx_count; } |
5164
|
486 |
|
487 idx_vector *get_idx (void) const { return idx; } |
|
488 |
|
489 void maybe_delete_elements (idx_vector& i); |
|
490 |
|
491 void maybe_delete_elements (idx_vector& i, idx_vector& j); |
|
492 |
|
493 void maybe_delete_elements (Array<idx_vector>& ra_idx); |
|
494 |
|
495 Sparse<T> value (void); |
|
496 |
|
497 Sparse<T> index (idx_vector& i, int resize_ok = 0) const; |
|
498 |
|
499 Sparse<T> index (idx_vector& i, idx_vector& j, int resize_ok = 0) const; |
|
500 |
|
501 Sparse<T> index (Array<idx_vector>& ra_idx, int resize_ok = 0) const; |
|
502 |
|
503 void print_info (std::ostream& os, const std::string& prefix) const; |
|
504 |
|
505 }; |
|
506 |
|
507 // NOTE: these functions should be friends of the Sparse<T> class and |
|
508 // Sparse<T>::dimensions should be protected, not public, but we can't |
|
509 // do that because of bugs in gcc prior to 3.3. |
|
510 |
|
511 template <class LT, class RT> |
|
512 /* friend */ int |
|
513 assign (Sparse<LT>& lhs, const Sparse<RT>& rhs); |
|
514 |
|
515 template <class LT, class RT> |
|
516 /* friend */ int |
|
517 assign1 (Sparse<LT>& lhs, const Sparse<RT>& rhs); |
|
518 |
|
519 #define INSTANTIATE_SPARSE_ASSIGN(LT, RT) \ |
|
520 template int assign (Sparse<LT>&, const Sparse<RT>&); \ |
|
521 template int assign1 (Sparse<LT>&, const Sparse<RT>&); |
|
522 |
|
523 #define INSTANTIATE_SPARSE(T) \ |
|
524 template class Sparse<T>; |
|
525 |
|
526 #define INSTANTIATE_SPARSE_AND_ASSIGN(T) \ |
|
527 INSTANTIATE_SPARSE (T); \ |
|
528 INSTANTIATE_SPARSE_ASSIGN (T, T) |
|
529 |
|
530 #endif |
|
531 |
|
532 /* |
|
533 ;;; Local Variables: *** |
|
534 ;;; mode: C++ *** |
|
535 ;;; End: *** |
|
536 */ |