b2vector_operation.H

//------------------------------------------------------------------------
// b2vector_operation.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_OPERATION_H__
#define __B2VECTOR_OPERATION_H__
 
#include "b2linear_algebra_def.H"
 
namespace b2000 { namespace b2linalg {
 
template <typename T1, typename T2>
struct Sum {
    typedef Sum const_base;
};
 
template <typename T, typename STORAGE1, typename STORAGE2>
class Vector<T, Sum<STORAGE1, STORAGE2> > {
public:
    Vector(const Vector<T, STORAGE1>& v1_, const Vector<T, STORAGE2>& v2_, T scale_ = T(1))
        : v1(v1_), v2(v2_), scale(scale_) {
        if (v1.size() != v2.size()) {
            Exception() << "The two vectors do not have the same size." << THROW;
        }
    }
 
    size_t size() const { return v1.size(); }
 
    const T& operator[](size_t i) const { return scale * (v1[i] + v2[i]); }
 
    Vector& operator*(T t) {
        scale *= t;
        return *this;
    }
 
private:
    Vector<T, STORAGE1> v1;
    Vector<T, STORAGE2> v2;
    T scale;
    MVFRIEND;
};
 
template <typename T1, typename T2>
struct MVProd {
    typedef MVProd const_base;
};
 
template <typename T, typename STORAGE1, typename STORAGE2>
class Vector<T, MVProd<STORAGE1, STORAGE2> > {
public:
    Vector(const Matrix<T, STORAGE1>& m_, const Vector<T, STORAGE2>& v_, T scale_ = T(1))
        : m(m_), v(v_), scale(scale_) {
        if (m.size2() != v.size()) {
            Exception() << "The number of columns of the matrix (=" << m.size2()
                        << ") differ from the size of the vector (=" << v.size() << ")." << THROW;
        }
    }
 
    size_t size() const { return m.size1(); }
 
    T operator()(size_t i) const { return scale * transpose(m)[i] * v; }
 
