b2vector_index.H

//------------------------------------------------------------------------
// b2vector_index.H --
//
//     A framework for generic linear algebraic (matrix) computations.
//
// written by Mathias Doreille
//
// Copyright (c) 2004-2012 SMR Engineering & Development SA
//                         2502 Bienne, Switzerland
//
// All Rights Reserved.  Proprietary source code.  The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR.
//------------------------------------------------------------------------
 
#ifndef __B2VECTOR_INDEX_H__
#define __B2VECTOR_INDEX_H__
 
#include "b2linear_algebra_def.H"
 
namespace b2000 { namespace b2linalg {
 
struct Vindex_constref {
    typedef Vindex_constref base;
    typedef Vindex_scale_constref const_base;
    typedef Vdense copy;
};
 
template <typename T>
class Vector<T, Vindex_constref> {
public:
    Vector() : s(0), v(0), index(0) {}
 
    Vector(size_t s_, const T* v_, const size_t* index_) : s(s_), v(v_), index(index_) {}
 
    Vector(const Vector& v_) : s(v_.s), v(v_.v), index(v_.index) {}
 
    Vector(const Vector<T, Vindex_ref>& v_) : s(v_.s), v(v_.v), index(v_.index) {}
 
    bool is_null() const { return v == 0; }
 
    size_t size() const { return s; }
 
    const T& operator[](size_t i) const { return v[index[i]]; }
 
    Vector operator()(const Interval& i) const { return Vector(i.size(), v, index + i[0]); }
 
    static Vector null;
 
protected:
    size_t s;
    const T* v;
    const size_t* index;
    MVFRIEND;
};
 
template <typename T>
Vector<T, Vindex_constref> Vector<T, Vindex_constref>::null;
 
struct Vindex1_constref {};
 
template <typename T>
class Vector<T, Vindex1_constref> : public Vector<T, Vindex_constref> {
public:
    Vector(T default_value_ = T(0))
        : Vector<T, Vindex_constref>(), s_v(0), default_value(default_value_) {}
 
    Vector(size_t s_, const T* v_, size_t s_v_, const size_t* index_)
        : Vector<T, Vindex_constref>(s_, v_, index_), s_v(s_v_), default_value(T(0)) {}
 
    Vector(const Vector& v_)
        : Vector<T, Vindex_constref>(v_), s_v(v_.s_v), default_value(v_.default_value) {}
 
    void set_default_value(T default_value_) { default_value = default_value_; }
 
    const T& operator[](size_t i) const {
        if (s_v == 0 || Vector<T, Vindex_constref>::index[i] < s_v) {
            return Vector<T, Vindex_constref>::v[Vector<T, Vindex_constref>::index[i]];
        }
        return default_value;
    }
 
private:
    size_t s_v;
    T default_value;
};
 
struct Vindex_ref {
    typedef Vindex_ref base;
    typedef Vindex_scale_constref const_base;
    typedef Vdense copy;
};
 
template <typename T>
class Vector<T, Vindex_ref> {
public:
    Vector() : s(0), v(0), index(0) {}
 
    Vector(size_t s_ = 0, T* v_ = 0, const size_t* index_ = 0) : s(s_), v(v_), index(index_) {}
 
    Vector(const Vector& v_) : s(v_.s), v(v_.v), index(v_.index) {}
 
    bool is_null() const { return v == 0; }
 
    size_t size() const { return s; }
 
    void set_zero() {
        for (size_t i = 0; i != s; ++i) { v[i] = T(0); }
    }
 
    const T& operator[](size_t i) const { return v[index[i]]; }
 
    T& operator[](size_t i) { return v[index[i]]; }
 
    Vector operator()(const Interval& i) const { return Vector(i.size(), v, index + i[0]); }
 
    Vector& operator=(const Vector& v_) {
        if (s != v_.s) { Exception() << "The two vectors do not have the same size." << THROW; }
        for (size_t i = 0; i != s; ++i) { v[index[i]] = v_[v_.index[i]]; }
        return *this;
    }
 
    Vector& operator=(const Vector<T, Vdense>& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size." << THROW;
        }
        for (size_t i = 0; i != s; ++i) { v[index[i]] = v_[i]; }
        return *this;
    }
 
