IT++ Logo
smat.h
Go to the documentation of this file.
1 
29 #ifndef SMAT_H
30 #define SMAT_H
31 
32 #include <itpp/base/svec.h>
33 #include <itpp/itexports.h>
34 
35 namespace itpp
36 {
37 
38 // Declaration of class Vec
39 template <class T> class Vec;
40 // Declaration of class Mat
41 template <class T> class Mat;
42 // Declaration of class Sparse_Vec
43 template <class T> class Sparse_Vec;
44 // Declaration of class Sparse_Mat
45 template <class T> class Sparse_Mat;
46 
47 // ------------------------ Sparse_Mat Friends -------------------------------------
48 
50 template <class T>
51 Sparse_Mat<T> operator+(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
52 
54 template <class T>
55 Sparse_Mat<T> operator*(const T &c, const Sparse_Mat<T> &m);
56 
58 template <class T>
59 Sparse_Mat<T> operator*(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
60 
62 template <class T>
63 Sparse_Vec<T> operator*(const Sparse_Mat<T> &m, const Sparse_Vec<T> &v);
64 
66 template <class T>
67 Vec<T> operator*(const Sparse_Mat<T> &m, const Vec<T> &v);
68 
70 template <class T>
71 Vec<T> operator*(const Vec<T> &v, const Sparse_Mat<T> &m);
72 
74 template <class T>
75 Mat<T> trans_mult(const Sparse_Mat<T> &m);
76 
78 template <class T>
79 Sparse_Mat<T> trans_mult_s(const Sparse_Mat<T> &m);
80 
82 template <class T>
83 Sparse_Mat<T> trans_mult(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
84 
86 template <class T>
87 Vec<T> trans_mult(const Sparse_Mat<T> &m, const Vec<T> &v);
88 
90 template <class T>
91 Sparse_Mat<T> mult_trans(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
92 
106 template <class T>
108 {
109 public:
110 
112  Sparse_Mat();
113 
124  Sparse_Mat(int rows, int cols, int row_data_init = 200);
125 
127  Sparse_Mat(const Sparse_Mat<T> &m);
128 
130  Sparse_Mat(const Mat<T> &m);
131 
137  Sparse_Mat(const Mat<T> &m, T epsilon);
138 
140  ~Sparse_Mat();
141 
152  void set_size(int rows, int cols, int row_data_init = -1);
153 
155  int rows() const { return n_rows; }
156 
158  int cols() const { return n_cols; }
159 
161  int nnz();
162 
164  double density();
165 
167  void compact();
168 
170  void full(Mat<T> &m) const;
171 
173  Mat<T> full() const;
174 
176  T operator()(int r, int c) const;
177 
179  void set(int r, int c, T v);
180 
182  void set_new(int r, int c, T v);
183 
185  void add_elem(const int r, const int c, const T v);
186 
188  void zeros();
189 
191  void zero_elem(const int r, const int c);
192 
194  void clear();
195 
197  void clear_elem(const int r, const int c);
198 
200  void set_submatrix(int r1, int r2, int c1, int c2, const Mat<T> &m);
201 
203  void set_submatrix(int r, int c, const Mat<T>& m);
204 
206  Sparse_Mat<T> get_submatrix(int r1, int r2, int c1, int c2) const;
207 
209  Sparse_Mat<T> get_submatrix_cols(int c1, int c2) const;
210 
212  void get_col(int c, Sparse_Vec<T> &v) const;
213 
215  Sparse_Vec<T> get_col(int c) const;
216 
218  void set_col(int c, const Sparse_Vec<T> &v);
219 
224  void transpose(Sparse_Mat<T> &m) const;
225 
230  Sparse_Mat<T> transpose() const;
231 
236  // Sparse_Mat<T> T() const { return this->transpose(); };
237 
239  void operator=(const Sparse_Mat<T> &m);
240 
242  void operator=(const Mat<T> &m);
243 
245  Sparse_Mat<T> operator-() const;
246 
248  bool operator==(const Sparse_Mat<T> &m) const;
249 
251  void operator+=(const Sparse_Mat<T> &v);
252 
254  void operator+=(const Mat<T> &v);
255 
257  void operator-=(const Sparse_Mat<T> &v);
258 
260  void operator-=(const Mat<T> &v);
261 
263  void operator*=(const T &v);
264 
266  void operator/=(const T &v);
267 
269  friend Sparse_Mat<T> operator+<>(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
270 
272  friend Sparse_Mat<T> operator*<>(const T &c, const Sparse_Mat<T> &m);
273 
275  friend Sparse_Mat<T> operator*<>(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
276 
278  friend Sparse_Vec<T> operator*<>(const Sparse_Mat<T> &m, const Sparse_Vec<T> &v);
279 
281  friend Vec<T> operator*<>(const Sparse_Mat<T> &m, const Vec<T> &v);
282 
284  friend Vec<T> operator*<>(const Vec<T> &v, const Sparse_Mat<T> &m);
285 
287  friend Mat<T> trans_mult <>(const Sparse_Mat<T> &m);
288 
290  friend Sparse_Mat<T> trans_mult_s <>(const Sparse_Mat<T> &m);
291 
293  friend Sparse_Mat<T> trans_mult <>(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
294 
296  friend Vec<T> trans_mult <>(const Sparse_Mat<T> &m, const Vec<T> &v);
297 
299  friend Sparse_Mat<T> mult_trans <>(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
300 
301 private:
302  void init();
303  void alloc_empty();
304  void alloc(int row_data_size = 200);
305  void free();
306 
307  int n_rows, n_cols;
308  Sparse_Vec<T> *col;
309 };
310 
315 typedef Sparse_Mat<int> sparse_imat;
316 
321 typedef Sparse_Mat<double> sparse_mat;
322 
327 typedef Sparse_Mat<std::complex<double> > sparse_cmat;
328 
329 //---------------------- Implementation starts here --------------------------------
330 
331 template <class T>
332 void Sparse_Mat<T>::init()
333 {
334  n_rows = 0;
335  n_cols = 0;
336  col = 0;
337 }
338 
339 template <class T>
340 void Sparse_Mat<T>::alloc_empty()
341 {
342  if (n_cols == 0)
343  col = 0;
344  else
345  col = new Sparse_Vec<T>[n_cols];
346 }
347 
348 template <class T>
349 void Sparse_Mat<T>::alloc(int row_data_init)
350 {
351  if (n_cols == 0)
352  col = 0;
353  else
354  col = new Sparse_Vec<T>[n_cols];
355  for (int c = 0; c < n_cols; c++)
356  col[c].set_size(n_rows, row_data_init);
357 }
358 
359 template <class T>
360 void Sparse_Mat<T>::free()
361 {
362  delete [] col;
363  col = 0;
364 }
365 
366 template <class T>
368 {
369  init();
370 }
371 
372 template <class T>
373 Sparse_Mat<T>::Sparse_Mat(int rows, int cols, int row_data_init)
374 {
375  init();
376  n_rows = rows;
377  n_cols = cols;
378  alloc(row_data_init);
379 }
380 
381 template <class T>
383 {
384  init();
385  n_rows = m.n_rows;
386  n_cols = m.n_cols;
387  alloc_empty();
388 
389  for (int c = 0; c < n_cols; c++)
390  col[c] = m.col[c];
391 }
392 
393 template <class T>
395 {
396  init();
397  n_rows = m.rows();
398  n_cols = m.cols();
399  alloc();
400 
401  for (int c = 0; c < n_cols; c++) {
402  for (int r = 0; r < n_rows; r++) {
403  //if (abs(m(r,c)) != T(0))
404  if (m(r, c) != T(0))
405  col[c].set_new(r, m(r, c));
406  }
407  col[c].compact();
408  }
409 }
410 
411 template <class T>
412 Sparse_Mat<T>::Sparse_Mat(const Mat<T> &m, T epsilon)
413 {
414  init();
415  n_rows = m.rows();
416  n_cols = m.cols();
417  alloc();
418 
419  for (int c = 0; c < n_cols; c++) {
420  for (int r = 0; r < n_rows; r++) {
421  if (std::abs(m(r, c)) > std::abs(epsilon))
422  col[c].set_new(r, m(r, c));
423  }
424  col[c].compact();
425  }
426 }
427 
428 template <class T>
430 {
431  free();
432 }
433 
434 template <class T>
435 void Sparse_Mat<T>::set_size(int rows, int cols, int row_data_init)
436 {
437  n_rows = rows;
438 
439  //Allocate new memory for data if the number of columns has changed or if row_data_init != -1
440  if (cols != n_cols || row_data_init != -1) {
441  n_cols = cols;
442  free();
443  alloc(row_data_init);
444  }
445 }
446 
447 template <class T>
449 {
450  int n = 0;
451  for (int c = 0; c < n_cols; c++)
452  n += col[c].nnz();
453 
454  return n;
455 }
456 
457 template <class T>
459 {
460  //return static_cast<double>(nnz())/(n_rows*n_cols);
461  return double(nnz()) / (n_rows*n_cols);
462 }
463 
464 template <class T>
466 {
467  for (int c = 0; c < n_cols; c++)
468  col[c].compact();
469 }
470 
471 template <class T>
473 {
474  m.set_size(n_rows, n_cols);
475  m = T(0);
476  for (int c = 0; c < n_cols; c++) {
477  for (int p = 0; p < col[c].nnz(); p++)
478  m(col[c].get_nz_index(p), c) = col[c].get_nz_data(p);
479  }
480 }
481 
482 template <class T>
484 {
485  Mat<T> r(n_rows, n_cols);
486  full(r);
487  return r;
488 }
489 
490 template <class T>
491 T Sparse_Mat<T>::operator()(int r, int c) const
492 {
493  it_assert_debug(r >= 0 && r<n_rows && c >= 0 && c < n_cols, "Incorrect input indexes given");
494  return col[c](r);
495 }
496 
497 template <class T>
498 void Sparse_Mat<T>::set(int r, int c, T v)
499 {
500  it_assert_debug(r >= 0 && r<n_rows && c >= 0 && c < n_cols, "Incorrect input indexes given");
501  col[c].set(r, v);
502 }
503 
504 template <class T>
505 void Sparse_Mat<T>::set_new(int r, int c, T v)
506 {
507  it_assert_debug(r >= 0 && r<n_rows && c >= 0 && c < n_cols, "Incorrect input indexes given");
508  col[c].set_new(r, v);
509 }
510 
511 template <class T>
512 void Sparse_Mat<T>::add_elem(int r, int c, T v)
513 {
514  it_assert_debug(r >= 0 && r<n_rows && c >= 0 && c < n_cols, "Incorrect input indexes given");
515  col[c].add_elem(r, v);
516 }
517 
518 template <class T>
520 {
521  for (int c = 0; c < n_cols; c++)
522  col[c].zeros();
523 }
524 
525 template <class T>
526 void Sparse_Mat<T>::zero_elem(const int r, const int c)
527 {
528  it_assert_debug(r >= 0 && r<n_rows && c >= 0 && c < n_cols, "Incorrect input indexes given");
529  col[c].zero_elem(r);
530 }
531 
532 template <class T>
534 {
535  for (int c = 0; c < n_cols; c++)
536  col[c].clear();
537 }
538 
539 template <class T>
540 void Sparse_Mat<T>::clear_elem(const int r, const int c)
541 {
542  it_assert_debug(r >= 0 && r<n_rows && c >= 0 && c < n_cols, "Incorrect input indexes given");
543  col[c].clear_elem(r);
544 }
545 
546 template <class T>
547 void Sparse_Mat<T>::set_submatrix(int r1, int r2, int c1, int c2, const Mat<T>& m)
548 {
549  if (r1 == -1) r1 = n_rows - 1;
550  if (r2 == -1) r2 = n_rows - 1;
551  if (c1 == -1) c1 = n_cols - 1;
552  if (c2 == -1) c2 = n_cols - 1;
553 
554  it_assert_debug(r1 >= 0 && r2 >= 0 && r1 < n_rows && r2 < n_rows &&
555  c1 >= 0 && c2 >= 0 && c1 < n_cols && c2 < n_cols, "Sparse_Mat<Num_T>::set_submatrix(): index out of range");
556 
557  it_assert_debug(r2 >= r1 && c2 >= c1, "Sparse_Mat<Num_T>::set_submatrix: r2<r1 or c2<c1");
558  it_assert_debug(m.rows() == r2 - r1 + 1 && m.cols() == c2 - c1 + 1, "Mat<Num_T>::set_submatrix(): sizes don't match");
559 
560  for (int i = 0 ; i < m.rows() ; i++) {
561  for (int j = 0 ; j < m.cols() ; j++) {
562  set(r1 + i, c1 + j, m(i, j));
563  }
564  }
565 }
566 
567 template <class T>
568 void Sparse_Mat<T>::set_submatrix(int r, int c, const Mat<T>& m)
569 {
570  it_assert_debug(r >= 0 && r + m.rows() <= n_rows &&
571  c >= 0 && c + m.cols() <= n_cols, "Sparse_Mat<Num_T>::set_submatrix(): index out of range");
572 
573  for (int i = 0 ; i < m.rows() ; i++) {
574  for (int j = 0 ; j < m.cols() ; j++) {
575  set(r + i, c + j, m(i, j));
576  }
577  }
578 }
579 
580 template <class T>
581 Sparse_Mat<T> Sparse_Mat<T>::get_submatrix(int r1, int r2, int c1, int c2) const
582 {
583  it_assert_debug(r1 <= r2 && r1 >= 0 && r1 < n_rows && c1 <= c2 && c1 >= 0 && c1 < n_cols,
584  "Sparse_Mat<T>::get_submatrix(): illegal input variables");
585 
586  Sparse_Mat<T> r(r2 - r1 + 1, c2 - c1 + 1);
587 
588  for (int c = c1; c <= c2; c++)
589  r.col[c-c1] = col[c].get_subvector(r1, r2);
590  r.compact();
591 
592  return r;
593 }
594 
595 template <class T>
597 {
598  it_assert_debug(c1 <= c2 && c1 >= 0 && c1 < n_cols, "Sparse_Mat<T>::get_submatrix_cols()");
599  Sparse_Mat<T> r(n_rows, c2 - c1 + 1, 0);
600 
601  for (int c = c1; c <= c2; c++)
602  r.col[c-c1] = col[c];
603  r.compact();
604 
605  return r;
606 }
607 
608 template <class T>
610 {
611  it_assert(c >= 0 && c < n_cols, "Sparse_Mat<T>::get_col()");
612  v = col[c];
613 }
614 
615 template <class T>
617 {
618  it_assert(c >= 0 && c < n_cols, "Sparse_Mat<T>::get_col()");
619  return col[c];
620 }
621 
622 template <class T>
624 {
625  it_assert(c >= 0 && c < n_cols, "Sparse_Mat<T>::set_col()");
626  col[c] = v;
627 }
628 
629 template <class T>
631 {
632  m.set_size(n_cols, n_rows);
633  for (int c = 0; c < n_cols; c++) {
634  for (int p = 0; p < col[c].nnz(); p++)
635  m.col[col[c].get_nz_index(p)].set_new(c, col[c].get_nz_data(p));
636  }
637 }
638 
639 template <class T>
641 {
642  Sparse_Mat<T> m;
643  transpose(m);
644  return m;
645 }
646 
647 template <class T>
649 {
650  free();
651  n_rows = m.n_rows;
652  n_cols = m.n_cols;
653  alloc_empty();
654 
655  for (int c = 0; c < n_cols; c++)
656  col[c] = m.col[c];
657 }
658 
659 template <class T>
661 {
662  free();
663  n_rows = m.rows();
664  n_cols = m.cols();
665  alloc();
666 
667  for (int c = 0; c < n_cols; c++) {
668  for (int r = 0; r < n_rows; r++) {
669  if (m(r, c) != T(0))
670  col[c].set_new(r, m(r, c));
671  }
672  col[c].compact();
673  }
674 }
675 
676 template <class T>
678 {
679  Sparse_Mat r(n_rows, n_cols, 0);
680 
681  for (int c = 0; c < n_cols; c++) {
682  r.col[c].resize_data(col[c].nnz());
683  for (int p = 0; p < col[c].nnz(); p++)
684  r.col[c].set_new(col[c].get_nz_index(p), -col[c].get_nz_data(p));
685  }
686 
687  return r;
688 }
689 
690 template <class T>
692 {
693  if (n_rows != m.n_rows || n_cols != m.n_cols)
694  return false;
695  for (int c = 0; c < n_cols; c++) {
696  if (!(col[c] == m.col[c]))
697  return false;
698  }
699  // If they passed all tests, they must be equal
700  return true;
701 }
702 
703 template <class T>
705 {
706  it_assert_debug(m.rows() == n_rows && m.cols() == n_cols, "Addition of unequal sized matrices is not allowed");
707 
708  Sparse_Vec<T> v;
709  for (int c = 0; c < n_cols; c++) {
710  m.get_col(c, v);
711  col[c] += v;
712  }
713 }
714 
715 template <class T>
717 {
718  it_assert_debug(m.rows() == n_rows && m.cols() == n_cols, "Addition of unequal sized matrices is not allowed");
719 
720  for (int c = 0; c < n_cols; c++)
721  col[c] += (m.get_col(c));
722 }
723 
724 template <class T>
726 {
727  it_assert_debug(m.rows() == n_rows && m.cols() == n_cols, "Subtraction of unequal sized matrices is not allowed");
728 
729  Sparse_Vec<T> v;
730  for (int c = 0; c < n_cols; c++) {
731  m.get_col(c, v);
732  col[c] -= v;
733  }
734 }
735 
736 template <class T>
738 {
739  it_assert_debug(m.rows() == n_rows && m.cols() == n_cols, "Subtraction of unequal sized matrices is not allowed");
740 
741  for (int c = 0; c < n_cols; c++)
742  col[c] -= (m.get_col(c));
743 }
744 
745 template <class T>
747 {
748  for (int c = 0; c < n_cols; c++)
749  col[c] *= m;
750 }
751 
752 template <class T>
754 {
755  for (int c = 0; c < n_cols; c++)
756  col[c] /= m;
757 }
758 
759 template <class T>
761 {
762  it_assert_debug(m1.n_cols == m2.n_cols && m1.n_rows == m2.n_rows , "Sparse_Mat<T> + Sparse_Mat<T>");
763 
764  Sparse_Mat<T> m(m1.n_rows, m1.n_cols, 0);
765 
766  for (int c = 0; c < m.n_cols; c++)
767  m.col[c] = m1.col[c] + m2.col[c];
768 
769  return m;
770 }
771 
772 // This function added by EGL, May'05
773 template <class T>
774 Sparse_Mat<T> operator*(const T &c, const Sparse_Mat<T> &m)
775 {
776  int i, j;
777  Sparse_Mat<T> ret(m.n_rows, m.n_cols);
778  for (j = 0; j < m.n_cols; j++) {
779  for (i = 0; i < m.col[j].nnz(); i++) {
780  T x = c * m.col[j].get_nz_data(i);
781  int k = m.col[j].get_nz_index(i);
782  ret.set_new(k, j, x);
783  }
784  }
785  return ret;
786 }
787 
788 template <class T>
790 {
791  it_assert_debug(m1.n_cols == m2.n_rows, "Sparse_Mat<T> * Sparse_Mat<T>");
792 
793  Sparse_Mat<T> ret(m1.n_rows, m2.n_cols);
794 
795  for (int c = 0; c < m2.n_cols; c++) {
796  Sparse_Vec<T> &m2colc = m2.col[c];
797  for (int p2 = 0; p2 < m2colc.nnz(); p2++) {
798  Sparse_Vec<T> &mcol = m1.col[m2colc.get_nz_index(p2)];
799  T x = m2colc.get_nz_data(p2);
800  for (int p1 = 0; p1 < mcol.nnz(); p1++) {
801  int r = mcol.get_nz_index(p1);
802  T inc = mcol.get_nz_data(p1) * x;
803  ret.col[c].add_elem(r, inc);
804  }
805  }
806  }
807  // old code
808  /* for (int c=0; c<m2.n_cols; c++) { */
809  /* for (int p2=0; p2<m2.col[c].nnz(); p2++) { */
810  /* Sparse_Vec<T> &mcol = m1.col[m2.col[c].get_nz_index(p2)]; */
811  /* for (int p1=0; p1<mcol.nnz(); p1++) { */
812  /* int r = mcol.get_nz_index(p1); */
813  /* T inc = mcol.get_nz_data(p1) * m2.col[c].get_nz_data(p2); */
814  /* ret.col[c].add_elem(r,inc); */
815  /* } */
816  /* } */
817  /* } */
818  ret.compact();
819  return ret;
820 }
821 
822 
823 // This is apparently buggy.
824 /* template <class T> */
825 /* Sparse_Mat<T> operator*(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2) */
826 /* { */
827 /* it_assert_debug(m1.n_cols == m2.n_rows, "Sparse_Mat<T> * Sparse_Mat<T>"); */
828 
829 /* Sparse_Mat<T> ret(m1.n_rows, m2.n_cols); */
830 /* ivec occupied_by(ret.n_rows), pos(ret.n_rows); */
831 /* for (int rp=0; rp<m1.n_rows; rp++) */
832 /* occupied_by[rp] = -1; */
833 /* for (int c=0; c<ret.n_cols; c++) { */
834 /* Sparse_Vec<T> &m2col = m2.col[c]; */
835 /* for (int p2=0; p2<m2col.nnz(); p2++) { */
836 /* Sparse_Vec<T> &m1col = m1.col[m2col.get_nz_index(p2)]; */
837 /* for (int p1=0; p1<m1col.nnz(); p1++) { */
838 /* int r = m1col.get_nz_index(p1); */
839 /* T inc = m1col.get_nz_data(p1) * m2col.get_nz_data(p2); */
840 /* if (occupied_by[r] == c) { */
841 /* int index=ret.col[c].get_nz_index(pos[r]); */
842 /* ret.col[c].add_elem(index,inc); */
843 /* } */
844 /* else { */
845 /* occupied_by[r] = c; */
846 /* pos[r] = ret.col[c].nnz(); */
847 /* ret.col[c].set_new(r, inc); */
848 /* } */
849 /* } */
850 /* } */
851 /* } */
852 /* ret.compact(); */
853 
854 /* return ret; */
855 /* } */
856 
857 
858 // This function added by EGL, May'05
859 template <class T>
861 {
862  it_assert_debug(m.n_cols == v.size(), "Sparse_Mat<T> * Sparse_Vec<T>");
863 
864  Sparse_Vec<T> ret(m.n_rows);
865 
866  /* The two lines below added because the input parameter "v" is
867  declared const, but the some functions (e.g., nnz()) change
868  the vector... Is there a better workaround? */
869  Sparse_Vec<T> vv = v;
870 
871  for (int p2 = 0; p2 < vv.nnz(); p2++) {
872  Sparse_Vec<T> &mcol = m.col[vv.get_nz_index(p2)];
873  T x = vv.get_nz_data(p2);
874  for (int p1 = 0; p1 < mcol.nnz(); p1++) {
875  int r = mcol.get_nz_index(p1);
876  T inc = mcol.get_nz_data(p1) * x;
877  ret.add_elem(r, inc);
878  }
879  }
880  ret.compact();
881  return ret;
882 }
883 
884 
885 template <class T>
887 {
888  it_assert_debug(m.n_cols == v.size(), "Sparse_Mat<T> * Vec<T>");
889 
890  Vec<T> r(m.n_rows);
891  r.clear();
892 
893  for (int c = 0; c < m.n_cols; c++) {
894  for (int p = 0; p < m.col[c].nnz(); p++)
895  r(m.col[c].get_nz_index(p)) += m.col[c].get_nz_data(p) * v(c);
896  }
897 
898  return r;
899 }
900 
901 template <class T>
903 {
904  it_assert_debug(v.size() == m.n_rows, "Vec<T> * Sparse_Mat<T>");
905 
906  Vec<T> r(m.n_cols);
907  r.clear();
908 
909  for (int c = 0; c < m.n_cols; c++)
910  r[c] = v * m.col[c];
911 
912  return r;
913 }
914 
915 template <class T>
917 {
918  Mat<T> ret(m.n_cols, m.n_cols);
919  Vec<T> col;
920  for (int c = 0; c < ret.cols(); c++) {
921  m.col[c].full(col);
922  for (int r = 0; r < c; r++) {
923  T tmp = m.col[r] * col;
924  ret(r, c) = tmp;
925  ret(c, r) = tmp;
926  }
927  ret(c, c) = m.col[c].sqr();
928  }
929 
930  return ret;
931 }
932 
933 template <class T>
935 {
936  Sparse_Mat<T> ret(m.n_cols, m.n_cols);
937  Vec<T> col;
938  T tmp;
939  for (int c = 0; c < ret.n_cols; c++) {
940  m.col[c].full(col);
941  for (int r = 0; r < c; r++) {
942  tmp = m.col[r] * col;
943  if (tmp != T(0)) {
944  ret.col[c].set_new(r, tmp);
945  ret.col[r].set_new(c, tmp);
946  }
947  }
948  tmp = m.col[c].sqr();
949  if (tmp != T(0))
950  ret.col[c].set_new(c, tmp);
951  }
952 
953  return ret;
954 }
955 
956 template <class T>
958 {
959  it_assert_debug(m1.n_rows == m2.n_rows, "trans_mult()");
960 
961  Sparse_Mat<T> ret(m1.n_cols, m2.n_cols);
962  Vec<T> col;
963  for (int c = 0; c < ret.n_cols; c++) {
964  m2.col[c].full(col);
965  for (int r = 0; r < ret.n_rows; r++)
966  ret.col[c].set_new(r, m1.col[r] * col);
967  }
968 
969  return ret;
970 }
971 
972 template <class T>
974 {
975  Vec<T> r(m.n_cols);
976  for (int c = 0; c < m.n_cols; c++)
977  r(c) = m.col[c] * v;
978 
979  return r;
980 }
981 
982 template <class T>
984 {
985  return trans_mult(m1.transpose(), m2.transpose());
986 }
987 
989 template <class T>
990 inline Sparse_Mat<T> sparse(const Mat<T> &m, T epsilon)
991 {
992  Sparse_Mat<T> s(m, epsilon);
993  return s;
994 }
995 
997 template <class T>
998 inline Mat<T> full(const Sparse_Mat<T> &s)
999 {
1000  Mat<T> m;
1001  s.full(m);
1002  return m;
1003 }
1004 
1006 template <class T>
1008 {
1009  Sparse_Mat<T> m;
1010  s.transpose(m);
1011  return m;
1012 }
1013 
1015 
1016 // ---------------------------------------------------------------------
1017 // Instantiations
1018 // ---------------------------------------------------------------------
1019 
1020 ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Sparse_Mat<int>;
1021 ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Sparse_Mat<double>;
1022 ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Sparse_Mat<std::complex<double> >;
1023 
1024 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_imat operator+(const sparse_imat &, const sparse_imat &);
1025 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_mat operator+(const sparse_mat &, const sparse_mat &);
1026 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_cmat operator+(const sparse_cmat &, const sparse_cmat &);
1027 
1028 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_imat operator*(const sparse_imat &, const sparse_imat &);
1029 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_mat operator*(const sparse_mat &, const sparse_mat &);
1030 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_cmat operator*(const sparse_cmat &, const sparse_cmat &);
1031 
1032 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec operator*(const ivec &, const sparse_imat &);
1033 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec operator*(const vec &, const sparse_mat &);
1034 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec operator*(const cvec &, const sparse_cmat &);
1035 
1036 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec operator*(const sparse_imat &, const ivec &);
1037 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec operator*(const sparse_mat &, const vec &);
1038 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec operator*(const sparse_cmat &, const cvec &);
1039 
1040 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT imat trans_mult(const sparse_imat &);
1041 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT mat trans_mult(const sparse_mat &);
1042 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cmat trans_mult(const sparse_cmat &);
1043 
1044 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_imat trans_mult_s(const sparse_imat &);
1045 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_mat trans_mult_s(const sparse_mat &);
1046 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_cmat trans_mult_s(const sparse_cmat &);
1047 
1048 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_imat trans_mult(const sparse_imat &, const sparse_imat &);
1049 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_mat trans_mult(const sparse_mat &, const sparse_mat &);
1050 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_cmat trans_mult(const sparse_cmat &, const sparse_cmat &);
1051 
1052 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec trans_mult(const sparse_imat &, const ivec &);
1053 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec trans_mult(const sparse_mat &, const vec &);
1054 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec trans_mult(const sparse_cmat &, const cvec &);
1055 
1056 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_imat mult_trans(const sparse_imat &, const sparse_imat &);
1057 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_mat mult_trans(const sparse_mat &, const sparse_mat &);
1058 ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_cmat mult_trans(const sparse_cmat &, const sparse_cmat &);
1059 
1061 
1062 } // namespace itpp
1063 
1064 #endif // #ifndef SMAT_H
1065 
SourceForge Logo

Generated on Sat Jul 6 2013 10:54:21 for IT++ by Doxygen 1.8.2