IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
svec.h
Go to the documentation of this file.
1 
29 #ifndef SVEC_H
30 #define SVEC_H
31 
32 #include <itpp/base/vec.h>
33 #include <itpp/base/math/min_max.h>
34 #include <cstdlib>
35 
36 
37 namespace itpp
38 {
39 
40 // Declaration of class Vec
41 template <class T> class Vec;
42 // Declaration of class Sparse_Vec
43 template <class T> class Sparse_Vec;
44 
45 // ----------------------- Sparse_Vec Friends -------------------------------
46 
48 template <class T>
49 Sparse_Vec<T> operator+(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
50 
52 template <class T>
53 T operator*(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
54 
56 template <class T>
57 T operator*(const Sparse_Vec<T> &v1, const Vec<T> &v2);
58 
60 template <class T>
61 T operator*(const Vec<T> &v1, const Sparse_Vec<T> &v2);
62 
64 template <class T>
65 Sparse_Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
66 
68 template <class T>
69 Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Vec<T> &v2);
70 
72 template <class T>
73 Sparse_Vec<T> elem_mult_s(const Sparse_Vec<T> &v1, const Vec<T> &v2);
74 
76 template <class T>
77 Vec<T> elem_mult(const Vec<T> &v1, const Sparse_Vec<T> &v2);
78 
80 template <class T>
81 Sparse_Vec<T> elem_mult_s(const Vec<T> &v1, const Sparse_Vec<T> &v2);
82 
93 template <class T>
95 {
96 public:
97 
99  Sparse_Vec();
100 
107  Sparse_Vec(int sz, int data_init = 200);
108 
114  Sparse_Vec(const Sparse_Vec<T> &v);
115 
121  Sparse_Vec(const Vec<T> &v);
122 
128  Sparse_Vec(const Vec<T> &v, T epsilon);
129 
131  ~Sparse_Vec();
132 
139  void set_size(int sz, int data_init = -1);
140 
142  int size() const { return v_size; }
143 
145  inline int nnz() {
146  if (check_small_elems_flag) {
148  }
149  return used_size;
150  }
151 
153  double density();
154 
156  void set_small_element(const T& epsilon);
157 
163  void remove_small_elements();
164 
170  void resize_data(int new_size);
171 
173  void compact();
174 
176  void full(Vec<T> &v) const;
177 
179  Vec<T> full() const;
180 
182  T operator()(int i) const;
183 
185  void set(int i, T v);
186 
188  void set(const ivec &index_vec, const Vec<T> &v);
189 
191  void set_new(int i, T v);
192 
194  void set_new(const ivec &index_vec, const Vec<T> &v);
195 
197  void add_elem(const int i, const T v);
198 
200  void add(const ivec &index_vec, const Vec<T> &v);
201 
203  void zeros();
204 
206  void zero_elem(const int i);
207 
209  void clear();
210 
212  void clear_elem(const int i);
213 
217  inline void get_nz_data(int p, T& data_out) {
218  if (check_small_elems_flag) {
220  }
221  data_out = data[p];
222  }
223 
225  inline T get_nz_data(int p) {
226  if (check_small_elems_flag) {
228  }
229  return data[p];
230  }
231 
233  inline int get_nz_index(int p) {
234  if (check_small_elems_flag) {
236  }
237  return index[p];
238  }
239 
241  inline void get_nz(int p, int &idx, T &dat) {
242  if (check_small_elems_flag) {
244  }
245  idx = index[p];
246  dat = data[p];
247  }
248 
250  ivec get_nz_indices();
251 
253  Sparse_Vec<T> get_subvector(int i1, int i2) const;
254 
256  T sqr() const;
257 
259  void operator=(const Sparse_Vec<T> &v);
261  void operator=(const Vec<T> &v);
262 
264  Sparse_Vec<T> operator-() const;
265 
267  bool operator==(const Sparse_Vec<T> &v);
268 
270  void operator+=(const Sparse_Vec<T> &v);
271 
273  void operator+=(const Vec<T> &v);
274 
276  void operator-=(const Sparse_Vec<T> &v);
277 
279  void operator-=(const Vec<T> &v);
280 
282  void operator*=(const T &v);
283 
285  void operator/=(const T &v);
286 
288  friend Sparse_Vec<T> operator+<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
290  friend T operator*<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
292  friend T operator*<>(const Sparse_Vec<T> &v1, const Vec<T> &v2);
294  friend T operator*<>(const Vec<T> &v1, const Sparse_Vec<T> &v2);
295 
297  friend Sparse_Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
298 
300  friend Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Vec<T> &v2);
301 
303  friend Sparse_Vec<T> elem_mult_s <>(const Sparse_Vec<T> &v1, const Vec<T> &v2);
304 
306  friend Vec<T> elem_mult <>(const Vec<T> &v1, const Sparse_Vec<T> &v2);
307 
309  friend Sparse_Vec<T> elem_mult_s <>(const Vec<T> &v1, const Sparse_Vec<T> &v2);
310 
311 private:
312  void init();
313  void alloc();
314  void free();
315 
316  int v_size, used_size, data_size;
317  T *data;
318  int *index;
319  T eps;
320  bool check_small_elems_flag;
321 };
322 
323 
329 
335 
341 
342 // ----------------------- Implementation starts here --------------------------------
343 
344 template <class T>
345 void Sparse_Vec<T>::init()
346 {
347  v_size = 0;
348  used_size = 0;
349  data_size = 0;
350  data = 0;
351  index = 0;
352  eps = 0;
353  check_small_elems_flag = true;
354 }
355 
356 template <class T>
358 {
359  if (data_size != 0) {
360  data = new T[data_size];
361  index = new int[data_size];
362  }
363 }
364 
365 template <class T>
366 void Sparse_Vec<T>::free()
367 {
368  delete [] data;
369  data = 0;
370  delete [] index;
371  index = 0;
372 }
373 
374 template <class T>
376 {
377  init();
378 }
379 
380 template <class T>
381 Sparse_Vec<T>::Sparse_Vec(int sz, int data_init)
382 {
383  init();
384  v_size = sz;
385  used_size = 0;
386  data_size = data_init;
387  alloc();
388 }
389 
390 template <class T>
392 {
393  init();
394  v_size = v.v_size;
395  used_size = v.used_size;
396  data_size = v.data_size;
397  eps = v.eps;
398  check_small_elems_flag = v.check_small_elems_flag;
399  alloc();
400 
401  for (int i = 0; i < used_size; i++) {
402  data[i] = v.data[i];
403  index[i] = v.index[i];
404  }
405 }
406 
407 template <class T>
409 {
410  init();
411  v_size = v.size();
412  used_size = 0;
413  data_size = std::min(v.size(), 10000);
414  alloc();
415 
416  for (int i = 0; i < v_size; i++) {
417  if (v(i) != T(0)) {
418  if (used_size == data_size)
419  resize_data(data_size*2);
420  data[used_size] = v(i);
421  index[used_size] = i;
422  used_size++;
423  }
424  }
425  compact();
426 }
427 
428 template <class T>
429 Sparse_Vec<T>::Sparse_Vec(const Vec<T> &v, T epsilon)
430 {
431  init();
432  v_size = v.size();
433  used_size = 0;
434  data_size = std::min(v.size(), 10000);
435  eps = epsilon;
436  alloc();
437 
438  double e = std::abs(epsilon);
439  for (int i = 0; i < v_size; i++) {
440  if (std::abs(v(i)) > e) {
441  if (used_size == data_size)
442  resize_data(data_size*2);
443  data[used_size] = v(i);
444  index[used_size] = i;
445  used_size++;
446  }
447  }
448  compact();
449 }
450 
451 template <class T>
453 {
454  free();
455 }
456 
457 template <class T>
458 void Sparse_Vec<T>::set_size(int new_size, int data_init)
459 {
460  v_size = new_size;
461  used_size = 0;
462  if (data_init != -1) {
463  free();
464  data_size = data_init;
465  alloc();
466  }
467 }
468 
469 template <class T>
471 {
472  if (check_small_elems_flag) {
473  remove_small_elements();
474  }
475  //return static_cast<double>(used_size) / v_size;
476  return double(used_size) / v_size;
477 }
478 
479 template <class T>
480 void Sparse_Vec<T>::set_small_element(const T& epsilon)
481 {
482  eps = epsilon;
483  remove_small_elements();
484 }
485 
486 template <class T>
488 {
489  int i;
490  int nrof_removed_elements = 0;
491  double e;
492 
493  //Remove small elements
494  e = std::abs(eps);
495  for (i = 0;i < used_size;i++) {
496  if (std::abs(data[i]) <= e) {
497  nrof_removed_elements++;
498  }
499  else if (nrof_removed_elements > 0) {
500  data[i-nrof_removed_elements] = data[i];
501  index[i-nrof_removed_elements] = index[i];
502  }
503  }
504 
505  //Set new size after small elements have been removed
506  used_size -= nrof_removed_elements;
507 
508  //Set the flag to indicate that all small elements have been removed
509  check_small_elems_flag = false;
510 }
511 
512 
513 template <class T>
514 void Sparse_Vec<T>::resize_data(int new_size)
515 {
516  it_assert(new_size >= used_size, "Sparse_Vec<T>::resize_data(int new_size): New size is to small");
517 
518  if (new_size != data_size) {
519  if (new_size == 0)
520  free();
521  else {
522  T *tmp_data = data;
523  int *tmp_pos = index;
524  data_size = new_size;
525  alloc();
526  for (int p = 0; p < used_size; p++) {
527  data[p] = tmp_data[p];
528  index[p] = tmp_pos[p];
529  }
530  delete [] tmp_data;
531  delete [] tmp_pos;
532  }
533  }
534 }
535 
536 template <class T>
538 {
539  if (check_small_elems_flag) {
540  remove_small_elements();
541  }
542  resize_data(used_size);
543 }
544 
545 template <class T>
547 {
548  v.set_size(v_size);
549 
550  v = T(0);
551  for (int p = 0; p < used_size; p++)
552  v(index[p]) = data[p];
553 }
554 
555 template <class T>
557 {
558  Vec<T> r(v_size);
559  full(r);
560  return r;
561 }
562 
563 // This is slow. Implement a better search
564 template <class T>
566 {
567  it_assert_debug(i >= 0 && i < v_size, "The index of the element is out of range");
568 
569  bool found = false;
570  int p;
571  for (p = 0; p < used_size; p++) {
572  if (index[p] == i) {
573  found = true;
574  break;
575  }
576  }
577  return found ? data[p] : T(0);
578 }
579 
580 template <class T>
581 void Sparse_Vec<T>::set(int i, T v)
582 {
583  it_assert_debug(i >= 0 && i < v_size, "The index of the element is out of range");
584 
585  bool found = false;
586  bool larger_than_eps;
587  int p;
588 
589  for (p = 0; p < used_size; p++) {
590  if (index[p] == i) {
591  found = true;
592  break;
593  }
594  }
595 
596  larger_than_eps = (std::abs(v) > std::abs(eps));
597 
598  if (found && larger_than_eps)
599  data[p] = v;
600  else if (larger_than_eps) {
601  if (used_size == data_size)
602  resize_data(data_size*2 + 100);
603  data[used_size] = v;
604  index[used_size] = i;
605  used_size++;
606  }
607 
608  //Check if the stored element is smaller than eps. In that case it should be removed.
609  if (std::abs(v) <= std::abs(eps)) {
610  remove_small_elements();
611  }
612 
613 }
614 
615 template <class T>
616 void Sparse_Vec<T>::set_new(int i, T v)
617 {
618  it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector");
619 
620  //Check that the new element is larger than eps!
621  if (std::abs(v) > std::abs(eps)) {
622  if (used_size == data_size)
623  resize_data(data_size*2 + 100);
624  data[used_size] = v;
625  index[used_size] = i;
626  used_size++;
627  }
628 }
629 
630 template <class T>
631 void Sparse_Vec<T>::add_elem(const int i, const T v)
632 {
633  bool found = false;
634  int p;
635 
636  it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector");
637 
638  for (p = 0; p < used_size; p++) {
639  if (index[p] == i) {
640  found = true;
641  break;
642  }
643  }
644  if (found)
645  data[p] += v;
646  else {
647  if (used_size == data_size)
648  resize_data(data_size*2 + 100);
649  data[used_size] = v;
650  index[used_size] = i;
651  used_size++;
652  }
653 
654  check_small_elems_flag = true;
655 
656 }
657 
658 template <class T>
659 void Sparse_Vec<T>::add(const ivec& index_vec, const Vec<T>& v)
660 {
661  bool found = false;
662  int i, p, q;
663  int nrof_nz = v.size();
664 
665  it_assert_debug(v_size > max(index_vec), "The indices exceeds the size of the sparse vector");
666 
667  //Elements are added if they have identical indices
668  for (q = 0; q < nrof_nz; q++) {
669  i = index_vec(q);
670  for (p = 0; p < used_size; p++) {
671  if (index[p] == i) {
672  found = true;
673  break;
674  }
675  }
676  if (found)
677  data[p] += v(q);
678  else {
679  if (used_size == data_size)
680  resize_data(data_size*2 + 100);
681  data[used_size] = v(q);
682  index[used_size] = i;
683  used_size++;
684  }
685  found = false;
686  }
687 
688  check_small_elems_flag = true;
689 
690 }
691 
692 template <class T>
694 {
695  used_size = 0;
696  check_small_elems_flag = false;
697 }
698 
699 template <class T>
700 void Sparse_Vec<T>::zero_elem(const int i)
701 {
702  bool found = false;
703  int p;
704 
705  it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector");
706 
707  for (p = 0; p < used_size; p++) {
708  if (index[p] == i) {
709  found = true;
710  break;
711  }
712  }
713  if (found) {
714  data[p] = data[used_size-1];
715  index[p] = index[used_size-1];
716  used_size--;
717  }
718 }
719 
720 template <class T>
722 {
723  used_size = 0;
724  check_small_elems_flag = false;
725 }
726 
727 template <class T>
728 void Sparse_Vec<T>::clear_elem(const int i)
729 {
730  bool found = false;
731  int p;
732 
733  it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector");
734 
735  for (p = 0; p < used_size; p++) {
736  if (index[p] == i) {
737  found = true;
738  break;
739  }
740  }
741  if (found) {
742  data[p] = data[used_size-1];
743  index[p] = index[used_size-1];
744  used_size--;
745  }
746 }
747 
748 template <class T>
749 void Sparse_Vec<T>::set(const ivec& index_vec, const Vec<T>& v)
750 {
751  it_assert_debug(v_size > max(index_vec), "The indices exceeds the size of the sparse vector");
752 
753  //Clear all old non-zero elements
754  clear();
755 
756  //Add the new non-zero elements
757  add(index_vec, v);
758 }
759 
760 template <class T>
761 void Sparse_Vec<T>::set_new(const ivec& index_vec, const Vec<T>& v)
762 {
763  int q;
764  int nrof_nz = v.size();
765 
766  it_assert_debug(v_size > max(index_vec), "The indices exceeds the size of the sparse vector");
767 
768  //Clear all old non-zero elements
769  clear();
770 
771  for (q = 0; q < nrof_nz; q++) {
772  if (std::abs(v[q]) > std::abs(eps)) {
773  if (used_size == data_size)
774  resize_data(data_size*2 + 100);
775  data[used_size] = v(q);
776  index[used_size] = index_vec(q);
777  used_size++;
778  }
779  }
780 }
781 
782 template <class T>
784 {
785  int n = nnz();
786  ivec r(n);
787  for (int i = 0; i < n; i++) {
788  r(i) = get_nz_index(i);
789  }
790  return r;
791 }
792 
793 template <class T>
795 {
796  it_assert_debug(v_size > i1 && v_size > i2 && i1 <= i2 && i1 >= 0, "The index of the element exceeds the size of the sparse vector");
797 
798  Sparse_Vec<T> r(i2 - i1 + 1);
799 
800  for (int p = 0; p < used_size; p++) {
801  if (index[p] >= i1 && index[p] <= i2) {
802  if (r.used_size == r.data_size)
803  r.resize_data(r.data_size*2 + 100);
804  r.data[r.used_size] = data[p];
805  r.index[r.used_size] = index[p] - i1;
806  r.used_size++;
807  }
808  }
809  r.eps = eps;
810  r.check_small_elems_flag = check_small_elems_flag;
811  r.compact();
812 
813  return r;
814 }
815 
816 template <class T>
818 {
819  T sum(0);
820  for (int p = 0; p < used_size; p++)
821  sum += data[p] * data[p];
822 
823  return sum;
824 }
825 
826 template <class T>
828 {
829  free();
830  v_size = v.v_size;
831  used_size = v.used_size;
832  data_size = v.data_size;
833  eps = v.eps;
834  check_small_elems_flag = v.check_small_elems_flag;
835  alloc();
836 
837  for (int i = 0; i < used_size; i++) {
838  data[i] = v.data[i];
839  index[i] = v.index[i];
840  }
841 }
842 
843 template <class T>
845 {
846  free();
847  v_size = v.size();
848  used_size = 0;
849  data_size = std::min(v.size(), 10000);
850  eps = T(0);
851  check_small_elems_flag = false;
852  alloc();
853 
854  for (int i = 0; i < v_size; i++) {
855  if (v(i) != T(0)) {
856  if (used_size == data_size)
857  resize_data(data_size*2);
858  data[used_size] = v(i);
859  index[used_size] = i;
860  used_size++;
861  }
862  }
863  compact();
864 }
865 
866 template <class T>
868 {
869  Sparse_Vec r(v_size, used_size);
870 
871  for (int p = 0; p < used_size; p++) {
872  r.data[p] = -data[p];
873  r.index[p] = index[p];
874  }
875  r.used_size = used_size;
876 
877  return r;
878 }
879 
880 template <class T>
882 {
883  int p, q;
884  bool found = false;
885 
886  //Remove small elements before comparing the two sparse_vectors
887  if (check_small_elems_flag)
888  remove_small_elements();
889 
890  if (v_size != v.v_size) {
891  //Return false if vector sizes are unequal
892  return false;
893  }
894  else {
895  for (p = 0;p < used_size;p++) {
896  for (q = 0;q < v.used_size;q++) {
897  if (index[p] == v.index[q]) {
898  found = true;
899  break;
900  }
901  }
902  if (found == false)
903  //Return false if non-zero element not found, or if elements are unequal
904  return false;
905  else if (data[p] != v.data[q])
906  //Return false if non-zero element not found, or if elements are unequal
907  return false;
908  else
909  //Check next non-zero element
910  found = false;
911  }
912  }
913 
914  /*Special handling if sizes do not match.
915  Required since v may need to do remove_small_elements() for true comparison*/
916  if (used_size != v.used_size) {
917  if (used_size > v.used_size) {
918  //Return false if number of non-zero elements is less in v
919  return false;
920  }
921  else {
922  //Ensure that the remaining non-zero elements in v are smaller than v.eps
923  int nrof_small_elems = 0;
924  for (q = 0;q < v.used_size;q++) {
925  if (std::abs(v.data[q]) <= std::abs(v.eps))
926  nrof_small_elems++;
927  }
928  if (v.used_size - nrof_small_elems != used_size)
929  //Return false if the number of "true" non-zero elements are unequal
930  return false;
931  }
932  }
933 
934  //All elements checks => return true
935  return true;
936 }
937 
938 template <class T>
940 {
941  int i, p;
942  T tmp_data;
943  int nrof_nz_v = v.used_size;
944 
945  it_assert_debug(v_size == v.size(), "Attempted addition of unequal sized sparse vectors");
946 
947  for (p = 0; p < nrof_nz_v; p++) {
948  i = v.index[p];
949  tmp_data = v.data[p];
950  //get_nz(p,i,tmp_data);
951  add_elem(i, tmp_data);
952  }
953 
954  check_small_elems_flag = true;
955 }
956 
957 template <class T>
959 {
960  int i;
961 
962  it_assert_debug(v_size == v.size(), "Attempted addition of unequal sized sparse vectors");
963 
964  for (i = 0; i < v.size(); i++)
965  if (v(i) != T(0))
966  add_elem(i, v(i));
967 
968  check_small_elems_flag = true;
969 }
970 
971 
972 template <class T>
974 {
975  int i, p;
976  T tmp_data;
977  int nrof_nz_v = v.used_size;
978 
979  it_assert_debug(v_size == v.size(), "Attempted subtraction of unequal sized sparse vectors");
980 
981  for (p = 0; p < nrof_nz_v; p++) {
982  i = v.index[p];
983  tmp_data = v.data[p];
984  //v.get_nz(p,i,tmp_data);
985  add_elem(i, -tmp_data);
986  }
987 
988  check_small_elems_flag = true;
989 }
990 
991 template <class T>
993 {
994  int i;
995 
996  it_assert_debug(v_size == v.size(), "Attempted subtraction of unequal sized sparse vectors");
997 
998  for (i = 0; i < v.size(); i++)
999  if (v(i) != T(0))
1000  add_elem(i, -v(i));
1001 
1002  check_small_elems_flag = true;
1003 }
1004 
1005 template <class T>
1007 {
1008  int p;
1009 
1010  for (p = 0; p < used_size; p++) {
1011  data[p] *= v;
1012  }
1013 
1014  check_small_elems_flag = true;
1015 }
1016 
1017 template <class T>
1019 {
1020  int p;
1021  for (p = 0; p < used_size; p++) {
1022  data[p] /= v;
1023  }
1024 
1025  if (std::abs(eps) > 0) {
1026  check_small_elems_flag = true;
1027  }
1028 }
1029 
1030 template <class T>
1031 T operator*(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2)
1032 {
1033  it_assert_debug(v1.v_size == v2.v_size, "Sparse_Vec<T> * Sparse_Vec<T>");
1034 
1035  T sum(0);
1036  Vec<T> v1f(v1.v_size);
1037  v1.full(v1f);
1038  for (int p = 0; p < v2.used_size; p++) {
1039  if (v1f[v2.index[p]] != T(0))
1040  sum += v1f[v2.index[p]] * v2.data[p];
1041  }
1042 
1043  return sum;
1044 }
1045 
1046 template <class T>
1047 T operator*(const Sparse_Vec<T> &v1, const Vec<T> &v2)
1048 {
1049  it_assert_debug(v1.size() == v2.size(), "Multiplication of unequal sized vectors attempted");
1050 
1051  T sum(0);
1052  for (int p1 = 0; p1 < v1.used_size; p1++)
1053  sum += v1.data[p1] * v2[v1.index[p1]];
1054 
1055  return sum;
1056 }
1057 
1058 template <class T>
1059 T operator*(const Vec<T> &v1, const Sparse_Vec<T> &v2)
1060 {
1061  it_assert_debug(v1.size() == v2.size(), "Multiplication of unequal sized vectors attempted");
1062 
1063  T sum(0);
1064  for (int p2 = 0; p2 < v2.used_size; p2++)
1065  sum += v1[v2.index[p2]] * v2.data[p2];
1066 
1067  return sum;
1068 }
1069 
1070 template <class T>
1072 {
1073  it_assert_debug(v1.v_size == v2.v_size, "elem_mult(Sparse_Vec<T>, Sparse_Vec<T>)");
1074 
1075  Sparse_Vec<T> r(v1.v_size);
1076  ivec pos(v1.v_size);
1077  pos = -1;
1078  for (int p1 = 0; p1 < v1.used_size; p1++)
1079  pos[v1.index[p1]] = p1;
1080  for (int p2 = 0; p2 < v2.used_size; p2++) {
1081  if (pos[v2.index[p2]] != -1) {
1082  if (r.used_size == r.data_size)
1083  r.resize_data(r.used_size*2 + 100);
1084  r.data[r.used_size] = v1.data[pos[v2.index[p2]]] * v2.data[p2];
1085  r.index[r.used_size] = v2.index[p2];
1086  r.used_size++;
1087  }
1088  }
1089  r.compact();
1090 
1091  return r;
1092 }
1093 
1094 template <class T>
1095 Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Vec<T> &v2)
1096 {
1097  it_assert_debug(v1.v_size == v2.size(), "elem_mult(Sparse_Vec<T>, Vec<T>)");
1098 
1099  Vec<T> r(v1.v_size);
1100  r = T(0);
1101  for (int p1 = 0; p1 < v1.used_size; p1++)
1102  r[v1.index[p1]] = v1.data[p1] * v2[v1.index[p1]];
1103 
1104  return r;
1105 }
1106 
1107 template <class T>
1109 {
1110  it_assert_debug(v1.v_size == v2.size(), "elem_mult(Sparse_Vec<T>, Vec<T>)");
1111 
1112  Sparse_Vec<T> r(v1.v_size);
1113  for (int p1 = 0; p1 < v1.used_size; p1++) {
1114  if (v2[v1.index[p1]] != T(0)) {
1115  if (r.used_size == r.data_size)
1116  r.resize_data(r.used_size*2 + 100);
1117  r.data[r.used_size] = v1.data[p1] * v2[v1.index[p1]];
1118  r.index[r.used_size] = v1.index[p1];
1119  r.used_size++;
1120  }
1121  }
1122  r.compact();
1123 
1124  return r;
1125 }
1126 
1127 template <class T>
1128 Vec<T> elem_mult(const Vec<T> &v1, const Sparse_Vec<T> &v2)
1129 {
1130  it_assert_debug(v1.size() == v2.v_size, "elem_mult(Vec<T>, Sparse_Vec<T>)");
1131 
1132  Vec<T> r(v2.v_size);
1133  r = T(0);
1134  for (int p2 = 0; p2 < v2.used_size; p2++)
1135  r[v2.index[p2]] = v1[v2.index[p2]] * v2.data[p2];
1136 
1137  return r;
1138 }
1139 
1140 template <class T>
1142 {
1143  it_assert_debug(v1.size() == v2.v_size, "elem_mult(Vec<T>, Sparse_Vec<T>)");
1144 
1145  Sparse_Vec<T> r(v2.v_size);
1146  for (int p2 = 0; p2 < v2.used_size; p2++) {
1147  if (v1[v2.index[p2]] != T(0)) {
1148  if (r.used_size == r.data_size)
1149  r.resize_data(r.used_size*2 + 100);
1150  r.data[r.used_size] = v1[v2.index[p2]] * v2.data[p2];
1151  r.index[r.used_size] = v2.index[p2];
1152  r.used_size++;
1153  }
1154  }
1155  r.compact();
1156 
1157  return r;
1158 }
1159 
1160 template <class T>
1162 {
1163  it_assert_debug(v1.v_size == v2.v_size, "Sparse_Vec<T> + Sparse_Vec<T>");
1164 
1165  Sparse_Vec<T> r(v1);
1166  ivec pos(v1.v_size);
1167  pos = -1;
1168  for (int p1 = 0; p1 < v1.used_size; p1++)
1169  pos[v1.index[p1]] = p1;
1170  for (int p2 = 0; p2 < v2.used_size; p2++) {
1171  if (pos[v2.index[p2]] == -1) {// A new entry
1172  if (r.used_size == r.data_size)
1173  r.resize_data(r.used_size*2 + 100);
1174  r.data[r.used_size] = v2.data[p2];
1175  r.index[r.used_size] = v2.index[p2];
1176  r.used_size++;
1177  }
1178  else
1179  r.data[pos[v2.index[p2]]] += v2.data[p2];
1180  }
1181  r.check_small_elems_flag = true; // added dec 7, 2006
1182  r.compact();
1183 
1184  return r;
1185 }
1186 
1188 template <class T>
1189 inline Sparse_Vec<T> sparse(const Vec<T> &v)
1190 {
1191  Sparse_Vec<T> s(v);
1192  return s;
1193 }
1194 
1196 template <class T>
1197 inline Sparse_Vec<T> sparse(const Vec<T> &v, T epsilon)
1198 {
1199  Sparse_Vec<T> s(v, epsilon);
1200  return s;
1201 }
1202 
1204 template <class T>
1205 inline Vec<T> full(const Sparse_Vec<T> &s)
1206 {
1207  Vec<T> v;
1208  s.full(v);
1209  return v;
1210 }
1211 
1213 
1214 // ---------------------------------------------------------------------
1215 // Instantiations
1216 // ---------------------------------------------------------------------
1217 
1218 #ifndef _MSC_VER
1219 
1220 extern template class Sparse_Vec<int>;
1221 extern template class Sparse_Vec<double>;
1222 extern template class Sparse_Vec<std::complex<double> >;
1223 
1224 extern template sparse_ivec operator+(const sparse_ivec &,
1225  const sparse_ivec &);
1226 extern template sparse_vec operator+(const sparse_vec &,
1227  const sparse_vec &);
1228 extern template sparse_cvec operator+(const sparse_cvec &,
1229  const sparse_cvec &);
1230 
1231 extern template int operator*(const sparse_ivec &, const sparse_ivec &);
1232 extern template double operator*(const sparse_vec &, const sparse_vec &);
1233 extern template std::complex<double> operator*(const sparse_cvec &,
1234  const sparse_cvec &);
1235 
1236 extern template int operator*(const sparse_ivec &, const ivec &);
1237 extern template double operator*(const sparse_vec &, const vec &);
1238 extern template std::complex<double> operator*(const sparse_cvec &,
1239  const cvec &);
1240 
1241 extern template int operator*(const ivec &, const sparse_ivec &);
1242 extern template double operator*(const vec &, const sparse_vec &);
1243 extern template std::complex<double> operator*(const cvec &,
1244  const sparse_cvec &);
1245 
1246 extern template sparse_ivec elem_mult(const sparse_ivec &,
1247  const sparse_ivec &);
1248 extern template sparse_vec elem_mult(const sparse_vec &, const sparse_vec &);
1249 extern template sparse_cvec elem_mult(const sparse_cvec &,
1250  const sparse_cvec &);
1251 
1252 extern template ivec elem_mult(const sparse_ivec &, const ivec &);
1253 extern template vec elem_mult(const sparse_vec &, const vec &);
1254 extern template cvec elem_mult(const sparse_cvec &, const cvec &);
1255 
1256 extern template sparse_ivec elem_mult_s(const sparse_ivec &, const ivec &);
1257 extern template sparse_vec elem_mult_s(const sparse_vec &, const vec &);
1258 extern template sparse_cvec elem_mult_s(const sparse_cvec &, const cvec &);
1259 
1260 extern template ivec elem_mult(const ivec &, const sparse_ivec &);
1261 extern template vec elem_mult(const vec &, const sparse_vec &);
1262 extern template cvec elem_mult(const cvec &, const sparse_cvec &);
1263 
1264 extern template sparse_ivec elem_mult_s(const ivec &, const sparse_ivec &);
1265 extern template sparse_vec elem_mult_s(const vec &, const sparse_vec &);
1266 extern template sparse_cvec elem_mult_s(const cvec &, const sparse_cvec &);
1267 
1268 #endif // _MSC_VER
1269 
1271 
1272 } // namespace itpp
1273 
1274 #endif // #ifndef SVEC_H
1275 
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:13 for IT++ by Doxygen 1.8.1.2