41 template <
class T>
class Vec;
43 template <
class T>
class Sparse_Vec;
49 Sparse_Vec<T>
operator+(
const Sparse_Vec<T> &v1,
const Sparse_Vec<T> &v2);
53 T
operator*(
const Sparse_Vec<T> &v1,
const Sparse_Vec<T> &v2);
57 T
operator*(
const Sparse_Vec<T> &v1,
const Vec<T> &v2);
61 T
operator*(
const Vec<T> &v1,
const Sparse_Vec<T> &v2);
65 Sparse_Vec<T>
elem_mult(
const Sparse_Vec<T> &v1,
const Sparse_Vec<T> &v2);
69 Vec<T>
elem_mult(
const Sparse_Vec<T> &v1,
const Vec<T> &v2);
73 Sparse_Vec<T>
elem_mult_s(
const Sparse_Vec<T> &v1,
const Vec<T> &v2);
77 Vec<T>
elem_mult(
const Vec<T> &v1,
const Sparse_Vec<T> &v2);
81 Sparse_Vec<T>
elem_mult_s(
const Vec<T> &v1,
const Sparse_Vec<T> &v2);
139 void set_size(
int sz,
int data_init = -1);
142 int size()
const {
return v_size; }
146 if (check_small_elems_flag) {
185 void set(
int i, T v);
188 void set(
const ivec &index_vec,
const Vec<T> &v);
197 void add_elem(
const int i,
const T v);
200 void add(
const ivec &index_vec,
const Vec<T> &v);
218 if (check_small_elems_flag) {
226 if (check_small_elems_flag) {
234 if (check_small_elems_flag) {
241 inline void get_nz(
int p,
int &idx, T &dat) {
242 if (check_small_elems_flag) {
316 int v_size, used_size, data_size;
320 bool check_small_elems_flag;
353 check_small_elems_flag =
true;
359 if (data_size != 0) {
360 data =
new T[data_size];
361 index =
new int[data_size];
366 void Sparse_Vec<T>::free()
386 data_size = data_init;
395 used_size = v.used_size;
396 data_size = v.data_size;
398 check_small_elems_flag = v.check_small_elems_flag;
401 for (
int i = 0; i < used_size; i++) {
403 index[i] = v.index[i];
416 for (
int i = 0; i < v_size; i++) {
418 if (used_size == data_size)
419 resize_data(data_size*2);
420 data[used_size] = v(i);
421 index[used_size] = i;
439 for (
int i = 0; i < v_size; i++) {
441 if (used_size == data_size)
442 resize_data(data_size*2);
443 data[used_size] = v(i);
444 index[used_size] = i;
462 if (data_init != -1) {
464 data_size = data_init;
472 if (check_small_elems_flag) {
473 remove_small_elements();
476 return double(used_size) / v_size;
483 remove_small_elements();
490 int nrof_removed_elements = 0;
495 for (i = 0;i < used_size;i++) {
497 nrof_removed_elements++;
499 else if (nrof_removed_elements > 0) {
500 data[i-nrof_removed_elements] = data[i];
501 index[i-nrof_removed_elements] = index[i];
506 used_size -= nrof_removed_elements;
509 check_small_elems_flag =
false;
516 it_assert(new_size >= used_size,
"Sparse_Vec<T>::resize_data(int new_size): New size is to small");
518 if (new_size != data_size) {
523 int *tmp_pos = index;
524 data_size = new_size;
526 for (
int p = 0; p < used_size; p++) {
527 data[p] = tmp_data[p];
528 index[p] = tmp_pos[p];
539 if (check_small_elems_flag) {
540 remove_small_elements();
542 resize_data(used_size);
551 for (
int p = 0; p < used_size; p++)
552 v(index[p]) = data[p];
567 it_assert_debug(i >= 0 && i < v_size,
"The index of the element is out of range");
571 for (p = 0; p < used_size; p++) {
577 return found ? data[p] : T(0);
583 it_assert_debug(i >= 0 && i < v_size,
"The index of the element is out of range");
586 bool larger_than_eps;
589 for (p = 0; p < used_size; p++) {
598 if (found && larger_than_eps)
600 else if (larger_than_eps) {
601 if (used_size == data_size)
602 resize_data(data_size*2 + 100);
604 index[used_size] = i;
610 remove_small_elements();
618 it_assert_debug(v_size > i,
"The index of the element exceeds the size of the sparse vector");
622 if (used_size == data_size)
623 resize_data(data_size*2 + 100);
625 index[used_size] = i;
636 it_assert_debug(v_size > i,
"The index of the element exceeds the size of the sparse vector");
638 for (p = 0; p < used_size; p++) {
647 if (used_size == data_size)
648 resize_data(data_size*2 + 100);
650 index[used_size] = i;
654 check_small_elems_flag =
true;
663 int nrof_nz = v.
size();
665 it_assert_debug(v_size >
max(index_vec),
"The indices exceeds the size of the sparse vector");
668 for (q = 0; q < nrof_nz; q++) {
670 for (p = 0; p < used_size; p++) {
679 if (used_size == data_size)
680 resize_data(data_size*2 + 100);
681 data[used_size] = v(q);
682 index[used_size] = i;
688 check_small_elems_flag =
true;
696 check_small_elems_flag =
false;
705 it_assert_debug(v_size > i,
"The index of the element exceeds the size of the sparse vector");
707 for (p = 0; p < used_size; p++) {
714 data[p] = data[used_size-1];
715 index[p] = index[used_size-1];
724 check_small_elems_flag =
false;
733 it_assert_debug(v_size > i,
"The index of the element exceeds the size of the sparse vector");
735 for (p = 0; p < used_size; p++) {
742 data[p] = data[used_size-1];
743 index[p] = index[used_size-1];
751 it_assert_debug(v_size >
max(index_vec),
"The indices exceeds the size of the sparse vector");
764 int nrof_nz = v.
size();
766 it_assert_debug(v_size >
max(index_vec),
"The indices exceeds the size of the sparse vector");
771 for (q = 0; q < nrof_nz; q++) {
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);
787 for (
int i = 0; i < n; i++) {
788 r(i) = get_nz_index(i);
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");
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)
804 r.data[r.used_size] = data[p];
805 r.index[r.used_size] = index[p] - i1;
810 r.check_small_elems_flag = check_small_elems_flag;
820 for (
int p = 0; p < used_size; p++)
821 sum += data[p] * data[p];
831 used_size = v.used_size;
832 data_size = v.data_size;
834 check_small_elems_flag = v.check_small_elems_flag;
837 for (
int i = 0; i < used_size; i++) {
839 index[i] = v.index[i];
851 check_small_elems_flag =
false;
854 for (
int i = 0; i < v_size; i++) {
856 if (used_size == data_size)
857 resize_data(data_size*2);
858 data[used_size] = v(i);
859 index[used_size] = i;
871 for (
int p = 0; p < used_size; p++) {
872 r.data[p] = -data[p];
873 r.index[p] = index[p];
875 r.used_size = used_size;
887 if (check_small_elems_flag)
888 remove_small_elements();
890 if (v_size != v.v_size) {
895 for (p = 0;p < used_size;p++) {
896 for (q = 0;q < v.used_size;q++) {
897 if (index[p] == v.index[q]) {
905 else if (data[p] != v.data[q])
916 if (used_size != v.used_size) {
917 if (used_size > v.used_size) {
923 int nrof_small_elems = 0;
924 for (q = 0;q < v.used_size;q++) {
928 if (v.used_size - nrof_small_elems != used_size)
943 int nrof_nz_v = v.used_size;
947 for (p = 0; p < nrof_nz_v; p++) {
949 tmp_data = v.data[p];
951 add_elem(i, tmp_data);
954 check_small_elems_flag =
true;
964 for (i = 0; i < v.
size(); i++)
968 check_small_elems_flag =
true;
977 int nrof_nz_v = v.used_size;
979 it_assert_debug(v_size == v.
size(),
"Attempted subtraction of unequal sized sparse vectors");
981 for (p = 0; p < nrof_nz_v; p++) {
983 tmp_data = v.data[p];
985 add_elem(i, -tmp_data);
988 check_small_elems_flag =
true;
996 it_assert_debug(v_size == v.
size(),
"Attempted subtraction of unequal sized sparse vectors");
998 for (i = 0; i < v.
size(); i++)
1002 check_small_elems_flag =
true;
1010 for (p = 0; p < used_size; p++) {
1014 check_small_elems_flag =
true;
1021 for (p = 0; p < used_size; p++) {
1026 check_small_elems_flag =
true;
1033 it_assert_debug(v1.v_size == v2.v_size,
"Sparse_Vec<T> * Sparse_Vec<T>");
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];
1052 for (
int p1 = 0; p1 < v1.used_size; p1++)
1053 sum += v1.data[p1] * v2[v1.index[p1]];
1064 for (
int p2 = 0; p2 < v2.used_size; p2++)
1065 sum += v1[v2.index[p2]] * v2.data[p2];
1073 it_assert_debug(v1.v_size == v2.v_size,
"elem_mult(Sparse_Vec<T>, Sparse_Vec<T>)");
1076 ivec pos(v1.v_size);
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)
1084 r.data[r.used_size] = v1.data[pos[v2.index[p2]]] * v2.data[p2];
1085 r.index[r.used_size] = v2.index[p2];
1101 for (
int p1 = 0; p1 < v1.used_size; p1++)
1102 r[v1.index[p1]] = v1.data[p1] * v2[v1.index[p1]];
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];
1134 for (
int p2 = 0; p2 < v2.used_size; p2++)
1135 r[v2.index[p2]] = v1[v2.index[p2]] * v2.data[p2];
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];
1163 it_assert_debug(v1.v_size == v2.v_size,
"Sparse_Vec<T> + Sparse_Vec<T>");
1166 ivec pos(v1.v_size);
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) {
1172 if (r.used_size == r.data_size)
1174 r.data[r.used_size] = v2.data[p2];
1175 r.index[r.used_size] = v2.index[p2];
1179 r.data[pos[v2.index[p2]]] += v2.data[p2];
1181 r.check_small_elems_flag =
true;
1220 extern template class Sparse_Vec<int>;
1221 extern template class Sparse_Vec<double>;
1222 extern template class Sparse_Vec<std::complex<double> >;
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 &);
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 &);
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 &,
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 &);
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 &);
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 &);
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 &);
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 &);
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 &);
1274 #endif // #ifndef SVEC_H