    Vector& operator=(const Vector<T, Vindex_scale_constref>& v_) {
        if (s != v_.s) { Exception() << "The two vectors do not have the same size." << THROW; }
        if (v_.v == v && v_.index == index && v_.scale == T(1)) { return *this; }
        for (size_t i = 0; i != s; ++i) { v[index[i]] = v_[i]; }
        return *this;
    }
 
    Vector& operator+=(const Vector<T, Vindex_scale_constref>& v_) {
        if (s != v_.s) { Exception() << "The two vectors do not have the same size." << THROW; }
        for (size_t i = 0; i != s; ++i) { v[index[i]] += v_[i]; }
        return *this;
    }
 
    Vector& operator=(const Vector<T, Vincrement_scale_constref>& v_) {
        if (s != v_.s) { Exception() << "The two vectors do not have the same size." << THROW; }
        for (size_t i = 0; i != s; ++i) { v[index[i]] = v_[i]; }
        return *this;
    }
 
    Vector& operator+=(const Vector<T, Vincrement_scale_constref>& v_) {
        if (s != v_.s) { Exception() << "The two vectors do not have the same size." << THROW; }
        for (size_t i = 0; i != s; ++i) { v[index[i]] += v_[i]; }
        return *this;
    }
 
    Vector& operator+=(const Vector<T, Sum<Vindex_scale_constref, Vindex_scale_constref> >& v_) {
        if (s != v_.s) { Exception() << "The two vectors do not have the same size." << THROW; }
        for (size_t i = 0; i != s; ++i) { v[index[i]] += v_[i]; }
        return *this;
    }
 
    Vector& operator=(
          const Vector<T, MVProd<Mrectangle_increment_st_constref, Vindex_scale_constref> >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        Vector<T, Vdense> tmp1(s);
        Vector<T, Vdense> tmp2(v_.v);
        tmp1 = v_.scale * (v_.m * tmp2);
        *this = tmp1;
        return *this;
    }
 
    Vector& operator+=(
          const Vector<T, MVProd<Mrectangle_increment_st_constref, Vindex_scale_constref> >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        Vector<T, Vdense> tmp1(s);
        Vector<T, Vdense> tmp2(v_.v);
        tmp1 = v_.scale * (v_.m * tmp2);
        *this += tmp1;
        return *this;
    }
 
    Vector& operator=(
          const Vector<
                T, Sum<Vindex_scale_constref,
                       MVProd<Mrectangle_increment_st_constref, Vindex_scale_constref> > >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        Vector<T, Vdense> tmp1(s);
        Vector<T, Vdense> tmp2(v_.v2.v);
        tmp1 = (v_.scale * v_.v2.scale) * (v_.v2.m * tmp2);
        *this = v_.scale * v_.v1;
        *this += tmp1;
        return *this;
    }
 
    Vector& operator=(const Vector<T, MVProd<Mpacked_st_constref, Vindex_scale_constref> >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        Vector<T, Vdense> tmp1(s);
        Vector<T, Vdense> tmp2(v_.v);
        tmp1 = v_.scale * (v_.m * tmp2);
        *this = tmp1;
        return *this;
    }
 
    Vector& operator+=(const Vector<T, MVProd<Mpacked_st_constref, Vindex_scale_constref> >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        Vector<T, Vdense> tmp1(s);
        Vector<T, Vdense> tmp2(v_.v);
        tmp1 = v_.scale * (v_.m * tmp2);
        *this += tmp1;
        return *this;
    }
 
    Vector& operator=(const Vector<
                      T, Sum<Vindex_scale_constref,
                             MVProd<Mpacked_st_constref, Vindex_scale_constref> > >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        Vector<T, Vdense> tmp1(s);
        Vector<T, Vdense> tmp2(v_.v2.v);
        tmp1 = (v_.scale * v_.v2.scale) * (v_.v2.m * tmp2);
        *this = v_.scale * v_.v1;
        *this += tmp1;
        return *this;
    }
 
