b2vector_increment.H

//------------------------------------------------------------------------
// b2vector_increment.H --
//
//     A framework for generic linear algebraic (matrix) computations.
//
// written by Mathias Doreille
//
// Copyright (c) 2004-2012,2017
//     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_INCREMENT_H__
#define __B2VECTOR_INCREMENT_H__
 
#include "b2linear_algebra_def.H"
 
namespace b2000 { namespace b2linalg {
 
struct Vincrement_constref {
    typedef Vincrement_scale_constref base;
    typedef Vincrement_scale_constref const_base;
    typedef Vdense copy;
};
 
template <typename T>
class Vector<T, Vincrement_constref> {
public:
    Vector() : s(0), v(0), incr(0) {}
 
    Vector(size_t s_, const T* v_, size_t incr_ = 1) : s(s_), v(v_), incr(incr_) {}
 
    Vector(const Vector& v_) : s(v_.s), v(v_.v), incr(v_.incr) {}
 
    Vector(const Vector<T, Vdense>& v_) : s(v_.size()), v(&v_[0]), incr(1) {}
 
    Vector(const Vector<T, Vdense_constref>& v_) : s(v_.s), v(v_.v), incr(1) {}
 
    Vector(const std::vector<T>& v_) : s(v_.size()), v(&v_[0]), incr(1) {}
 
    bool is_null() const { return v == 0; }
 
    size_t size() const { return s; }
 
    const T& operator[](size_t i) const { return v[i * incr]; }
 
    Vector<T, Vincrement_constref> operator()(const Interval& i) const {
        return Vector<T, Vincrement_constref>(i.size(), &v[i[0] * incr], incr);
    }
 
    Vector<T, Vincrement_constref> operator()(const Slice& i) const {
        return Vector<T, Vincrement_constref>(i.size(), &v[i[0] * incr], i.increment() * incr);
    }
 
    T norm2() const { return blas::nrm2(s, v, incr); }
 
    static Vector null;
 
    friend logging::Logger& operator<<(logging::Logger& l, const Vector& v) {
        l << "vector ";
        l.write(v.s, v.v, v.incr);
        return l;
    }
 
private:
    size_t s;
    const T* v;
    size_t incr;
    MVFRIEND;
};
 
template <typename T>
Vector<T, Vincrement_constref> Vector<T, Vincrement_constref>::null;
 
struct Vincrement_ref {
    typedef Vincrement_ref base;
    typedef Vincrement_scale_constref const_base;
    typedef Vdense copy;
};
 
template <typename T>
class Vector<T, Vincrement_ref> {
public:
    Vector() : s(0), v(0), incr(0) {}
 
    Vector(size_t s_, T* v_, size_t incr_ = 1) : s(s_), v(v_), incr(incr_) {}
 
    Vector(const Vector& v_) : s(v_.s), v(v_.v), incr(v_.incr) {}
 
    Vector(Vector<T, Vdense>& v_) : s(v_.size()), v(&v_[0]), incr(1) {}
 
    Vector(Vector<T, Vdense_ref>& v_) : s(v_.size()), v(v_.v), incr(1) {}
 
    bool is_null() const { return v == 0; }
 
    size_t size() const { return s; }
 
    void set_zero() {
        T const* vv_end = v + s * incr;
        for (T* vv = v; vv != vv_end; vv += incr) { *vv = T(0); }
    }
 
    const T& operator[](size_t i) const { return v[i * incr]; }
 
    T& operator[](size_t i) { return v[i * incr]; }
 
    Vector<T, Vincrement_ref> operator()(const Interval& i) {
        return Vector<T, Vincrement_ref>(i.size(), &v[i[0] * incr], incr);
    }
 
    Vector<T, Vincrement_ref> operator()(const Slice& i) {
        return Vector<T, Vincrement_ref>(i.size(), &v[i[0] * incr], i.increment() * incr);
    }
 
    template <typename STORAGE>
    Vector& operator+=(const Vector<T, STORAGE>& v_) {
        *this = *this + v_;
        return *this;
    }
 