    Vector& operator*(T t) {
        scale *= t;
        return *this;
    }
 
private:
    Matrix<T, STORAGE1> m;
    Vector<T, STORAGE2> v;
    T scale;
    MVFRIEND;
};
 
template <typename T1, typename S1, typename T2, typename S2>
inline void copy_v(const Vector<T1, S1>& a, Vector<T2, S2>& b) {
    b.resize(a.size(), T2(0));
    for (size_t i = 0; i != a.size(); ++i) { b[i] = a[i]; }
}
 
template <typename T, typename STORAGE1, typename STORAGE2>
Vector<T, Sum<typename STORAGE1::const_base, typename STORAGE2::const_base> > operator+(
      const Vector<T, STORAGE1>& v1, const Vector<T, STORAGE2>& v2) {
    return Vector<T, Sum<typename STORAGE1::const_base, typename STORAGE2::const_base> >(v1, v2);
}
 
template <typename T, typename STORAGE1, typename STORAGE2>
Vector<T, Sum<typename STORAGE1::const_base, typename STORAGE2::const_base> > operator-(
      const Vector<T, STORAGE1>& v1, const Vector<T, STORAGE2>& v2) {
    return Vector<T, Sum<typename STORAGE1::const_base, typename STORAGE2::const_base> >(v1, -v2);
}
 
template <typename T, typename STORAGE1, typename STORAGE2>
T operator*(const Vector<T, STORAGE1>& v1, const Vector<T, STORAGE2>& v2) {
    return Vector<T, typename STORAGE1::const_base>(v1)
           * Vector<T, typename STORAGE2::const_base>(v2);
}
 
template <typename T, typename STORAGE1, typename STORAGE2>
Vector<T, MVProd<typename STORAGE1::const_base, typename STORAGE2::const_base> > operator*(
      const Matrix<T, STORAGE1>& m, const Vector<T, STORAGE2>& v) {
    return Vector<T, MVProd<typename STORAGE1::const_base, typename STORAGE2::const_base> >(m, v);
}
 
template <typename T, typename STORAGE1, typename STORAGE2>
Vector<T, MVProd<typename STORAGE1::const_base, typename STORAGE2::const_base> > operator*(
      const Vector<T, STORAGE2>& v, const Matrix<T, STORAGE1>& m) {
    return Vector<T, MVProd<typename STORAGE1::const_base, typename STORAGE2::const_base> >(
          transposed(m), v);
}
 
template <typename T, typename STORAGE>
Vector<T, typename STORAGE::const_base> operator*(T t, const Vector<T, STORAGE>& v) {
    return Vector<T, typename STORAGE::const_base>(v) * t;
}
 
template <typename T, typename STORAGE>
Vector<T, typename STORAGE::const_base> operator*(const Vector<T, STORAGE>& v, T t) {
    return Vector<T, typename STORAGE::const_base>(v) * t;
}
 
template <typename T, typename STORAGE>
Vector<T, typename STORAGE::const_base> operator-(const Vector<T, STORAGE>& v) {
    return Vector<T, typename STORAGE::const_base>(v) * T(-1);
}
 
template <typename T, typename TYPE>
std::ostream& operator<<(std::ostream& out, const Vector<T, TYPE>& v) {
    out << "( ";
    for (size_t i = 0; i != v.size(); ++i) { out << v[i] << ", "; }
    out << ")";
    return out;
}
 
template <typename T>
void copy(const Vector<T, Vdense_constref>& v1, Vector<T, Vdense>& v2) {
    v2 = v1;
}
 
template <typename T>
void copy(const Vector<b2000::csda<T>, Vdense_constref>& v1, Vector<T, Vdense>& v2) {
    v2.resize(v1.size());
    for (size_t i = 0; i != v2.size(); ++i) { v2[i] = v1[i].real(); }
}
 
template <typename T>
void copy(const Vector<std::complex<T>, Vdense_constref>& v1, Vector<T, Vdense>& v2) {
    v2.resize(v1.size());
    for (size_t i = 0; i != v2.size(); ++i) { v2[i] = v1[i].real(); }
}
 
template <typename T>
void copy(const Vector<T, Vdense>& v1, Vector<T, Vdense>& v2) {
    v2 = v1;
}
 
template <typename T>
void copy(const Vector<T, Vdense>& v1, Vector<b2000::csda<T>, Vdense>& v2) {
    v2.resize(v1.size());
    for (size_t i = 0; i != v2.size(); ++i) { v2[i] = b2000::csda<T>(v1[i]); }
}
 
template <typename T>
void copy(const Vector<T, Vdense>& v1, Vector<std::complex<T>, Vdense>& v2) {
    v2.resize(v1.size());
    for (size_t i = 0; i != v2.size(); ++i) { v2[i] = std::complex<T>(v1[i]); }
}
 
template <typename T>
inline Vector<T, Vincrement_ref> real(std::vector<b2000::csda<T> >& v) {
    return Vector<T, Vincrement_ref>(v.size(), reinterpret_cast<T*>(&v[0]), 2);
}
 
template <typename T>
inline Vector<T, Vincrement_ref> real(std::vector<std::complex<T> >& v) {
    return Vector<T, Vincrement_ref>(v.size(), reinterpret_cast<T*>(&v[0]), 2);
}
 
template <typename T>
inline Vector<T, Vincrement_ref> imag(std::vector<b2000::csda<T> >& v) {
    return Vector<T, Vincrement_ref>(v.size(), reinterpret_cast<T*>(&v[0]) + 1, 2);
}
 
template <typename T>
inline Vector<T, Vincrement_ref> imag(std::vector<std::complex<T> >& v) {
    return Vector<T, Vincrement_ref>(v.size(), reinterpret_cast<T*>(&v[0]) + 1, 2);
}
 
template <typename T>
inline Vector<T, Vincrement_constref> real(const std::vector<b2000::csda<T> >& v) {
    return Vector<T, Vincrement_ref>(v.size(), reinterpret_cast<T*>(&v[0]), 2);
}
 
template <typename T>
inline Vector<T, Vincrement_constref> real(const std::vector<std::complex<T> >& v) {
    return Vector<T, Vincrement_ref>(v.size(), reinterpret_cast<T*>(&v[0]), 2);
}
 
template <typename T>
inline Vector<T, Vincrement_constref> imag(const std::vector<b2000::csda<T> >& v) {
    return Vector<T, Vincrement_ref>(v.size(), reinterpret_cast<T*>(&v[0]) + 1, 2);
}
 
template <typename T>
inline Vector<T, Vincrement_constref> imag(const std::vector<std::complex<T> >& v) {
    return Vector<T, Vincrement_ref>(v.size(), reinterpret_cast<T*>(&v[0]) + 1, 2);
}
 
}}  // namespace b2000::b2linalg
 
#endif