b2vector_compressed.H

//------------------------------------------------------------------------
// b2vector_compressed.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_COMPRESSED_H__
#define __B2VECTOR_COMPRESSED_H__
 
#include "b2linear_algebra_def.H"
 
namespace b2000 { namespace b2linalg {
 
struct Vcompressed {
    typedef Vcompressed base;
    typedef Vcompressed_scale_constref const_base;
    typedef Vcompressed copy;
};
 
template <typename T>
class Vector<T, Vcompressed> {
public:
    Vector(size_t s_ = 0) : s(s_) {}
 
    Vector(const Vector& v_) : s(v_.s), index(v_.index), v(v_.v) {}
 
    template <typename T1>
    Vector(const Vector<T1, Vcompressed>& v_)
        : s(v_.s), index(v_.index), v(v_.v.begin(), v_.v.end()) {}
 
    template <typename T1>
    Vector& operator=(const Vector<T1, Vcompressed>& v_) {
        s = v_.s;
        index = v_.index;
        v.assign(v_.v.begin(), v_.v.end());
        return *this;
    }
 
    Vector& operator=(const Vector<T, Vdense>& v_) {
        s = v_.size();
        size_t nb_nz = 0;
        for (typename Vector<T, Vdense>::const_iterator i = v_.begin(); i != v_.end(); ++i) {
            if (*i != T(0)) { ++nb_nz; }
        }
        index.resize(nb_nz);
        v.resize(nb_nz);
        nb_nz = 0;
        for (typename Vector<T, Vdense>::const_iterator i = v_.begin(); i != v_.end(); ++i) {
            if (*i != T(0)) {
                index[nb_nz] = i - v_.begin();
                v[nb_nz] = *i;
                ++nb_nz;
            }
        }
        return *this;
    }
 
    size_t size() const { return s; }
 
    size_t get_nb_nonzero() const { return index.size(); }
 
    bool is_null() const { return this == &null; }
 
    void set_zero() {
        index.clear();
        v.clear();
    }
 
    void resize(size_t s_, size_t snn) {
        s = s_;
        index.resize(snn);
        v.resize(snn);
    }
 
    void get_value(size_t*& index_, T*& value_) {
        index_ = &index[0];
        value_ = &v[0];
    }
 
    auto& get_index() const { return index; }
 
    Vector& operator=(
          const Vector<T, MVProd<Mcompressed_col_st_constref, Vcompressed_scale_constref> >& v_) {
        if (v_.m.trans) { UnimplementedError() << THROW; }
 
        std::vector<size_t> vtmp_index;
        std::vector<T> vtmp_value;
        if (v_.v.v == &v[0]) {
            vtmp_index.swap(index);
            vtmp_value.swap(v);
        } else {
            resize(v_.size(), 0);
        }
 
        size_t* tmp_index = TemporaryBuffer<size_t>::get(v_.m.size1());
        const size_t end_list_flag = v_.m.size1();
        T* tmp_value = TemporaryBuffer<T>::get(v_.m.size1());
        size_t end_list = end_list_flag;
 
        T scale = v_.scale * v_.m.scale * v_.v.scale;
 
        const size_t* b_rowind_begin = v_.v.index;
        const size_t* const b_rowind_end = b_rowind_begin + v_.v.snn;
        const T* b_value_begin = v_.v.v;
        while (b_rowind_begin != b_rowind_end) {
            const size_t* a_rowind_begin = v_.m.index + v_.m.si[*b_rowind_begin];
            const size_t* const a_rowind_end = v_.m.index + v_.m.si[*b_rowind_begin + 1];
            const T* a_value_begin = v_.m.m + v_.m.si[*b_rowind_begin++];
            const T b_value_begin_v = *b_value_begin++;
            while (a_rowind_begin != a_rowind_end) {
                tmp_value[*a_rowind_begin] += *a_value_begin++ * b_value_begin_v;
                if (!(std::numeric_limits<size_t>::max() - tmp_index[*a_rowind_begin])) {
                    tmp_index[*a_rowind_begin] = end_list;
                    end_list = *a_rowind_begin;
                }
                ++a_rowind_begin;
            }
        }
        std::vector<std::pair<size_t, T> > res_tmp;
        while (end_list != end_list_flag) {
            res_tmp.push_back(std::pair<size_t, T>(end_list, tmp_value[end_list]));
            const size_t end_list_next = tmp_index[end_list];
            tmp_index[end_list] = std::numeric_limits<size_t>::max();
            tmp_value[end_list] = T(0);
            end_list = end_list_next;
        }
        std::sort(res_tmp.begin(), res_tmp.end(), CompareFirstOfPair());
 
        v.reserve(res_tmp.size());
        for (typename std::vector<std::pair<size_t, T> >::iterator i = res_tmp.begin();
             i != res_tmp.end(); ++i) {
            index.push_back(i->first);
            v.push_back(i->second * scale);
        }
 
        return *this;
    }
 