    template <typename STORAGE>
    Vector& operator-=(const Vector<T, STORAGE>& v_) {
        *this = *this - v_;
        return *this;
    }
 
    Vector& operator=(const Vector& v_) {
        if (s != v_.s) {
            Exception() << "The two vectors do not have the same size, " << s << " and " << v_.s
                        << "." << THROW;
        }
        if (v_.v != v || v_.incr != incr) { blas::copy(s, v_.v, v_.incr, v, incr); }
        return *this;
    }
 
    template <typename T1>
    Vector& operator=(const Vector<T1, Vdense>& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        if (incr == 1) {
            std::copy(&v_[0], &v_[0] + s, v);
        } else {
            T* out = v;
            T* out_end = out + s * incr;
            const T1* in = &v_[0];
            for (; out < out_end; ++in, out += incr) { *out = *in; }
        }
        return *this;
    }
 
    template <typename T1>
    Vector& operator=(const Vector<T1, Vincrement_ref>& v_) {
        if (s != v_.s) {
            Exception() << "The two vectors do not have the same size, " << s << " and " << v_.s
                        << "." << THROW;
        }
        if (v_.incr == 1 && incr == 1) {
            std::copy(v_.v, v_.v + s, v);
        } else {
            T* out = v;
            T* out_end = out + s * incr;
            const T1* in = v_.v;
            for (; out < out_end; in += v_.incr, out += incr) { *out = *in; }
        }
        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, " << s << " and " << v_.s
                        << "." << THROW;
        }
        if (v_.v != v || v_.incr != incr) { blas::copy(s, v_.v, v_.incr, v, incr); }
        if (v_.scale != T(1)) { blas::scal(s, v_.scale, v, incr); }
        return *this;
    }
 
    template <typename T1>
    Vector& operator=(const Vector<T1, Vincrement_scale_constref>& v_) {
        if (s != v_.s) {
            Exception() << "The two vectors do not have the same size, " << s << " and " << v_.s
                        << "." << THROW;
        }
        if (v_.incr == 1 && incr == 1) {
            std::copy(v_.v, v_.v + s, v);
        } else {
            T* out = v;
            T* out_end = out + s * incr;
            const T1* in = v_.v;
            for (; out < out_end; in += v_.incr, out += incr) { *out = *in; }
        }
        if (v_.scale != T(1)) { blas::scal(s, v_.scale, v, incr); }
        return *this;
    }
 
    Vector& operator=(
          const Vector<T, Sum<Vincrement_scale_constref, Vincrement_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;
        blas::axpy(s, v_.v2.scale * v_.scale, v_.v2.v, v_.v2.incr, v, incr);
        return *this;
    }
 
    Vector& operator=(
          const Vector<T, Sum<Vincrement_scale_constref, Vcompressed_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* index = v_.v2.index;
        const size_t* index_end = index + v_.v2.snn;
        const T* value = v_.v2.v;
        T scale = v_.scale * v_.v2.scale;
        for (; index != index_end; ++index, ++value) { v[*index] += scale * *value; }
        return *this;
    }
 
    Vector& operator=(
          const Vector<T, MVProd<Mrectangle_increment_st_constref, Vincrement_scale_constref> >&
                v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        blas::gemv(
              v_.m.trans ? 't' : 'n', v_.m.trans ? v_.m.size2() : v_.m.size1(),
              v_.m.trans ? v_.m.size1() : v_.m.size2(), v_.scale * v_.v.scale * v_.m.scale, v_.m.m,
              v_.m.ldm, v_.v.v, v_.v.incr, 0, v, incr);
        return *this;
    }
 