    Vector& operator=(
          const Vector<T, MVProd<Mcompressed_col_st_constref, Vindex_scale_constref> >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
 
        const size_t* colind_begin = v_.m.si;
        const size_t* rowind_begin = v_.m.index + *colind_begin;
        const T* value_begin = v_.m.m + *colind_begin;
        T scale = v_.scale * v_.v.scale * v_.m.scale;
        if (v_.m.trans) {
            for (size_t i = 0; i != s; ++i) {
                const size_t* const rowind_end = v_.m.index + *++colind_begin;
                T val = T(0);
                while (rowind_begin != rowind_end) {
                    val += *value_begin++ * v_.v[*rowind_begin++];
                }
                (*this)[i] = scale * val;
            }
        } else {
            set_zero();
            for (size_t i = 0; i != v_.v.s; ++i) {
                const size_t* const rowind_end = v_.m.index + *++colind_begin;
                const T val = v_.v[i];
                while (rowind_begin != rowind_end) {
                    (*this)[*rowind_begin++] += scale * *value_begin++ * val;
                }
            }
        }
        return *this;
    }
 
    Vector& operator=(const Vector<
                      T, Sum<Vindex_scale_constref,
                             MVProd<Mcompressed_col_st_constref, Vindex_scale_constref> > >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
 
        (*this) = v_.scale * v_.v1;
        const size_t* colind_begin = v_.v2.m.si;
        const size_t* rowind_begin = v_.v2.m.index + *colind_begin;
        const T* value_begin = v_.v2.m.m + *colind_begin;
        T scale = v_.scale * v_.v.scale * v_.m.scale;
        if (v_.m.trans) {
            for (size_t i = 0; i != s; ++i) {
                const size_t* const rowind_end = v_.v2.m.index + *++colind_begin;
                T val = T(0);
                while (rowind_begin != rowind_end) {
                    val += *value_begin++ * v_.v2.v[*rowind_begin++];
                }
                (*this)[i] += scale * val;
            }
        } else {
            for (size_t i = 0; i != v_.v2.v.s; ++i) {
                const size_t* const rowind_end = v_.v2.m.index + *++colind_begin;
                const T val = v_.v2[i];
                while (rowind_begin != rowind_end) {
                    (*this)[*rowind_begin++] += scale * *value_begin++ * val;
                }
            }
        }
        return *this;
    }
 
    static Vector null;
 
private:
    size_t s;
    T* v;
    const size_t* index;
    MVFRIEND;
};
 
template <typename T>
Vector<T, Vindex_ref> Vector<T, Vindex_ref>::null;
 
struct Vindex_scale_constref {
    typedef Vindex_scale_constref base;
    typedef Vindex_scale_constref const_base;
    typedef Vdense copy;
};
 
template <typename T>
class Vector<T, Vindex_scale_constref> {
public:
    Vector() : s(0), v(0), index(0), scale(0) {}
 
    Vector(size_t s_, const T* v_, const size_t* index_, T scale_ = T(1))
        : s(s_), v(v_), index(index_), scale(scale_) {}
 
    Vector(const Vector& v_) : s(v_.s), v(v_.v), index(v_.index), scale(v_.scale) {}
 
    Vector(const Vector<T, Vindex_ref>& v_) : s(v_.s), v(v_.v), index(v_.index), scale(1) {}
 
    Vector(const Vector<T, Vindex_constref>& v_) : s(v_.s), v(v_.v), index(v_.index), scale(1) {}
 
    bool is_null() const { return v == 0; }
 
    size_t size() const { return s; }
 
    const T operator[](size_t i) const { return scale * v[index[i]]; }
 
    Vector operator()(const Interval& i) const { return Vector(i.size(), v, index + i[0], scale); }
 
    Vector& operator*(T scale_) {
        scale *= scale_;
        return *this;
    }
 
    T operator*(const Vector<T, Vincrement_scale_constref>& v_) const {
        if (s != v_.s) { Exception() << "The two vectors do not have the same size." << THROW; }
        T tmp = 0;
        for (size_t i = 0; i != s; ++i) { tmp += v[index[i]] * v_.v[v_.index[i]]; }
        return tmp * scale * v_.scale;
    }
 
    static Vector null;
 
private:
    size_t s;
    const T* v;
    const size_t* index;
    T scale;
    MVFRIEND;
};
 
template <typename T>
Vector<T, Vindex_scale_constref> Vector<T, Vindex_scale_constref>::null;
 
}}  // namespace b2000::b2linalg
 
#endif