Mercurial > octave-nkf

comparison liboctave/Sparse.h @ 10312:cbc402e64d83

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
untabify liboctave header files
author John W. Eaton <jwe@octave.org>
date Thu, 11 Feb 2010 12:14:48 -0500
parents 4c0cdbe0acca
children a3635bc1ea19
comparison
equal deleted inserted replaced
10311:a217e1d74353 10312:cbc402e64d83
67 octave_idx_type nrows; 67 octave_idx_type nrows;
68 octave_idx_type ncols; 68 octave_idx_type ncols;
69 int count; 69 int count;
70 70
71 SparseRep (void) : d (0), r (0), c (new octave_idx_type [1]), nzmx (0), nrows (0), 71 SparseRep (void) : d (0), r (0), c (new octave_idx_type [1]), nzmx (0), nrows (0),
72 ncols (0), count (1) { c[0] = 0; } 72 ncols (0), count (1) { c[0] = 0; }
73 73
74 SparseRep (octave_idx_type n) : d (0), r (0), c (new octave_idx_type [n+1]), nzmx (0), nrows (n), 74 SparseRep (octave_idx_type n) : d (0), r (0), c (new octave_idx_type [n+1]), nzmx (0), nrows (n),
75 ncols (n), count (1) 75 ncols (n), count (1)
76 { 76 {
77 for (octave_idx_type i = 0; i < n + 1; i++) 77 for (octave_idx_type i = 0; i < n + 1; i++)
78 c[i] = 0; 78 c[i] = 0;
79 } 79 }
80 80
81 SparseRep (octave_idx_type nr, octave_idx_type nc) : d (0), r (0), c (new octave_idx_type [nc+1]), nzmx (0), 81 SparseRep (octave_idx_type nr, octave_idx_type nc) : d (0), r (0), c (new octave_idx_type [nc+1]), nzmx (0),
82 nrows (nr), ncols (nc), count (1) 82 nrows (nr), ncols (nc), count (1)
83 { 83 {
84 for (octave_idx_type i = 0; i < nc + 1; i++) 84 for (octave_idx_type i = 0; i < nc + 1; i++)
85 c[i] = 0; 85 c[i] = 0;
86 } 86 }
87 87
88 SparseRep (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz) : d (new T [nz]), 88 SparseRep (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz) : d (new T [nz]),
89 r (new octave_idx_type [nz]), c (new octave_idx_type [nc+1]), nzmx (nz), nrows (nr), 89 r (new octave_idx_type [nz]), c (new octave_idx_type [nc+1]), nzmx (nz), nrows (nr),
90 ncols (nc), count (1) 90 ncols (nc), count (1)
91 { 91 {
92 for (octave_idx_type i = 0; i < nc + 1; i++) 92 for (octave_idx_type i = 0; i < nc + 1; i++)
93 c[i] = 0; 93 c[i] = 0;
94 } 94 }
95 95
96 SparseRep (const SparseRep& a) 96 SparseRep (const SparseRep& a)
97 : d (new T [a.nzmx]), r (new octave_idx_type [a.nzmx]), c (new octave_idx_type [a.ncols + 1]), 97 : d (new T [a.nzmx]), r (new octave_idx_type [a.nzmx]), c (new octave_idx_type [a.ncols + 1]),
98 nzmx (a.nzmx), nrows (a.nrows), ncols (a.ncols), count (1) 98 nzmx (a.nzmx), nrows (a.nrows), ncols (a.ncols), count (1)
99 { 99 {
100 for (octave_idx_type i = 0; i < nzmx; i++) 100 for (octave_idx_type i = 0; i < nzmx; i++)
101 { 101 {
102 d[i] = a.d[i]; 102 d[i] = a.d[i];
103 r[i] = a.r[i]; 103 r[i] = a.r[i];
104 } 104 }
105 for (octave_idx_type i = 0; i < ncols + 1; i++) 105 for (octave_idx_type i = 0; i < ncols + 1; i++)
106 c[i] = a.c[i]; 106 c[i] = a.c[i];
107 } 107 }
108 108
109 ~SparseRep (void) { delete [] d; delete [] r; delete [] c; } 109 ~SparseRep (void) { delete [] d; delete [] r; delete [] c; }
110 110
111 octave_idx_type length (void) const { return nzmx; } 111 octave_idx_type length (void) const { return nzmx; }
144 //-------------------------------------------------------------------- 144 //--------------------------------------------------------------------
145 145
146 void make_unique (void) 146 void make_unique (void)
147 { 147 {
148 if (rep->count > 1) 148 if (rep->count > 1)
149 { 149 {
150 --rep->count; 150 --rep->count;
151 rep = new SparseRep (*rep); 151 rep = new SparseRep (*rep);
152 } 152 }
153 } 153 }
154 154
155 public: 155 public:
156 156
157 // !!! WARNING !!! -- these should be protected, not public. You 157 // !!! WARNING !!! -- these should be protected, not public. You
168 private: 168 private:
169 169
170 typename Sparse<T>::SparseRep *nil_rep (void) const 170 typename Sparse<T>::SparseRep *nil_rep (void) const
171 { 171 {
172 static typename Sparse<T>::SparseRep *nr 172 static typename Sparse<T>::SparseRep *nr
173 = new typename Sparse<T>::SparseRep (); 173 = new typename Sparse<T>::SparseRep ();
174 174
175 nr->count++; 175 nr->count++;
176 176
177 return nr; 177 return nr;
178 } 178 }
216 Sparse (const dim_vector& dv); 216 Sparse (const dim_vector& dv);
217 217
218 Sparse (const Sparse<T>& a, const dim_vector& dv); 218 Sparse (const Sparse<T>& a, const dim_vector& dv);
219 219
220 Sparse (const Array<T>& a, const Array<octave_idx_type>& r, const Array<octave_idx_type>& c, 220 Sparse (const Array<T>& a, const Array<octave_idx_type>& r, const Array<octave_idx_type>& c,
221 octave_idx_type nr, octave_idx_type nc, bool sum_terms); 221 octave_idx_type nr, octave_idx_type nc, bool sum_terms);
222 222
223 Sparse (const Array<T>& a, const Array<double>& r, const Array<double>& c, 223 Sparse (const Array<T>& a, const Array<double>& r, const Array<double>& c,
224 octave_idx_type nr, octave_idx_type nc, bool sum_terms); 224 octave_idx_type nr, octave_idx_type nc, bool sum_terms);
225 225
226 // Sparsify a normal matrix 226 // Sparsify a normal matrix
227 Sparse (const Array2<T>& a); 227 Sparse (const Array2<T>& a);
228 Sparse (const Array<T>& a); 228 Sparse (const Array<T>& a);
229 229
311 // clean way to do it. 311 // clean way to do it.
312 312
313 T& checkelem (octave_idx_type n) 313 T& checkelem (octave_idx_type n)
314 { 314 {
315 if (n < 0 || n >= numel ()) 315 if (n < 0 || n >= numel ())
316 return range_error ("T& Sparse<T>::checkelem", n); 316 return range_error ("T& Sparse<T>::checkelem", n);
317 else 317 else
318 { 318 {
319 make_unique (); 319 make_unique ();
320 return xelem (n); 320 return xelem (n);
321 } 321 }
322 } 322 }
323 323
324 T& checkelem (octave_idx_type i, octave_idx_type j) 324 T& checkelem (octave_idx_type i, octave_idx_type j)
325 { 325 {
326 if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ()) 326 if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ())
327 return range_error ("T& Sparse<T>::checkelem", i, j); 327 return range_error ("T& Sparse<T>::checkelem", i, j);
328 else 328 else
329 { 329 {
330 make_unique (); 330 make_unique ();
331 return xelem (i, j); 331 return xelem (i, j);
332 } 332 }
333 } 333 }
334 334
335 T& checkelem (const Array<octave_idx_type>& ra_idx) 335 T& checkelem (const Array<octave_idx_type>& ra_idx)
336 { 336 {
337 octave_idx_type i = compute_index (ra_idx); 337 octave_idx_type i = compute_index (ra_idx);
338 338
339 if (i < 0) 339 if (i < 0)
340 return range_error ("T& Sparse<T>::checkelem", ra_idx); 340 return range_error ("T& Sparse<T>::checkelem", ra_idx);
341 else 341 else
342 return elem (i); 342 return elem (i);
343 } 343 }
344 344
345 T& elem (octave_idx_type n) 345 T& elem (octave_idx_type n)
346 { 346 {
347 make_unique (); 347 make_unique ();
368 #endif 368 #endif
369 369
370 T checkelem (octave_idx_type n) const 370 T checkelem (octave_idx_type n) const
371 { 371 {
372 if (n < 0 || n >= numel ()) 372 if (n < 0 || n >= numel ())
373 return range_error ("T Sparse<T>::checkelem", n); 373 return range_error ("T Sparse<T>::checkelem", n);
374 else 374 else
375 return xelem (n); 375 return xelem (n);
376 } 376 }
377 377
378 T checkelem (octave_idx_type i, octave_idx_type j) const 378 T checkelem (octave_idx_type i, octave_idx_type j) const
379 { 379 {
380 if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ()) 380 if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ())
381 return range_error ("T Sparse<T>::checkelem", i, j); 381 return range_error ("T Sparse<T>::checkelem", i, j);
382 else 382 else
383 return xelem (i, j); 383 return xelem (i, j);
384 } 384 }
385 385
386 T checkelem (const Array<octave_idx_type>& ra_idx) const 386 T checkelem (const Array<octave_idx_type>& ra_idx) const
387 { 387 {
388 octave_idx_type i = compute_index (ra_idx); 388 octave_idx_type i = compute_index (ra_idx);
389 389
390 if (i < 0) 390 if (i < 0)
391 return range_error ("T Sparse<T>::checkelem", ra_idx); 391 return range_error ("T Sparse<T>::checkelem", ra_idx);
392 else 392 else
393 return Sparse<T>::elem (i); 393 return Sparse<T>::elem (i);
394 } 394 }
395 395
396 T elem (octave_idx_type n) const { return xelem (n); } 396 T elem (octave_idx_type n) const { return xelem (n); }
397 397
398 T elem (octave_idx_type i, octave_idx_type j) const { return xelem (i, j); } 398 T elem (octave_idx_type i, octave_idx_type j) const { return xelem (i, j); }
508 508
509 octave_idx_type *mex_get_jc (void) const { return const_cast<octave_idx_type *> (cidx ()); } 509 octave_idx_type *mex_get_jc (void) const { return const_cast<octave_idx_type *> (cidx ()); }
510 510
511 Sparse<T> sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const; 511 Sparse<T> sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const;
512 Sparse<T> sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0, 512 Sparse<T> sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
513 sortmode mode = ASCENDING) const; 513 sortmode mode = ASCENDING) const;
514 514
515 Sparse<T> diag (octave_idx_type k = 0) const; 515 Sparse<T> diag (octave_idx_type k = 0) const;
516 516
517 template <class U, class F> 517 template <class U, class F>
518 Sparse<U> 518 Sparse<U>
521 Sparse<U> result; 521 Sparse<U> result;
522 U f_zero = fcn (0.); 522 U f_zero = fcn (0.);
523 523
524 if (f_zero != 0.) 524 if (f_zero != 0.)
525 { 525 {
526 octave_idx_type nr = rows (); 526 octave_idx_type nr = rows ();
527 octave_idx_type nc = cols (); 527 octave_idx_type nc = cols ();
528 528
529 result = Sparse<U> (nr, nc, f_zero); 529 result = Sparse<U> (nr, nc, f_zero);
530 530
531 for (octave_idx_type j = 0; j < nc; j++) 531 for (octave_idx_type j = 0; j < nc; j++)
532 for (octave_idx_type i = cidx(j); i < cidx (j+1); i++) 532 for (octave_idx_type i = cidx(j); i < cidx (j+1); i++)
533 { 533 {
534 octave_quit (); 534 octave_quit ();
535 /* Use data instead of elem for better performance. */ 535 /* Use data instead of elem for better performance. */
536 result.data (ridx (i) + j * nr) = fcn (data(i)); 536 result.data (ridx (i) + j * nr) = fcn (data(i));
537 } 537 }
538 538
539 result.maybe_compress (true); 539 result.maybe_compress (true);
540 } 540 }
541 else 541 else
542 { 542 {
543 octave_idx_type nz = nnz (); 543 octave_idx_type nz = nnz ();
544 octave_idx_type nr = rows (); 544 octave_idx_type nr = rows ();
545 octave_idx_type nc = cols (); 545 octave_idx_type nc = cols ();
546 546
547 result = Sparse<U> (nr, nc, nz); 547 result = Sparse<U> (nr, nc, nz);
548 octave_idx_type ii = 0; 548 octave_idx_type ii = 0;
549 result.cidx (ii) = 0; 549 result.cidx (ii) = 0;
550 550
551 for (octave_idx_type j = 0; j < nc; j++) 551 for (octave_idx_type j = 0; j < nc; j++)
552 { 552 {
553 for (octave_idx_type i = cidx(j); i < cidx (j+1); i++) 553 for (octave_idx_type i = cidx(j); i < cidx (j+1); i++)
554 { 554 {
555 U val = fcn (data (i)); 555 U val = fcn (data (i));
556 if (val != 0.0) 556 if (val != 0.0)
557 { 557 {
558 result.data (ii) = val; 558 result.data (ii) = val;
559 result.ridx (ii++) = ridx (i); 559 result.ridx (ii++) = ridx (i);
560 } 560 }
561 octave_quit (); 561 octave_quit ();
562 } 562 }
563 result.cidx (j+1) = ii; 563 result.cidx (j+1) = ii;
564 } 564 }
565 565
566 result.maybe_compress (false); 566 result.maybe_compress (false);
567 } 567 }
568 568
569 return result; 569 return result;
570 } 570 }
571 571
596 assign1 (Sparse<LT>& lhs, const Sparse<RT>& rhs); 596 assign1 (Sparse<LT>& lhs, const Sparse<RT>& rhs);
597 597
598 template<typename T> 598 template<typename T>
599 std::istream& 599 std::istream&
600 read_sparse_matrix (std::istream& is, Sparse<T>& a, 600 read_sparse_matrix (std::istream& is, Sparse<T>& a,
601 T (*read_fcn) (std::istream&)) 601 T (*read_fcn) (std::istream&))
602 { 602 {
603 octave_idx_type nr = a.rows (); 603 octave_idx_type nr = a.rows ();
604 octave_idx_type nc = a.cols (); 604 octave_idx_type nc = a.cols ();
605 octave_idx_type nz = a.nzmax (); 605 octave_idx_type nz = a.nzmax ();
606 606
613 octave_idx_type ii = 0; 613 octave_idx_type ii = 0;
614 T tmp; 614 T tmp;
615 615
616 a.cidx (0) = 0; 616 a.cidx (0) = 0;
617 for (octave_idx_type i = 0; i < nz; i++) 617 for (octave_idx_type i = 0; i < nz; i++)
618 { 618 {
619 itmp = 0; jtmp = 0; 619 itmp = 0; jtmp = 0;
620 is >> itmp; 620 is >> itmp;
621 itmp--; 621 itmp--;
622 622
623 is >> jtmp; 623 is >> jtmp;
624 jtmp--; 624 jtmp--;
625 625
626 if (itmp < 0 || itmp >= nr) 626 if (itmp < 0 || itmp >= nr)
627 { 627 {
628 (*current_liboctave_error_handler) 628 (*current_liboctave_error_handler)
629 ("invalid sparse matrix: row index = %d out of range", 629 ("invalid sparse matrix: row index = %d out of range",
630 itmp + 1); 630 itmp + 1);
631 is.setstate (std::ios::failbit); 631 is.setstate (std::ios::failbit);
632 goto done; 632 goto done;
633 } 633 }
634 634
635 if (jtmp < 0 || jtmp >= nc) 635 if (jtmp < 0 || jtmp >= nc)
636 { 636 {
637 (*current_liboctave_error_handler) 637 (*current_liboctave_error_handler)
638 ("invalid sparse matrix: column index = %d out of range", 638 ("invalid sparse matrix: column index = %d out of range",
639 jtmp + 1); 639 jtmp + 1);
640 is.setstate (std::ios::failbit); 640 is.setstate (std::ios::failbit);
641 goto done; 641 goto done;
642 } 642 }
643 643
644 if (jtmp < jold) 644 if (jtmp < jold)
645 { 645 {
646 (*current_liboctave_error_handler) 646 (*current_liboctave_error_handler)
647 ("invalid sparse matrix: column indices must appear in ascending order"); 647 ("invalid sparse matrix: column indices must appear in ascending order");
648 is.setstate (std::ios::failbit); 648 is.setstate (std::ios::failbit);
649 goto done; 649 goto done;
650 } 650 }
651 else if (jtmp > jold) 651 else if (jtmp > jold)
652 { 652 {
653 for (octave_idx_type j = jold; j < jtmp; j++) 653 for (octave_idx_type j = jold; j < jtmp; j++)
654 a.cidx(j+1) = ii; 654 a.cidx(j+1) = ii;
655 } 655 }
656 else if (itmp < iold) 656 else if (itmp < iold)
657 { 657 {
658 (*current_liboctave_error_handler) 658 (*current_liboctave_error_handler)
659 ("invalid sparse matrix: row indices must appear in ascending order in each column"); 659 ("invalid sparse matrix: row indices must appear in ascending order in each column");
660 is.setstate (std::ios::failbit); 660 is.setstate (std::ios::failbit);
661 goto done; 661 goto done;
662 } 662 }
663 663
664 iold = itmp; 664 iold = itmp;
665 jold = jtmp; 665 jold = jtmp;
666 666
667 tmp = read_fcn (is); 667 tmp = read_fcn (is);
668 668
669 if (is) 669 if (is)
670 { 670 {
671 a.data (ii) = tmp; 671 a.data (ii) = tmp;
672 a.ridx (ii++) = itmp; 672 a.ridx (ii++) = itmp;
673 } 673 }
674 else 674 else
675 goto done; 675 goto done;
676 } 676 }
677 677
678 for (octave_idx_type j = jold; j < nc; j++) 678 for (octave_idx_type j = jold; j < nc; j++)
679 a.cidx(j+1) = ii; 679 a.cidx(j+1) = ii;
680 } 680 }
681 681
682 done: 682 done:
683 683
684 return is; 684 return is;