    Vector& operator=(
          const Vector<
                T, Sum<Vincrement_scale_constref,
                       MVProd<Mrectangle_increment_st_constref, Vincrement_scale_constref> > >&
                v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        if (v_.v1.v != v || v_.v1.incr != incr) { blas::copy(s, v_.v1.v, v_.v1.incr, v, incr); }
        blas::gemv(
              v_.v2.m.trans ? 't' : 'n', v_.v2.m.trans ? v_.v2.m.size2() : v_.v2.m.size1(),
              v_.v2.m.trans ? v_.v2.m.size1() : v_.v2.m.size2(),
              v_.scale * v_.v2.v.scale * v_.v2.m.scale, v_.v2.m.m, v_.v2.m.ldm, v_.v2.v.v,
              v_.v2.v.incr, v_.scale * v_.v1.scale, v, incr);
        return *this;
    }
 
    Vector& operator=(
          const Vector<
                T, Sum<Vincrement_scale_constref,
                       MVProd<Mrectangle_increment_st_constref, Vcompressed_scale_constref> > >&
                v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
 
        const T scale = v_.v2.scale * v_.v2.v.scale * v_.v2.m.scale;
        if (v_.v2.m.trans) {
            const T* m = v_.v2.m.m;
            for (size_t i = 0; i != s; ++i) {
                T tmp = 0;
                for (size_t j = 0; j != v_.v2.v.snn; ++j) {
                    tmp += m[v_.v2.v.index[j]] * v_.v2.v.v[j];
                }
                v[i] = v_.scale * (v_.v1[i] + scale * tmp);
                m += v_.v2.m.ldm;
            }
        } else {
            UnimplementedError() << THROW;
        }
        return *this;
    }
 
    Vector& operator=(
          const Vector<T, MVProd<Mpacked_st_constref, Vincrement_scale_constref> >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        blas::spmv(
              v_.m.upper ? 'u' : 'l', v_.m.size1(), v_.scale * v_.v.scale * v_.m.scale, v_.m.m,
              v_.v.v, v_.v.incr, 0, v, incr);
        return *this;
    }
 
    Vector& operator=(const Vector<
                      T, Sum<Vincrement_scale_constref,
                             MVProd<Mpacked_st_constref, Vincrement_scale_constref> > >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        if (v_.v1.v != v || v_.v1.incr != incr) { blas::copy(s, v_.v1.v, v_.v1.incr, v, incr); }
        blas::spmv(
              v_.v2.m.upper ? 'u' : 'l', v_.v2.m.size1(), v_.scale * v_.v2.v.scale * v_.v2.m.scale,
              v_.v2.m.m, v_.v2.v.v, v_.v2.v.incr, v_.scale * v_.v1.scale, v, incr);
        return *this;
    }
 