    Vector& operator=(
          const Vector<
                T, MVProd<Mstructured_constref<Mpacked_st_constref>, Vcompressed_scale_constref> >&
                v_) {
        v.clear();
        const size_t* rowind1 = v_.m.index;
        size_t i1 = 0;
        const size_t i1_end = v_.m.m.s;
        size_t i2 = 0;
        const size_t i2_end = v_.v.snn;
        for (;;) {
            if (rowind1[i1] < rowind1[i2]) {
                if (++i1 == i1_end) { break; }
            } else if (rowind1[i1] > rowind1[i2]) {
                if (++i2 == i2_end) { break; }
            } else {
                v.resize(v_.m.m.s);
                const T v2 = v_.v[i2];
                for (size_t i = 0; i != v.size(); ++i) { v[i] += v_.m.m(i, rowind1[i1]) * v2; }
                if (++i1 == i1_end) { break; }
                if (++i2 == i2_end) { break; }
            }
        }
        if (!v.empty()) {
            s = v_.m.s;
            index.assign(v_.m.index, v_.m.index + v.size());
        }
        return *this;
    }
 
    Vector& operator=(
          const Vector<
                T, MVProd<Mstructured_constref<Mpacked_st_constref>, Vincrement_scale_constref> >&
                v_) {
        s = v_.m.s;
        v.resize(v_.m.m.s);
        index.assign(v_.m.index, v_.m.index + v.size());
        Vector<T> vb(v_.m.m.s);
        for (size_t i = 0; i != vb.size(); ++i) { vb[i] = v_.v[index[i]]; }
        Vector<T, Vincrement_ref> va(index.size(), &v[0]);
        va = v_.scale * (v_.m.m * vb);
        return *this;
    }
 
    T operator[](size_t i) const {
        const std::vector<size_t>::const_iterator ii = std::find(index.begin(), index.end(), i);
        if (ii == index.end()) { return 0; }
        return v[ii - index.begin()];
    }
 
    static Vector null;
 
    friend logging::Logger& operator<<(logging::Logger& l, const Vector& v) {
        l << "compressed vector ";
        l.write(v.s, &v.index[0], 1, "index");
        l.write(v.s, &v.v[0], 1, "value");
        return l;
    }
 
private:
    MVFRIEND;
    size_t s;
    std::vector<size_t> index;
    std::vector<T> v;
};
 
template <typename T>
Vector<T, Vcompressed> Vector<T, Vcompressed>::null;
 
struct Vcompressed_scale_constref {
    typedef Vcompressed_scale_constref base;
    typedef Vcompressed_scale_constref const_base;
    typedef Vcompressed copy;
};
 
template <typename T>
class Vector<T, Vcompressed_scale_constref> {
public:
    Vector() : s(0), snn(0), scale(1), index(0), v(0) {}
 
    Vector(const Vector& v_) : s(v_.s), snn(v_.snn), scale(v_.scale), index(v_.index), v(v_.v) {}
 
    Vector(const Vector<T, Vcompressed>& v_)
        : s(v_.s), snn(v_.index.size()), scale(1), index(&v_.index[0]), v(&v_.v[0]) {}
 
    Vector(size_t s_, size_t snn_, const size_t* index_, const T* value_, T scale_)
        : s(s_), snn(snn_), scale(scale_), index(index_), v(value_) {}
 
    size_t size() const { return s; }
 
    size_t get_nb_nonzero() const { return snn; }
 
    bool is_null() const { return this == &null; }
 
    T operator[](size_t i) const {
        const size_t* ii = std::find(index, index + snn, i);
        if (ii == index + snn) { return 0; }
        return v[ii - index];
    }
 
    Vector& operator*=(T scale_) {
        scale *= scale_;
        return *this;
    }
 
    Vector operator*(T scale_) { return Vector(s, snn, index, v, scale * scale_); }
 
    static Vector null;
 
private:
    MVFRIEND;
    size_t s;
    size_t snn;
    T scale;
    const size_t* index;
    const T* v;
};
 
template <typename T>
Vector<T, Vcompressed_scale_constref> Vector<T, Vcompressed_scale_constref>::null;
 
}}  // namespace b2000::b2linalg
 
#endif