/*---------------------------------------------------------------------------*
* IT++ *
*---------------------------------------------------------------------------*
* Copyright (c) 1995-2004 by Tony Ottosson, Thomas Eriksson, Pål Frenger, *
* Tobias Ringström, and Jonas Samuelsson. *
* *
* Permission to use, copy, modify, and distribute this software and its *
* documentation under the terms of the GNU General Public License is hereby *
* granted. No representations are made about the suitability of this *
* software for any purpose. It is provided "as is" without expressed or *
* implied warranty. See the GNU General Public License for more details. *
*---------------------------------------------------------------------------*/
/*!
\file
\brief Sparse Vector Class Definitions
\author Tony Ottosson, and Tobias Ringström
1.17
2003/06/30 07:42:44
*/
#ifndef __itpp_svec_h
#define __itpp_svec_h
#include <algorithm>
#include "itconfig.h"
#include "base/vec.h"
#include "base/matfunc.h"
#include "base/stat.h"
namespace itpp {
// Declaration of class Vec
template <class T> class Vec;
// Declaration of class Sparse_Vec
template <class T> class Sparse_Vec;
// ----------------------- Sparse_Vec Friends -------------------------------
//! v1+v2 where v1 and v2 are sparse vector
template <class T>
Sparse_Vec<T> operator+(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
//! v1*v2 where v1 and v2 are sparse vectors
template <class T>
T operator*(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
//! v1*v2 where v1 is a sparse vector and v2 is a dense vector
template <class T>
T operator*(const Sparse_Vec<T> &v1, const Vec<T> &v2);
//! v1*v2 where v1 is a dense vector and v2 is a sparse vector
template <class T>
T operator*(const Vec<T> &v1, const Sparse_Vec<T> &v2);
//! Elementwise multiplication of two sparse vectors returning a sparse vector
template <class T>
Sparse_Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
//! Elementwise multiplication of a sparse vector and a dense vector returning a dense vector
template <class T>
Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Vec<T> &v2);
//! Elementwise multiplication of a sparse vector and a dense vector returning a sparse vector
template <class T>
Sparse_Vec<T> elem_mult_s(const Sparse_Vec<T> &v1, const Vec<T> &v2);
//! Elementwise multiplication of a dense vector and a sparse vector returning a dense vector
template <class T>
Vec<T> elem_mult(const Vec<T> &v1, const Sparse_Vec<T> &v2);
//! Elementwise multiplication of a dense vector and a sparse vector returning a sparse vector
template <class T>
Sparse_Vec<T> elem_mult_s(const Vec<T> &v1, const Sparse_Vec<T> &v2);
/*!
\brief Templated sparse vector class
\author Tony Ottosson and Tobias Ringstrom
*/
template <class T>
class Sparse_Vec {
public:
//! Default constructor
Sparse_Vec();
/*!
\brief Initiate an empty sparse vector
\param sz Size of the sparse vector (i.e. maximum index is (\c sz - 1))
\param data_init Maximum number of non-zero elements in the sparse vector (default value 200)
*/
Sparse_Vec(int sz, int data_init=200);
/*!
\brief Initiate a new sparse vector.
\param v The elements of \c v are copied into the new sparse vector
*/
Sparse_Vec(const Sparse_Vec<T> &v);
/*!
\brief Initiate a new sparse vector from a dense vector.
\param v The elements of \c v are copied into the new sparse vector
*/
Sparse_Vec(const Vec<T> &v);
/*!
\brief Initiate a new sparse vector from a dense vector. Elements of \c v larger than \c epsilon are copied into the new sparse vector.
\note If the type T is complex<double>, then the elements of \c v larger than \c abs(epsilon) are copied into the new sparse vector.
*/
Sparse_Vec(const Vec<T> &v, T epsilon);
//! Destructor
~Sparse_Vec();
/*!
\brief Set the size \c sz of the sparse vector. Default value \c data_init=-1 \c => allocated size for the data is not changed.
\param sz Size of the sparse vector (i.e. maximum index is (\c sz - 1))
\param data_init Maximum number of non-zero elements in the sparse vector (default value -1 \c => allocated size for the data is not changed)
*/
void set_size(int sz, int data_init=-1);
//! Returns the size of the sparse vector
int size() const { return v_size; }
//! Number of non-zero elements in the sparse vector
inline int nnz()
{
if (check_small_elems_flag) {
remove_small_elements();
}
return used_size;
}
//! Returns the density of the sparse vector: (number of non-zero elements)/(size of the vector)
double density();
//! Set that all elements smaller than \a epsilon should be set to zero.
void set_small_element(const T& epsilon);
/*!
Removes all elements that are smaller than \a epsilon from the non-zero elements.
\note The small element \a epsilon can be set by the member function set_small_element. If no small value is set, the default value is always \c epsilon=0.
*/
void remove_small_elements();
/*!
\brief Set the maximum number of non-zero elements to \c new_size
\param new_size The new maximum number of non-zero elements.
*/
void resize_data(int new_size);
//! Set the maximum number of non-zero elements equal to the actual number of non-zero elements
void compact();
//! Returns a full, dense vector in \c v
void full(Vec<T> &v) const;
//! Returns a full, dense vector
Vec<T> full() const;
//! Returns the element with index \c i
T operator()(int i) const;
//! Set element \c i equal to \c v
void set(int i, T v);
//! Set the elements of the sparse vector with indices \c index_vec to the values in \c v
void set(const ivec &index_vec, const Vec<T> &v);
//! Set a new element with index \c i equal to \c v
void set_new(int i, T v);
//! Set new elements with indices \c index_vec equal to the values in \c v (no check whether the same index is used several times)
void set_new(const ivec &index_vec, const Vec<T> &v);
//! Add element \c i with \c v
void add_elem(const int i, const T v);
//! Add \c v to the elements specified by \c index_vec with \c v
void add(const ivec &index_vec, const Vec<T> &v);
//! Set the sparse vector to the all zero vector (removes all non-zero elements)
void zeros();
//! Set the i-th element to zero (i.e. clear that element if it contains a non-zero value)
void zero_elem(const int i);
//! Clear all non-zero elements of the sparse vector
void clear();
//! Clear the i-th element (if it contains a non-zero value)
void clear_elem(const int i);
/*!
\brief Extract the reference to the p-th non-zero data element
*/
inline void get_nz_data(int p, T& data_out)
{
if (check_small_elems_flag) {
remove_small_elements();
}
data_out = data[p];
}
//! Returns the p-th non-zero data element
inline T get_nz_data(int p)
{
if (check_small_elems_flag) {
remove_small_elements();
}
return data[p];
}
//! Returns the vector index of the p-th non-zero element
inline int get_nz_index(int p)
{
if (check_small_elems_flag) {
remove_small_elements();
}
return index[p];
}
//! Returns the p-th non-zero value in \c dat and the corresponding vector index in \c idx
inline void get_nz(int p, int &idx, T &dat)
{
if (check_small_elems_flag) {
remove_small_elements();
}
idx=index[p];
dat=data[p];
}
//! Return sparse subvector from index \c i1 to index \c i2
Sparse_Vec<T> get_subvector(int i1, int i2) const;
//! Returns the sum of all values squared
T sqr() const;
//! Assign sparse vector the value and length of the sparse vector \c v
void operator=(const Sparse_Vec<T> &v);
//! Assign sparse vector the value and length of the dense vector \c v
void operator=(const Vec<T> &v);
//! Returns the sign inverse of all elements in the sparse vector
Sparse_Vec<T> operator-() const;
//! Compare two sparse vectors. False if wrong sizes or different values
bool operator==(const Sparse_Vec<T> &v);
//! Add sparse vector \c v to all non-zero elements of the sparse vector
void operator+=(const Sparse_Vec<T> &v);
//! Add vector \c v to all non-zero elements of the sparse vector
void operator+=(const Vec<T> &v);
//! Subtract sparse vector \c v from all non-zero elements of the sparse vector
void operator-=(const Sparse_Vec<T> &v);
//! Subtract vector \c v from all non-zero elements of the sparse vector
void operator-=(const Vec<T> &v);
//! Multiply the scalar \c v to all non-zero elements of the sparse vector
void operator*=(const T &v);
//! Divide all non-zero elements of the sparse vector with the scalar \c v
void operator/=(const T &v);
//! Addition v1+v2 where v1 and v2 are sparse vector
friend Sparse_Vec<T> operator+<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
//! Scalar product v1*v2 where v1 and v2 are sparse vectors
friend T operator*<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
//! Scalar product v1*v2 where v1 is a sparse vector and v2 is a dense vector
friend T operator*<>(const Sparse_Vec<T> &v1, const Vec<T> &v2);
//! Scalar product v1*v2 where v1 is a dense vector and v2 is a sparse vector
friend T operator*<>(const Vec<T> &v1, const Sparse_Vec<T> &v2);
//! Element wise multiplication of two sparse vectors
friend Sparse_Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
//! Element wise multiplication of a sparse vector and a dense vector
friend Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Vec<T> &v2);
//! Element wise multiplication of a sparse vector and a dense vector returning a sparse vector
friend Sparse_Vec<T> elem_mult_s <>(const Sparse_Vec<T> &v1, const Vec<T> &v2);
//! Element wise multiplication of a a dense vector and a sparse vector
friend Vec<T> elem_mult <>(const Vec<T> &v1, const Sparse_Vec<T> &v2);
//! Element wise multiplication of a a dense vector and a sparse vector returning a sparse vector
friend Sparse_Vec<T> elem_mult_s <>(const Vec<T> &v1, const Sparse_Vec<T> &v2);
private:
void init();
void alloc();
void free();
int v_size, used_size, data_size;
T *data;
int *index;
T eps;
bool check_small_elems_flag;
};
/*!
\relates Sparse_Vec
\brief Type definition of an integer sparse vector
*/
typedef Sparse_Vec<int> sparse_ivec;
/*!
\relates Sparse_Vec
\brief Type definition of a double sparse vector
*/
typedef Sparse_Vec<double> sparse_vec;
/*!
\relates Sparse_Vec
\brief Type definition of a complex<double> sparse vector
*/
typedef Sparse_Vec<complex<double> > sparse_cvec;
// ----------------------- Implementation starts here --------------------------------
template <class T>
void Sparse_Vec<T>::init()
{
v_size = 0;
used_size = 0;
data_size = 0;
data = 0;
index = 0;
eps = 0;
check_small_elems_flag = true;
}
template <class T>
void Sparse_Vec<T>::alloc()
{
if (data_size != 0) {
data = new T[data_size];
index = new int[data_size];
}
}
template <class T>
void Sparse_Vec<T>::free()
{
delete [] data;
data = 0;
delete [] index;
index = 0;
}
template <class T>
Sparse_Vec<T>::Sparse_Vec()
{
init();
}
template <class T>
Sparse_Vec<T>::Sparse_Vec(int sz, int data_init)
{
init();
v_size = sz;
used_size = 0;
data_size = data_init;
alloc();
}
template <class T>
Sparse_Vec<T>::Sparse_Vec(const Sparse_Vec<T> &v)
{
init();
v_size = v.v_size;
used_size = v.used_size;
data_size = v.data_size;
eps = v.eps;
check_small_elems_flag = v.check_small_elems_flag;
alloc();
for (int i=0; i<used_size; i++) {
data[i] = v.data[i];
index[i] = v.index[i];
}
}
template <class T>
Sparse_Vec<T>::Sparse_Vec(const Vec<T> &v)
{
init();
v_size = v.size();
used_size = 0;
data_size = std::min(v.size(), 10000);
alloc();
for (int i=0; i<v_size; i++) {
if (v(i) != T(0)) {
if (used_size == data_size)
resize_data(data_size*2);
data[used_size] = v(i);
index[used_size] = i;
used_size++;
}
}
compact();
}
template <class T>
Sparse_Vec<T>::Sparse_Vec(const Vec<T> &v, T epsilon)
{
init();
v_size = v.size();
used_size = 0;
data_size = std::min(v.size(), 10000);
eps = epsilon;
alloc();
double e = std::abs(epsilon);
for (int i=0; i<v_size; i++) {
if (std::abs(v(i)) > e) {
if (used_size == data_size)
resize_data(data_size*2);
data[used_size] = v(i);
index[used_size] = i;
used_size++;
}
}
compact();
}
template <class T>
Sparse_Vec<T>::~Sparse_Vec()
{
free();
}
template <class T>
void Sparse_Vec<T>::set_size(int new_size, int data_init)
{
v_size = new_size;
used_size = 0;
if (data_init != -1){
free();
data_size = data_init;
alloc();
}
}
template <class T>
double Sparse_Vec<T>::density()
{
if (check_small_elems_flag) {
remove_small_elements();
}
//return static_cast<double>(used_size) / v_size;
return double(used_size) / v_size;
}
template <class T>
void Sparse_Vec<T>::set_small_element(const T& epsilon)
{
eps=epsilon;
remove_small_elements();
}
template <class T>
void Sparse_Vec<T>::remove_small_elements()
{
int i;
int nrof_removed_elements = 0;
double e;
//Remove small elements
e = std::abs(eps);
for (i=0;i<used_size;i++) {
if (std::abs(data[i]) <= e) {
nrof_removed_elements++;
}
else if (nrof_removed_elements > 0) {
data[i-nrof_removed_elements] = data[i];
index[i-nrof_removed_elements] = index[i];
}
}
//Set new size after small elements have been removed
used_size -= nrof_removed_elements;
//Set the flag to indicate that all small elements have been removed
check_small_elems_flag = false;
}
template <class T>
void Sparse_Vec<T>::resize_data(int new_size)
{
it_assert(new_size>=used_size, "Sparse_Vec<T>::resize_data(int new_size): New size is to small");
if (new_size != data_size) {
if (new_size == 0)
free();
else {
T *tmp_data = data;
int *tmp_pos = index;
data_size = new_size;
alloc();
for (int p=0; p<used_size; p++) {
data[p] = tmp_data[p];
index[p] = tmp_pos[p];
}
delete [] tmp_data;
delete [] tmp_pos;
}
}
}
template <class T>
void Sparse_Vec<T>::compact()
{
if (check_small_elems_flag) {
remove_small_elements();
}
resize_data(used_size);
}
template <class T>
void Sparse_Vec<T>::full(Vec<T> &v) const
{
v.set_size(v_size);
v = T(0);
for (int p=0; p<used_size; p++)
v(index[p]) = data[p];
}
template <class T>
Vec<T> Sparse_Vec<T>::full() const
{
Vec<T> r(v_size);
full(r);
return r;
}
template <class T>
T Sparse_Vec<T>::operator()(int i) const
{
it_assert0(v_size > i, "The index of the element exceeds the size of the sparse vector");
bool found = false;
int p;
for (p=0; p<used_size; p++) {
if (index[p] == i) {
found = true;
break;
}
}
return found ? data[p] : T(0);
}
template <class T>
void Sparse_Vec<T>::set(int i, T v)
{
it_assert0(v_size > i, "The index of the element exceeds the size of the sparse vector");
bool found = false;
bool larger_than_eps;
int p;
for (p=0; p<used_size; p++) {
if (index[p] == i) {
found = true;
break;
}
}
larger_than_eps = (std::abs(v) > std::abs(eps));
if (found && larger_than_eps)
data[p] = v;
else if (larger_than_eps) {
if (used_size == data_size)
resize_data(data_size*2+100);
data[used_size] = v;
index[used_size] = i;
used_size++;
}
//Check if the stored element is smaller than eps. In that case it should be removed.
if (std::abs(v) <= std::abs(eps)) {
remove_small_elements();
}
}
template <class T>
void Sparse_Vec<T>::set_new(int i, T v)
{
it_assert0(v_size > i, "The index of the element exceeds the size of the sparse vector");
//Check that the new element is larger than eps!
if (std::abs(v) > std::abs(eps)) {
if (used_size == data_size)
resize_data(data_size*2+100);
data[used_size] = v;
index[used_size] = i;
used_size++;
}
}
template <class T>
void Sparse_Vec<T>::add_elem(const int i, const T v)
{
bool found = false;
int p;
it_assert0(v_size > i, "The index of the element exceeds the size of the sparse vector");
for (p=0; p<used_size; p++) {
if (index[p] == i) {
found = true;
break;
}
}
if (found)
data[p] += v;
else {
if (used_size == data_size)
resize_data(data_size*2+100);
data[used_size] = v;
index[used_size] = i;
used_size++;
}
check_small_elems_flag = true;
}
template <class T>
void Sparse_Vec<T>::add(const ivec& index_vec, const Vec<T>& v)
{
bool found = false;
int i,p,q;
int nrof_nz=v.size();
it_assert0(v_size>max(index_vec),"The indices exceeds the size of the sparse vector");
//Elements are added if they have identical indices
for (q=0; q<nrof_nz; q++){
i=index_vec(q);
for (p=0; p<used_size; p++) {
if (index[p] == i) {
found = true;
break;
}
}
if (found)
data[p] += v(q);
else {
if (used_size == data_size)
resize_data(data_size*2+100);
data[used_size] = v(q);
index[used_size] = i;
used_size++;
}
found = false;
}
check_small_elems_flag = true;
}
template <class T>
void Sparse_Vec<T>::zeros()
{
used_size = 0;
check_small_elems_flag = false;
}
template <class T>
void Sparse_Vec<T>::zero_elem(const int i)
{
bool found = false;
int p;
it_assert0(v_size > i, "The index of the element exceeds the size of the sparse vector");
for (p=0; p<used_size; p++) {
if (index[p] == i) {
found = true;
break;
}
}
if (found) {
data[p] = data[used_size-1];
index[p] = index[used_size-1];
used_size--;
}
}
template <class T>
void Sparse_Vec<T>::clear()
{
used_size = 0;
check_small_elems_flag = false;
}
template <class T>
void Sparse_Vec<T>::clear_elem(const int i)
{
bool found = false;
int p;
it_assert0(v_size > i, "The index of the element exceeds the size of the sparse vector");
for (p=0; p<used_size; p++) {
if (index[p] == i) {
found = true;
break;
}
}
if (found) {
data[p] = data[used_size-1];
index[p] = index[used_size-1];
used_size--;
}
}
template <class T>
void Sparse_Vec<T>::set(const ivec& index_vec, const Vec<T>& v)
{
it_assert0(v_size>max(index_vec),"The indices exceeds the size of the sparse vector");
//Clear all old non-zero elements
clear();
//Add the new non-zero elements
add(index_vec,v);
}
template <class T>
void Sparse_Vec<T>::set_new(const ivec& index_vec, const Vec<T>& v)
{
int q;
int nrof_nz=v.size();
it_assert0(v_size>max(index_vec),"The indices exceeds the size of the sparse vector");
//Clear all old non-zero elements
clear();
for (q=0; q<nrof_nz; q++){
if (std::abs(v[q]) > std::abs(eps)) {
if (used_size == data_size)
resize_data(data_size*2+100);
data[used_size] = v(q);
index[used_size] = index_vec(q);
used_size++;
}
}
}
template <class T>
Sparse_Vec<T> Sparse_Vec<T>::get_subvector(int i1, int i2) const
{
it_assert0(v_size > i1 && v_size > i2 && i1<=i2 && i1>=0, "The index of the element exceeds the size of the sparse vector");
Sparse_Vec<T> r(i2-i1+1);
for (int p=0; p<used_size; p++) {
if (index[p] >= i1 && index[p] <= i2) {
if (r.used_size == r.data_size)
r.resize_data(r.data_size*2+100);
r.data[r.used_size] = data[p];
r.index[r.used_size] = index[p]-i1;
r.used_size++;
}
}
r.eps = eps;
r.check_small_elems_flag = check_small_elems_flag;
r.compact();
return r;
}
template <class T>
T Sparse_Vec<T>::sqr() const
{
T sum(0);
for (int p=0; p<used_size; p++)
sum += data[p] * data[p];
return sum;
}
template <class T>
void Sparse_Vec<T>::operator=(const Sparse_Vec<T> &v)
{
free();
v_size = v.v_size;
used_size = v.used_size;
data_size = v.data_size;
eps = v.eps;
check_small_elems_flag = v.check_small_elems_flag;
alloc();
for (int i=0; i<used_size; i++) {
data[i] = v.data[i];
index[i] = v.index[i];
}
}
template <class T>
void Sparse_Vec<T>::operator=(const Vec<T> &v)
{
free();
v_size = v.size();
used_size = 0;
data_size = std::min(v.size(), 10000);
eps = T(0);
check_small_elems_flag = false;
alloc();
for (int i=0; i<v_size; i++) {
if (v(i) != T(0)) {
if (used_size == data_size)
resize_data(data_size*2);
data[used_size] = v(i);
index[used_size] = i;
used_size++;
}
}
compact();
}
template <class T>
Sparse_Vec<T> Sparse_Vec<T>::operator-() const
{
Sparse_Vec r(v_size, used_size);
for (int p=0; p<used_size; p++) {
r.data[p] = -data[p];
r.index[p] = index[p];
}
r.used_size = used_size;
return r;
}
template <class T>
bool Sparse_Vec<T>::operator==(const Sparse_Vec<T> &v)
{
int p,q;
bool found=false;
//Remove small elements before comparing the two sparse_vectors
if (check_small_elems_flag)
remove_small_elements();
if (v_size!=v.v_size) {
//Return false if vector sizes are unequal
return false;
} else {
for(p=0;p<used_size;p++) {
for(q=0;q<v.used_size;q++) {
if (index[q] == p) {
found = true;
break;
}
}
if (found==false)
//Return false if non-zero element not found, or if elements are unequal
return false;
else if (data[p]!=v.data[q])
//Return false if non-zero element not found, or if elements are unequal
return false;
else
//Check next non-zero element
found = false;
}
}
/*Special handling if sizes do not match.
Required since v may need to do remove_small_elements() for true comparison*/
if (used_size!=v.used_size) {
if (used_size > v.used_size) {
//Return false if number of non-zero elements is less in v
return false;
}
else {
//Ensure that the remaining non-zero elements in v are smaller than v.eps
int nrof_small_elems = 0;
for(q=0;q<v.used_size;q++) {
if (std::abs(v.data[q]) <= std::abs(v.eps))
nrof_small_elems++;
}
if (v.used_size-nrof_small_elems != used_size)
//Return false if the number of "true" non-zero elements are unequal
return false;
}
}
//All elements checks => return true
return true;
}
template <class T>
void Sparse_Vec<T>::operator+=(const Sparse_Vec<T> &v)
{
int i,p;
T tmp_data;
int nrof_nz_v=v.used_size;
it_assert0(v_size == v.size(), "Attempted addition of unequal sized sparse vectors");
for (p=0; p<nrof_nz_v; p++) {
i = v.index[p];
tmp_data = v.data[p];
//get_nz(p,i,tmp_data);
add_elem(i,tmp_data);
}
check_small_elems_flag = true;
}
template <class T>
void Sparse_Vec<T>::operator+=(const Vec<T> &v)
{
int i;
it_assert0(v_size == v.size(), "Attempted addition of unequal sized sparse vectors");
for (i=0; i<v.size(); i++)
if (v(i)!=T(0))
add_elem(i,v(i));
check_small_elems_flag = true;
}
template <class T>
void Sparse_Vec<T>::operator-=(const Sparse_Vec<T> &v)
{
int i,p;
T tmp_data;
int nrof_nz_v=v.used_size;
it_assert0(v_size == v.size(), "Attempted subtraction of unequal sized sparse vectors");
for (p=0; p<nrof_nz_v; p++) {
i = v.index[p];
tmp_data = v.data[p];
//v.get_nz(p,i,tmp_data);
add_elem(i,-tmp_data);
}
check_small_elems_flag = true;
}
template <class T>
void Sparse_Vec<T>::operator-=(const Vec<T> &v)
{
int i;
it_assert0(v_size == v.size(), "Attempted subtraction of unequal sized sparse vectors");
for (i=0; i<v.size(); i++)
if (v(i)!=T(0))
add_elem(i,-v(i));
check_small_elems_flag = true;
}
template <class T>
void Sparse_Vec<T>::operator*=(const T &v)
{
int p;
for (p=0; p<used_size; p++) {
data[p]*=v;}
check_small_elems_flag = true;
}
template <class T>
void Sparse_Vec<T>::operator/=(const T &v)
{
int p;
for (p=0; p<used_size; p++) {
data[p]/=v;}
if (std::abs(eps) > 0) {
check_small_elems_flag = true;
}
}
template <class T>
T operator*(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2)
{
it_assert0(v1.v_size == v2.v_size, "Sparse_Vec<T> * Sparse_Vec<T>");
T sum(0);
Vec<T> v1f(v1.v_size);
v1.full(v1f);
for (int p=0; p<v2.used_size; p++) {
if (v1f[v2.index[p]] != T(0))
sum += v1f[v2.index[p]] * v2.data[p];
}
return sum;
}
template <class T>
T operator*(const Sparse_Vec<T> &v1, const Vec<T> &v2)
{
it_assert0(v1.size() == v2.size(), "Multiplication of unequal sized vectors attempted");
T sum(0);
for (int p1=0; p1<v1.used_size; p1++)
sum += v1.data[p1] * v2[v1.index[p1]];
return sum;
}
template <class T>
T operator*(const Vec<T> &v1, const Sparse_Vec<T> &v2)
{
it_assert0(v1.size() == v2.size(), "Multiplication of unequal sized vectors attempted");
T sum(0);
for (int p2=0; p2<v2.used_size; p2++)
sum += v1[v2.index[p2]] * v2.data[p2];
return sum;
}
template <class T>
Sparse_Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2)
{
it_assert0(v1.v_size == v2.v_size, "elem_mult(Sparse_Vec<T>, Sparse_Vec<T>)");
Sparse_Vec<T> r(v1.v_size);
ivec pos(v1.v_size);
pos = -1;
for (int p1=0; p1<v1.used_size; p1++)
pos[v1.index[p1]] = p1;
for (int p2=0; p2<v2.used_size; p2++) {
if (pos[v2.index[p2]] != -1) {
if (r.used_size == r.data_size)
r.resize_data(r.used_size*2+100);
r.data[r.used_size] = v1.data[pos[v2.index[p2]]] * v2.data[p2];
r.index[r.used_size] = v2.index[p2];
r.used_size++;
}
}
r.compact();
return r;
}
template <class T>
Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Vec<T> &v2)
{
it_assert0(v1.v_size == v2.size(), "elem_mult(Sparse_Vec<T>, Vec<T>)");
Vec<T> r(v1.v_size);
r = T(0);
for (int p1=0; p1<v1.used_size; p1++)
r[v1.index[p1]] = v1.data[p1] * v2[v1.index[p1]];
return r;
}
template <class T>
Sparse_Vec<T> elem_mult_s(const Sparse_Vec<T> &v1, const Vec<T> &v2)
{
it_assert0(v1.v_size == v2.size(), "elem_mult(Sparse_Vec<T>, Vec<T>)");
Sparse_Vec<T> r(v1.v_size);
for (int p1=0; p1<v1.used_size; p1++) {
if (v2[v1.index[p1]] != T(0)) {
if (r.used_size == r.data_size)
r.resize_data(r.used_size*2+100);
r.data[r.used_size] = v1.data[p1] * v2[v1.index[p1]];
r.index[r.used_size] = v1.index[p1];
r.used_size++;
}
}
r.compact();
return r;
}
template <class T>
Vec<T> elem_mult(const Vec<T> &v1, const Sparse_Vec<T> &v2)
{
it_assert0(v1.size() == v2.v_size, "elem_mult(Vec<T>, Sparse_Vec<T>)");
Vec<T> r(v2.v_size);
r = T(0);
for (int p2=0; p2<v2.used_size; p2++)
r[v2.index[p2]] = v1[v2.index[p2]] * v2.data[p2];
return r;
}
template <class T>
Sparse_Vec<T> elem_mult_s(const Vec<T> &v1, const Sparse_Vec<T> &v2)
{
it_assert0(v1.size() == v2.v_size, "elem_mult(Vec<T>, Sparse_Vec<T>)");
Sparse_Vec<T> r(v2.v_size);
for (int p2=0; p2<v2.used_size; p2++) {
if (v1[v2.index[p2]] != T(0)) {
if (r.used_size == r.data_size)
r.resize_data(r.used_size*2+100);
r.data[r.used_size] = v1[v2.index[p2]] * v2.data[p2];
r.index[r.used_size] = v2.index[p2];
r.used_size++;
}
}
r.compact();
return r;
}
template <class T>
Sparse_Vec<T> operator+(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2)
{
it_assert0(v1.v_size == v2.v_size, "Sparse_Vec<T> + Sparse_Vec<T>");
Sparse_Vec<T> r(v1);
ivec pos(v1.v_size);
pos = -1;
for (int p1=0; p1<v1.used_size; p1++)
pos[v1.index[p1]] = p1;
for (int p2=0; p2<v2.used_size; p2++) {
if (pos[v2.index[p2]] == -1) {// A new entry
if (r.used_size == r.data_size)
r.resize_data(r.used_size*2+100);
r.data[r.used_size] = v2.data[p2];
r.index[r.used_size] = v2.index[p2];
r.used_size++;
}
else
r.data[pos[v2.index[p2]]] += v2.data[p2];
}
r.compact();
return r;
}
template <class T>
inline Sparse_Vec<T> sparse(const Vec<T> &v)
{
Sparse_Vec<T> s(v);
return s;
}
template <class T>
inline Sparse_Vec<T> sparse(const Vec<T> &v, T epsilon)
{
Sparse_Vec<T> s(v, epsilon);
return s;
}
template <class T>
inline Vec<T> full(const Sparse_Vec<T> &s)
{
Vec<T> v;
s.full(v);
return v;
}
} //namespace itpp
#endif //__itpp_svec_h
syntax highlighted by Code2HTML, v. 0.9.1