    Vector& operator=(
          const Vector<T, MVProd<Mcompressed_col_st_constref, Vincrement_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) {
            T* result_begin = v;
            const T* const result_end = v + s * incr;
            while (result_begin != result_end) {
                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++];
                }
                *result_begin = scale * val;
                result_begin += incr;
            }
        } else {
            set_zero();
            const T* input_begin = v_.v.v;
            const T* const input_end = input_begin + v_.v.s * v_.v.incr;
            while (input_begin != input_end) {
                const size_t* const rowind_end = v_.m.index + *++colind_begin;
                const T val = scale * *input_begin;
                input_begin += v_.v.incr;
                while (rowind_begin != rowind_end) {
                    (*this)[*rowind_begin++] += *value_begin++ * val;
                }
            }
        }
        return *this;
    }
 
    Vector& operator=(
          const Vector<
                T, Sum<Vincrement_scale_constref,
                       MVProd<Mcompressed_col_st_constref, Vcompressed_scale_constref> > >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
 
        T scale = v_.scale * v_.v2.scale * v_.v2.v.scale * v_.v2.m.scale;
        if (v_.v2.m.trans) {
            UnimplementedError() << THROW;
        } else {
            *this = v_.scale * v_.v1;
            for (size_t i = 0; i != v_.v2.v.snn; ++i) {
                size_t j = v_.v2.m.si[v_.v2.v.index[i]];
                const size_t j_end = v_.v2.m.si[v_.v2.v.index[i] + 1];
                const T val = scale * v_.v2.v.v[i];
                for (; j != j_end; ++j) { (*this)[v_.v2.m.index[j]] += v_.v2.m.m[j] * val; }
            }
        }
        return *this;
    }
 
    Vector& operator=(const Vector<
                      T, MVProd<
                               MMProd<Mcompressed_col_st_constref, Mcompressed_col_st_constref>,
                               Vincrement_scale_constref> >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        if (v_.m.m1.trans || !v_.m.m2.trans) {
            Vector<T, Vdense> tmp;
            tmp = v_.m.m2 * v_.v;
            *this = (v_.scale * v_.m.scale) * (v_.m.m1 * tmp);
            return *this;
        }
 
        const size_t* rowind_begin_1 = v_.m.m1.index;
        const T* value_begin_1 = v_.m.m1.m;
        const size_t* colind_begin_1 = v_.m.m1.si;
        const size_t* rowind_begin_2 = v_.m.m2.index;
        const T* value_begin_2 = v_.m.m2.m;
        const size_t* colind_begin_2 = v_.m.m2.si;
        T scale = v_.scale * v_.v.scale * v_.m.scale * v_.m.m1.scale * v_.m.m2.scale;
        set_zero();
        for (size_t i = 0; i != size(); ++i) {
            T val = T(0);
            const size_t* rowind_end = v_.m.m2.index + *++colind_begin_2;
            while (rowind_begin_2 != rowind_end) {
                val += *value_begin_2++ * v_.v[*rowind_begin_2++];
            }
            val *= scale;
            rowind_end = v_.m.m1.index + *++colind_begin_1;
            while (rowind_begin_1 != rowind_end) {
                (*this)[*rowind_begin_1++] += *value_begin_1++ * val;
            }
        }
        return *this;
    }
 
    Vector& operator=(const Vector<
                      T, MVProd<
                               Msub_constref<Mcompressed_col_constref, Index, Index>,
                               Vincrement_scale_constref> >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        if (!v_.m.index1.is_all()) { UnimplementedError() << THROW; }
 
        const T* value_begin = v_.m.m.m;
        const size_t* colind_begin = &v_.m.index2[0];
        T scale = v_.scale * v_.v.scale;
        set_zero();
        const T* input_begin = v_.v.v;
        const T* const input_end = input_begin + v_.v.s * v_.v.incr;
        while (input_begin != input_end) {
            const size_t* rowind_begin = v_.m.m.index + v_.m.m.si[*colind_begin];
            const size_t* const rowind_end = v_.m.m.index + v_.m.m.si[*colind_begin++ + 1];
            const T val = scale * *input_begin;
            input_begin += v_.v.incr;
            while (rowind_begin != rowind_end) { (*this)[*rowind_begin++] += *value_begin++ * val; }
        }
        return *this;
    }
 
    Vector& operator=(
          const Vector<
                T, Sum<Vincrement_scale_constref,
                       MVProd<Mcompressed_col_st_constref, Vincrement_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_.v2.scale * v_.v2.v.scale * v_.v2.m.scale;
        if (v_.v2.m.trans) {
            T* result_begin = v;
            const T* const result_end = v + s * incr;
            while (result_begin != result_end) {
                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++];
                }
                *result_begin += scale * val;
                result_begin += incr;
            }
        } else {
            const T* input_begin = v_.v2.v.v;
            const T* const input_end = input_begin + v_.v2.v.s * v_.v2.v.incr;
            while (input_begin != input_end) {
                const size_t* const rowind_end = v_.v2.m.index + *++colind_begin;
                const T val = scale * *input_begin;
                input_begin += v_.v2.v.incr;
                while (rowind_begin != rowind_end) {
                    (*this)[*rowind_begin++] += *value_begin++ * val;
                }
            }
        }
        return *this;
    }
 
    Vector& operator=(
          const Vector<
                T, Sum<Vincrement_scale_constref,
                       MVProd<Msym_compressed_col_st_constref, Vincrement_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;
 
        T scale = v_.scale * v_.v2.scale * v_.v2.v.scale * v_.v2.m.scale;
        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* result_begin = v;
            const T* const result_end = v + s * incr;
            while (result_begin != result_end) {
                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++];
                }
                *result_begin += scale * val;
                result_begin += incr;
            }
        }
        colind_begin = v_.v2.m.si;
        rowind_begin = v_.v2.m.index + *colind_begin;
        value_begin = v_.v2.m.m + *colind_begin;
        {
            const T* input_begin = v_.v2.v.v;
            const T* const input_end = input_begin + v_.v2.v.s * v_.v2.v.incr;
            size_t i = 0;
            while (input_begin != input_end) {
                const size_t* const rowind_end = v_.v2.m.index + *++colind_begin;
                const T val = scale * *input_begin;
                input_begin += v_.v2.v.incr;
                if (rowind_begin < rowind_end && *rowind_begin == i) {
                    ++rowind_begin;
                    ++value_begin;
                }
                ++i;
                while (rowind_begin < rowind_end) {
                    (*this)[*rowind_begin++] += *value_begin++ * val;
                }
            }
        }
 
        return *this;
    }
 
    Vector& operator=(
          const Vector<T, MVProd<Msym_compressed_col_st_constref, Vincrement_scale_constref> >&
                v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
 
        T scale = v_.scale * v_.v.scale * v_.m.scale;
        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* result_begin = v;
            const T* const result_end = v + s * incr;
            while (result_begin != result_end) {
                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++];
                }
                *result_begin = scale * val;
                result_begin += incr;
            }
        }
        colind_begin = v_.m.si;
        rowind_begin = v_.m.index + *colind_begin;
        value_begin = v_.m.m + *colind_begin;
        {
            const T* input_begin = v_.v.v;
            const T* const input_end = input_begin + v_.v.s * v_.v.incr;
            size_t i = 0;
            while (input_begin != input_end) {
                const size_t* const rowind_end = v_.m.index + *++colind_begin;
                const T val = scale * *input_begin;
                input_begin += v_.v.incr;
                if (rowind_begin < rowind_end && *rowind_begin == i) {
                    ++rowind_begin;
                    ++value_begin;
                }
                ++i;
                while (rowind_begin < rowind_end) {
                    (*this)[*rowind_begin++] += *value_begin++ * val;
                }
            }
        }
        return *this;
    }
 
    Vector& operator=(
          const Vector<
                T, MVProd<
                         MMProd<Msym_compressed_col_st_constref, Msym_compressed_col_st_constref>,
                         Vincrement_scale_constref> >& v_) {
        UnimplementedError() << THROW;
        return *this;
    }
 
    Vector& operator=(const Vector<
                      T, MVProd<
                               Minverse_constref<Msym_compressed_col_update_inv>,
                               Vincrement_scale_constref> >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        if (incr != 1 || v_.scale != T(1)) { UnimplementedError() << THROW; }
        v_.m.m.resolve(s, 1, v_.v.v, s, v, s, ' ', v_.m.null_space, v_.m.max_null_space_vector);
        return *this;
    }
 
    Vector& operator=(
          const Vector<
                T,
                MVProd<Minverse_constref<Mcompressed_col_update_inv>, Vincrement_scale_constref> >&
                v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        if (incr != 1 || v_.scale != T(1)) { UnimplementedError() << THROW; }
        v_.m.m.resolve(s, 1, v_.v.v, s, v, s, ' ', v_.m.null_space, v_.m.max_null_space_vector);
        return *this;
    }
 
    Vector& operator=(const Vector<
                      T, MVProd<
                               Minverse_constref<Msym_compressed_col_update_inv_ext>,
                               Vincrement_scale_constref> >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        if (incr != 1 || v_.scale != T(1)) { UnimplementedError() << THROW; }
        v_.m.m.resolve(
              s, 1, v_.v.v, s, v, s, 0, 0, 0, ' ', v_.m.null_space, v_.m.max_null_space_vector);
        return *this;
    }
 
    Vector& operator=(const Vector<
                      T, MVProd<
                               Minverse_constref<Mcompressed_col_update_inv_ext>,
                               Vincrement_scale_constref> >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        if (incr != 1 || v_.scale != T(1)) { UnimplementedError() << THROW; }
        v_.m.m.resolve(
              s, 1, v_.v.v, s, v, s, 0, 0, 0, ' ', v_.m.null_space, v_.m.max_null_space_vector);
        return *this;
    }
 
    Vector& operator=(const Vector<
                      T, MVProd<
                               Minverse_ext_constref<Msym_compressed_col_update_inv_ext>,
                               Vincrement_scale_constref> >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        if (incr != 1 || v_.scale != T(1) || v_.v.scale != T(1)) { UnimplementedError() << THROW; }
        v_.m.m.resolve(
              s, 1, v_.v.v, s, v, s, v_.m.a, v_.m.b, v_.m.c, ' ', v_.m.null_space,
              v_.m.max_null_space_vector);
        return *this;
    }
 
    Vector& operator=(const Vector<
                      T, MVProd<
                               Minverse_ext_constref<Mcompressed_col_update_inv_ext>,
                               Vincrement_scale_constref> >& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        if (incr != 1 || v_.scale != T(1) || v_.v.scale != T(1)) { UnimplementedError() << THROW; }
        v_.m.m.resolve(
              s, 1, v_.v.v, s, v, s, v_.m.a, v_.m.b, v_.m.c, ' ', v_.m.null_space,
              v_.m.max_null_space_vector);
        return *this;
    }
 
    Vector& operator=(
          const Vector<
                T,
                MVProd<Minverse_constref<Mcompressed_col_st_constref>, Vincrement_scale_constref> >&
                v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        if (incr != 1 || v_.scale != T(1) || v_.v.scale != T(1)) { UnimplementedError() << THROW; }
        (*this) = v_.v;
        if (v_.m.trans) {
            for (size_t i = 0; i < size(); ++i) {
                size_t ii = v_.m.si[i + 1] - 1;
                T piv = ((*this)[i] /= v_.m.m[ii]);
                for (size_t j = v_.m.si[i]; j < ii; ++j) {
                    (*this)[v_.m.index[j]] -= v_.m.m[j] * piv;
                }
            }
        } else {
            for (size_t j = size(); j > 0;) {
                size_t jj = v_.m.si[j] - 1;
                T piv = ((*this)[j] /= v_.m.m[jj]);
                for (size_t i = v_.m.si[--j]; i < jj; ++i) {
                    (*this)[v_.m.index[i]] -= v_.m.m[i] * piv;
                }
            }
        }
        return *this;
    }
 
    Vector& operator=(const Vector<T, Vcompressed_scale_constref>& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        std::fill_n(v, s, 0);
        for (size_t i = 0; i != v_.snn; ++i) { v[v_.index[i]] = v_.scale * v_.v[i]; }
        return *this;
    }
 
    Vector& operator=(const Vector<T, Vindex_ref>& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        for (size_t i = 0; i != s; ++i) { v[i] = v_[i]; }
        return *this;
    }
 
    Vector& operator=(const Vector<T, Vindex_constref>& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        for (size_t i = 0; i != s; ++i) { v[i] = v_[i]; }
        return *this;
    }
 
    Vector& operator=(const Vector<T, Vindex_scale_constref>& v_) {
        if (s != v_.size()) {
            Exception() << "The two vectors do not have the same size, " << s << " and "
                        << v_.size() << "." << THROW;
        }
        for (size_t i = 0; i != s; ++i) { v[i] = v_[i]; }
        return *this;
    }
 
    Vector& operator=(
          const Vector<
                T, Sum<Vincrement_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> tmp(v_.v2.v);
        *this = v_.scale * (v_.v1 + v_.v2.scale * (v_.v2.m * tmp));
        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> tmp(v_.v2.v);
        *this = v_.scale * (v_.v1 + v_.v2.scale * (v_.v2.m * tmp));
        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> tmp(v_.v);
        *this = v_.scale * (v_.m * tmp);
        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> tmp(v_.v);
        *this += v_.scale * (v_.m * tmp);
        return *this;
    }
 
    Vector& operator=(const Vector<
                      T, Sum<Vincrement_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> tmp(v_.v2.v);
        *this = v_.scale * (v_.v1 + v_.v2.scale * (v_.v2.m * tmp));
        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> tmp(v_.v2.v);
        *this = v_.scale * (v_.v1 + v_.v2.scale * (v_.v2.m * tmp));
        return *this;
    }
 
    Vector& operator=(const Vector<T, MVProd<Mrectangle_inv, Vincrement_scale_constref> >& v_) {
        if (v_.scale != 1) { UnimplementedError() << THROW; }
        *this = v_.v;
        v_.m.resolve(*this);
        return *this;
    }
 
    static Vector null;
 
    friend logging::Logger& operator<<(logging::Logger& l, const Vector& v) {
        l << "vector ";
        l.write(v.s, v.v, v.incr);
        return l;
    }
 
