DenseVector.h

Go to the documentation of this file.

/*  Copyright (C) 2003-2021 Free Electron Organization
    Any use of this software requires a license.  If a valid license
    was not distributed with this file, visit freeelectron.org. */

/** @file */

#ifndef __solve_DenseVector_h__
#define __solve_DenseVector_h__

namespace fe
{
namespace ext
{

/**************************************************************************//**
    @brief Dense vector - size fixed at construction or reset

    @ingroup solve
*//***************************************************************************/
template<class T>
class DenseVector
{
    public:
                    DenseVector(U32 size=0);
                    DenseVector(const DenseVector<T> &other);
                    ~DenseVector(void);

        void        reset(U32 size);
        void        clear(void);

        T           &operator[](U32 index);
        T           operator[](U32 index) const;

        U32         size(void) const    { return m_size; }

        DenseVector<T>  &operator=(const DenseVector<T> &other);

    private:
        void        initialize(U32 size);

        U32 m_size;
        T*  m_pData;
};


template<class T>
inline DenseVector<T>::DenseVector(U32 size):
        m_pData(NULL)
{
    initialize(size);
}

template<class T>
inline void DenseVector<T>::initialize(U32 size)
{
    m_size=size;
    if(m_size)
    {
        m_pData=new T[m_size];
//      feLog("DenseVector<T>::initialize new m_pData %d\n",m_size);
    }

    clear();
}

template<class T>
inline DenseVector<T>::DenseVector(const DenseVector<T> &other):
        m_size(0),
        m_pData(NULL)
{
    operator=(other);
}

template<class T>
inline DenseVector<T>::~DenseVector(void)
{
    delete[] m_pData;
}

template<class T>
inline void DenseVector<T>::reset(U32 size)
{
    if(size==m_size)
    {
        clear();
    }
    else
    {
        delete[] m_pData;
        initialize(size);
    }
}

template<class T>
inline void DenseVector<T>::clear(void)
{
    for(U32 i=0; i<m_size; i++)
    {
        m_pData[i] = T(0);
    }
}

template<>
inline void DenseVector<F32>::clear(void)
{
    if(m_size)
    {
        memset(m_pData,0,m_size*sizeof(F32));
    }
}

template<>
inline void DenseVector<F64>::clear(void)
{
    if(m_size)
    {
        memset(m_pData,0,m_size*sizeof(F64));
    }
}

template<class T>
inline T &DenseVector<T>::operator[](U32 index)
{
#if FE_BOUNDSCHECK
    if(index >= m_size) { feX(e_invalidRange); }
#endif
    return m_pData[index];
}

template<class T>
inline T DenseVector<T>::operator[](U32 index) const
{
#if FE_BOUNDSCHECK
    if(index >= m_size) { feX(e_invalidRange); }
#endif
    return m_pData[index];
}

template<class T>
inline DenseVector<T> &DenseVector<T>::operator=(const DenseVector<T> &other)
{
    if(m_size!=other.m_size)
    {
        m_size=other.m_size;

        delete[] m_pData;
        m_pData=new T[m_size];
//      feLog("operator= new m_pData %d\n",m_size);
    }

    FEASSERT(m_pData);

    if(this != &other)
    {
#if FALSE
        for(U32 i = 0; i < m_size; i++)
        {
            m_pData[i] = other[i];
        }
#else
        memcpy(m_pData,other.m_pData,m_size*sizeof(T));
#endif
    }
    return *this;
}

/* non member operations */

/** Set all the elements to zero
    @relates DenseVector
    */
template<class T>
inline DenseVector<T>& set(DenseVector<T> &lhs)
{
    lhs.clear();
    return lhs;
}

/** Set all the elements to the given value
    @relates DenseVector
    */
template<class T>
inline DenseVector<T>& setAll(DenseVector<T> &lhs,const T value)
{
    const U32 size=lhs.size();
    for(U32 i = 0; i < size; i++)
    {
        lhs[i] = value;
    }
    return lhs;
}

/** Set the value at the index
    @relates DenseVector
    */
template<class T>
inline DenseVector<T> &setAt(DenseVector<T> &lhs, U32 index,const T value)
{
    lhs[index]=value;
    return lhs;
}

/** Add with scaling

    lhs = lhs + rhs * scalar

    @relates DenseVector
    */
template<class T,class U>
inline DenseVector<T> &addScaled(DenseVector<T> &lhs,U scalar,
        const DenseVector<T> &rhs)
{
    const U32 size=lhs.size();
#if FE_BOUNDSCHECK
    if(size != rhs.size()) { feX(e_invalidRange); }
#endif

    for(U32 i = 0; i < size; i++)
    {
        lhs[i] += rhs[i] * scalar;
    }
    return lhs;
}

/** Scale then add

    lhs = lhs * scalar + rhs

    @relates DenseVector
    */
template<class T,class U>
inline DenseVector<T> &scaleAndAdd(DenseVector<T> &lhs,U scalar,
        const DenseVector<T> &rhs)
{
    const U32 size=lhs.size();
#if FE_BOUNDSCHECK
    if(size != rhs.size()) { feX(e_invalidRange); }
#endif

    for(U32 i = 0; i < size; i++)
    {
        lhs[i] = lhs[i] * scalar + rhs[i];
    }
    return lhs;
}

/** In place add operator
    @relates DenseVector
    */
template<class T>
inline DenseVector<T> &operator+=(DenseVector<T> &lhs,
        const DenseVector<T> &rhs)
{
    const U32 size=lhs.size();
#if FE_BOUNDSCHECK
    if(size != rhs.size()) { feX(e_invalidRange); }
#endif

    for(U32 i = 0; i < size; i++)
    {
        lhs[i] += rhs[i];
    }
    return lhs;
}

/** In place subtract operator
    @relates DenseVector
    */
template<class T>
inline DenseVector<T> &operator-=(DenseVector<T> &lhs,
        const DenseVector<T> &rhs)
{
    const U32 size=lhs.size();
#if FE_BOUNDSCHECK
    if(size != rhs.size()) { feX(e_invalidRange); }
#endif

    for(U32 i = 0; i < size; i++)
    {
        lhs[i] -= rhs[i];
    }
    return lhs;
}

/** Negate operation
    @relates DenseVector
    */
template<class T>
inline DenseVector<T> operator-(const DenseVector<T> &rhs)
{
    const U32 size=rhs.size();

    DenseVector<T> v(size);
    for(U32 i = 0; i < size; i++)
    {
        v[i] = -rhs[i];
    }
    return v;
}

/** In place piecewise multiply operator
    @relates DenseVector
    */
template<class T>
inline DenseVector<T> &operator*=(DenseVector<T> &lhs,
        const DenseVector<T> &rhs)
{
    const U32 size=lhs.size();
#if FE_BOUNDSCHECK
    if(size != rhs.size()) { feX(e_invalidRange); }
#endif

    for(U32 i = 0; i < size; i++)
    {
        lhs[i] = lhs[i] * rhs[i];
    }
    return lhs;
}

/** In place piecewise scale operator
    @relates DenseVector
    */
template<class T,class U>
inline DenseVector<T> &operator*=(DenseVector<T> &lhs, U scale)
{
    const U32 size=lhs.size();

    for(U32 i = 0; i < size; i++)
    {
        lhs[i] *= scale;
    }
    return lhs;
}

/** Dot (inner) product
    @relates DenseVector
    */
template<class T>
inline T dot(const DenseVector<T> &lhs, const DenseVector<T> &rhs)
{
    const U32 size=lhs.size();
#if FE_BOUNDSCHECK
    if(size != rhs.size()) { feX(e_invalidRange); }
#endif

    T t = (T)0;
    for(U32 i = 0; i < size; i++)
    {
        t += lhs[i] * rhs[i];
//      feLog("dot %d %.6G*%.6G=%.6G %.6G\n",i,lhs[i],rhs[i],lhs[i]*rhs[i],t);
    }
    return t;
}

/** Frobenius norm operation
    @relates DenseVector
    */
template<class T>
inline T magnitude(const DenseVector<T> &rhs)
{
    const U32 size=rhs.size();

    T mag = 0.0;
    for(U32 i = 0; i < size; i++)
    {
        mag += (T)rhs[i] * (T)rhs[i];
    }
    return sqrt(mag);
}

/** Square of the length
    @relates DenseVector
    */
template<class T>
inline T magnitudeSquared(const DenseVector<T> &rhs)
{
    const U32 size=rhs.size();

    T mag = 0.0;
    for(U32 i = 0; i < size; i++)
    {
        mag += (T)rhs[i] * (T)rhs[i];
    }
    return mag;
}

/** Return normal
    @relates DenseVector
    */
template<class T>
inline DenseVector<T> normal(const DenseVector<T> &rhs)
{
    const U32 size=rhs.size();

    const T mag = magnitude(rhs);
    if(mag == 0.0)
    {
        feX("normal",
            "attempt to normalize zero magnitude vector");
    }
    const T inv = 1.0/mag;
    DenseVector<T> v(size);
    for(U32 i = 0; i < size; i++)
    {
        v[i] = inv * rhs[i];
    }
    return v;
}

/** In place normalize operator
    @relates DenseVector
    */
template<class T>
inline DenseVector<T> &normalize(DenseVector<T> &lhs)
{
    const U32 size=lhs.size();

    const T mag = magnitude(lhs);
    if(mag == 0.0)
    {
        feX("DenseVector<T>::normal",
            "attempt to normalize zero magnitude vector");
    }
    const T inv = 1.0/mag;
    for(U32 i = 0; i < size; i++)
    {
        lhs[i] *= inv;
    }
    return lhs;
}

/** add operation
    @relates DenseVector
    */
template<class T>
inline DenseVector<T> operator+(const DenseVector<T> &lhs,
        const DenseVector<T> &rhs)
{
    const U32 size=lhs.size();
#if FE_BOUNDSCHECK
    if(size != rhs.size()) { feX(e_invalidRange); }
#endif

    DenseVector<T> v(size);
    for(U32 i = 0; i < size; i++)
    {
        v[i] = lhs[i] + rhs[i];
    }
    return v;
}

/** subtractoperation
    @relates DenseVector
    */
template<class T>
inline DenseVector<T> operator-(const DenseVector<T> &lhs,
        const DenseVector<T> &rhs)
{
    const U32 size=lhs.size();
#if FE_BOUNDSCHECK
    if(size != rhs.size()) { feX(e_invalidRange); }
#endif

    DenseVector<T> v(size);
    for(U32 i = 0; i < size; i++)
    {
        v[i] = lhs[i] - rhs[i];
    }
    return v;
}

/** equality test
    @relates DenseVector
    */
template<class T>
inline bool operator==(const DenseVector<T> &lhs, const DenseVector<T> &rhs)
{
    const U32 size=lhs.size();
#if FE_BOUNDSCHECK
    if(size != rhs.size()) { feX(e_invalidRange); }
#endif

    for(U32 i = 0; i < size; i++)
    {
        if(rhs[i] != lhs[i])
        {
            return false;
        }
    }
    return true;
}

/** Equivalence test within the given tolerance @em margin
    @relates DenseVector
    */
template<class T>
inline bool equivalent(const DenseVector<T> &lhs, const DenseVector<T> &rhs,
        T margin)
{
    const U32 size=lhs.size();
#if FE_BOUNDSCHECK
    if(size != rhs.size()) { feX(e_invalidRange); }
#endif

    for(U32 i = 0; i < size; i++)
    {
        if(FE_INVALID_SCALAR(lhs[i]) ||
                FE_INVALID_SCALAR(rhs[i]) ||
                fabs(rhs[i] - lhs[i]) > margin)
        {
            return false;
        }
    }
    return true;
}

/** Piecewise multiply operation
    @relates DenseVector
    */
template<class T>
inline DenseVector<T> operator*(const DenseVector<T> &lhs,
        const DenseVector<T> &rhs)
{
    const U32 size=lhs.size();
#if FE_BOUNDSCHECK
    if(size != rhs.size()) { feX(e_invalidRange); }
#endif

    DenseVector<T> v(size);
    for(U32 i = 0; i < size; i++)
    {
        v[i] = lhs[i] * rhs[i];
    }
    return v;
}

/** Scale operation
    @relates DenseVector
    */
template<class T,class U>
inline DenseVector<T> operator*(const U lhs, const DenseVector<T> &rhs)
{
    const U32 size=rhs.size();

    DenseVector<T> v(size);
    for(U32 i = 0; i < size; i++)
    {
        v[i] = rhs[i] * lhs;
    }
    return v;
}

/** Scale operation
    @relates DenseVector
    */
template<class T,class U>
inline DenseVector<T> operator*(const DenseVector<T> &lhs, const U rhs)
{
    const U32 size=lhs.size();
#if FE_BOUNDSCHECK
    if(size != rhs.size()) { feX(e_invalidRange); }
#endif

    DenseVector<T> v(size);
    for(U32 i = 0; i < size; i++)
    {
        v[i] = lhs[i] * rhs;
    }
    return v;
}

/** Inverse Scale operation
    @relates DenseVector
    */
template<class T,class U>
inline DenseVector<T> operator/(const DenseVector<T> &lhs, const U rhs)
{
    if(rhs==0.0f)
        feX(e_unsolvable,"operator/(DenseVector<T>,T)","divide by zero");
    return lhs*(1.0f/rhs);
}

/** Return number of elements
    @relates DenseVector
    */
template<class T>
inline U32 size(const DenseVector<T> &lhs)
{
    return lhs.size();
}

/** Multiply each pair of components

    @relates DenseVector
    */
template<class T>
inline DenseVector<T>& componentMultiply(DenseVector<T> &result,
        const DenseVector<T> &lhs,const DenseVector<T> &rhs)
{
    const U32 size=lhs.size();
#if FE_BOUNDSCHECK
    if(size != rhs.size()) { feX(e_invalidRange); }
#endif

    for(U32 i = 0; i < size; i++)
    {
        result[i] = lhs[i] * rhs[i];
    }
    return result;
}
} /* namespace ext */
} /* namespace fe */

namespace fe
{

/** Print to a string
    @relates DenseVector
    */
template<class T>
inline String print(const ext::DenseVector<T> &rhs)
{
    const U32 size=rhs.size();

    String s = "";
    for(U32 i = 0; i < size; i++)
    {
        if(i)
        {
            s.cat(" ");
        }
        s.catf("%.6G", rhs[i]);
    }
    return s;
}

} /* namespace fe */

#endif /* __solve_DenseVector_h__ */