private:
    size_t s;
    T* v;
    size_t incr;
    MVFRIEND;
};
 
template <typename T>
Vector<T, Vincrement_ref> Vector<T, Vincrement_ref>::null;
 
struct Vincrement_scale_constref {
    typedef Vincrement_scale_constref base;
    typedef Vincrement_scale_constref const_base;
    typedef Vdense copy;
};
 
template <typename T>
class Vector<T, Vincrement_scale_constref> {
public:
    Vector() : s(0), v(0), incr(0), scale(0) {}
 
    Vector(size_t s_, const T* v_, size_t incr_ = 1, T scale_ = T(1))
        : s(s_), v(v_), incr(incr_), scale(scale_) {}
 
    Vector(const Vector& v_) : s(v_.s), v(v_.v), incr(v_.incr), scale(v_.scale) {}
 
    Vector(const Vector<T, Vdense>& v_) : s(v_.size()), v(&v_[0]), incr(1), scale(1) {}
 
    Vector(const Vector<T, Vdense_ref>& v_) : s(v_.s), v(v_.v), incr(1), scale(1) {}
 
    Vector(const Vector<T, Vdense_constref>& v_) : s(v_.s), v(v_.v), incr(1), scale(1) {}
 
    Vector(const Vector<T, Vincrement_ref>& v_) : s(v_.s), v(v_.v), incr(v_.incr), scale(1) {}
 
    Vector(const Vector<T, Vincrement_constref>& v_) : s(v_.s), v(v_.v), incr(v_.incr), scale(1) {}
 
    bool is_null() const { return v == 0; }
 
    size_t size() const { return s; }
 
    T operator[](size_t i) const { return scale * v[i * incr]; }
 
    Vector operator()(const Interval& i) const {
        return Vector(i.size(), &v[i[0] * incr], incr, scale);
    }
 
    Vector operator()(const Slice& i) const {
        return Vector(i.size(), &v[i[0] * incr], i.increment() * incr, scale);
    }
 
    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; }
        return blas::dot(s, v, incr, v_.v, v_.incr) * scale * v_.scale;
    }
 
    Vector& operator*(T scale_) {
        scale *= scale_;
        return *this;
    }
 
    static Vector null;
 
private:
    size_t s;
    const T* v;
    size_t incr;
    T scale;
    MVFRIEND;
};
 
template <typename T>
Vector<T, Vincrement_scale_constref> Vector<T, Vincrement_scale_constref>::null;
 
template <>
inline b2000::csda<double> Vector<b2000::csda<double>, Vincrement_constref>::norm2() const {
    b2000::csda<double> n = 0;
    for (size_t i = 0; i != size(); ++i) { n += std::norm(v[i]); }
    return n;
}
 
}}  // namespace b2000::b2linalg
 
#endif