SparseMatrix2x2.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 __SparseMatrix2x2_h__
#define __SparseMatrix2x2_h__

namespace fe
{
namespace ext
{

template <typename T>
class FE_DL_EXPORT SparseMatrix2x2 : public Component
{
    public:
typedef T                   t_real;
typedef Vector<2, t_real>   t_v2;
typedef Matrix<2,2,t_real>  t_matrix;

virtual ~SparseMatrix2x2(void) {}
virtual void multiply(std::vector<t_v2> &a_b,
            const std::vector<t_v2> &a_x) const                         = 0;
virtual void scale(t_real a_scale)                                      = 0;
};


template <typename T>
class FE_DL_EXPORT FullBCRS2 : public SparseMatrix2x2<T>
{
    public:
typedef T                   t_real;
typedef Vector<2, t_real>   t_v2;
typedef Matrix<2,2,t_real>  t_matrix;
        FullBCRS2(void);
virtual ~FullBCRS2(void);


virtual void multiply(std::vector<t_v2> &a_b,
            const std::vector<t_v2> &a_x) const;
virtual void scale(t_real a_scale);

        void                clear(void);
        void                startrow(void);
virtual void                done(void);
        t_matrix        &insert(unsigned int a_columnIndex,
                                const t_matrix &a_block);

        void                sparse(const sp<FullBCRS2> &a_other);
        void                project(const sp<FullBCRS2> &a_other);

    public:
        std::vector<t_matrix>   &blocks(void) { return m_blocks; }
        std::vector<unsigned int>   &colind(void) { return m_colind; }
        std::vector<unsigned int>   &rowptr(void) { return m_rowptr; }

    protected:
        std::vector<t_matrix>           m_blocks;
        std::vector<unsigned int>       m_colind;
        std::vector<unsigned int>       m_rowptr;
        unsigned int                    m_r;
        t_v2                            m_zeroVector;
};

template <typename T>
class FE_DL_EXPORT UpperTriangularBCRS2 : public FullBCRS2<T>
{
    public:
typedef T                   t_real;
typedef Vector<2, t_real>   t_v2;
typedef Matrix<2,2,t_real>  t_matrix;
        UpperTriangularBCRS2(void){}
virtual ~UpperTriangularBCRS2(void){}

virtual void multiply(std::vector<t_v2> &a_b,
            const std::vector<t_v2> &a_x) const;

        void incompleteSqrt(void);
        void incompleteSqrtOpt(void);
        void backSub(std::vector<t_v2> &a_x, const std::vector<t_v2> &a_b);
        void transposeForeSub(std::vector<t_v2> &a_x, const std::vector<t_v2> &a_b);
        void backSolve(std::vector<t_v2> &a_x, const std::vector<t_v2> &a_b);
virtual void done(void);

    private:
        std::vector<unsigned int>   m_ptr; // buffer for transposeForeSub
        std::vector<t_v2>           m_y; // buffer for backSolve

};



template <typename T>
inline FullBCRS2<T>::FullBCRS2(void)
{
    m_r = 0;
    m_zeroVector = t_v2(0.0,0.0);
}

template <typename T>
inline FullBCRS2<T>::~FullBCRS2(void)
{
}

template <typename T>
inline void FullBCRS2<T>::clear(void)
{
    m_blocks.clear();
    m_colind.clear();
    m_rowptr.clear();
    m_r = 0;
}

template <typename T>
inline void FullBCRS2<T>::startrow(void)
{
    m_rowptr.push_back(m_r);
}

template <typename T>
inline void FullBCRS2<T>::done(void)
{
    m_rowptr.push_back(m_r);
}

template <typename T>
inline void UpperTriangularBCRS2<T>::done(void)
{
    m_y.resize(FullBCRS2<T>::m_rowptr.size());
    FullBCRS2<T>::m_rowptr.push_back(FullBCRS2<T>::m_r);
    m_ptr.resize(FullBCRS2<T>::m_rowptr.size());
}

template <typename T>
inline Matrix<2,2,T> &FullBCRS2<T>::insert(unsigned int a_columnIndex, const t_matrix &a_block)
{
    m_blocks.push_back(a_block);
    m_colind.push_back(a_columnIndex);
    m_r++;
    return m_blocks.back();
}

template <typename T>
inline void FullBCRS2<T>::multiply(std::vector<t_v2> &a_b, const std::vector<t_v2> &a_x) const
{
    unsigned int n = (unsigned int)(m_rowptr.size() - 1);
    for(unsigned int i = 0; i < n; i++)
    {
        a_b[i] = m_zeroVector;
        for(unsigned int k = m_rowptr[i]; k < m_rowptr[i+1]; k++)
        {
            a_b[i] += fe::multiply(m_blocks[k], a_x[m_colind[k]]);
        }
    }
}

template <typename T>
inline void FullBCRS2<T>::scale(t_real a_scale)
{
    unsigned int n = (unsigned int)m_blocks.size();
    for(unsigned int i = 0; i < n; i++)
    {
        m_blocks[i] *= a_scale;
    }
}

template <typename T>
inline void UpperTriangularBCRS2<T>::multiply(std::vector<t_v2> &a_b, const std::vector<t_v2> &a_x) const
{
    unsigned int n = (unsigned int)(FullBCRS2<T>::m_rowptr.size() - 1);
    for(unsigned int i = 0; i < n; i++)
    {
        a_b[i] = FullBCRS2<T>::m_zeroVector;
    }
    for(unsigned int i = 0; i < n; i++)
    {
        unsigned int k = FullBCRS2<T>::m_rowptr[i];
        a_b[i] += fe::multiply(FullBCRS2<T>::m_blocks[k], a_x[i]);
        for(k=k+1; k < FullBCRS2<T>::m_rowptr[i+1]; k++)
        {
            a_b[i] += fe::multiply(FullBCRS2<T>::m_blocks[k], a_x[FullBCRS2<T>::m_colind[k]]);
            a_b[FullBCRS2<T>::m_colind[k]] += fe::transposeMultiply(FullBCRS2<T>::m_blocks[k], a_x[i]);
        }
    }
}

template <typename T>
inline void just_upper(Matrix<2,2,T> &a_m)
{
    a_m(1,0) = 0.0;
}


template <typename T>
inline void FullBCRS2<T>::sparse(const sp<FullBCRS2> &a_other)
{
    m_blocks.clear();
    m_colind.clear();
    m_rowptr.clear();
    m_r = 0;
    t_matrix zero;
    setAll(zero, 0.0);
    unsigned int n = a_other->m_rowptr.size()-1;
    for(unsigned int i = 0; i < n; i++)
    {
        startrow();
        int d = a_other->m_rowptr[i];
        for(int k=d; k < a_other->m_rowptr[i+1]; k++)
        {
            int c = a_other->m_colind[k];
            if(!(zero == a_other->m_blocks[k]))
            {
                insert(c, a_other->m_blocks[k]);
            }
        }
    }
    done();
}

template <typename T>
inline void FullBCRS2<T>::project(const sp<FullBCRS2> &a_other)
{
    unsigned int n = a_other->m_rowptr.size()-1;
    for(unsigned int i = 0; i < n; i++)
    {
        int d = a_other->m_rowptr[i];
        int d_self = m_rowptr[i];
        for(int k=d; k < a_other->m_rowptr[i+1]; k++)
        {
            int c = a_other->m_colind[k];
            for(int k_self=d_self; k_self < m_rowptr[i+1]; k_self++)
            {
                int c_self = m_colind[k_self];
                if(c_self == c)
                {
                    m_blocks[k_self] = a_other->m_blocks[k];
                }
            }
        }
    }
}

template <typename T>
inline void UpperTriangularBCRS2<T>::backSub(std::vector<t_v2> &a_x, const std::vector<t_v2> &a_b)
{
#if 0
t_matrix zero;
setAll(zero, 0.0);
#endif
    unsigned int n = (unsigned int)(FullBCRS2<T>::m_rowptr.size()-1);
//fe_fprintf(stderr, "root %d %d\n", n-1, m_colind[m_rowptr[n-1]]);
    fe::backSub(FullBCRS2<T>::m_blocks[FullBCRS2<T>::m_rowptr[n-1]], a_x[n-1], a_b[n-1]);
    t_v2 tmp;

    for(int i = n-2; i >= 0; i--)
    {
//fe_fprintf(stderr, "bs i %d\n", i);
        a_x[i] = a_b[i];
        int d = FullBCRS2<T>::m_rowptr[i];
        for(unsigned int k=d+1; k < FullBCRS2<T>::m_rowptr[i+1]; k++)
        {
            int c = FullBCRS2<T>::m_colind[k];
            fe::multiply(tmp, FullBCRS2<T>::m_blocks[k], a_x[c]);
            a_x[i] -= tmp;
#if 0
if(! (zero == m_blocks[k]))
{
fe_fprintf(stderr, "bs c %d -- %d <- T %d\n", c, i, k);
fprint(stderr, m_blocks[k]);
fe_fprintf(stderr, "\n");
}
#endif
            //a_x[i] -= fe::transposeMultiply(m_blocks[k], a_x[c]);
        }
        fe::backSub(FullBCRS2<T>::m_blocks[d], a_x[i], a_x[i]);
//fprint(stderr, m_blocks[d]);
//fe_fprintf(stderr, "\n");
    }
}

// operates on upper tri as if it is a transpose of the upper and is lower
template <typename T>
inline void UpperTriangularBCRS2<T>::transposeForeSub(std::vector<t_v2> &a_x, const std::vector<t_v2> &a_b)
{
//t_matrix zero;
//setAll(zero, 0.0);
    unsigned int n = (unsigned int)(FullBCRS2<T>::m_rowptr.size()-1);
    //std::vector<unsigned int> ptr = FullBCRS2<T>::m_rowptr;
    for(unsigned int i = 0 ; i < FullBCRS2<T>::m_rowptr.size(); i++)
    {
        m_ptr[i] = FullBCRS2<T>::m_rowptr[i];
    }
    fe::transposeForeSub(FullBCRS2<T>::m_blocks[FullBCRS2<T>::m_rowptr[0]], a_x[0], a_b[0]);
    t_v2 tmp;

    for(unsigned int i = 1; i < n; i++)
    {
//fe_fprintf(stderr, "fs i %d\n", i);
        a_x[i] = a_b[i];
        for(unsigned int m = 0; m < i; m++)
        {
            unsigned int c = FullBCRS2<T>::m_colind[m_ptr[m]];
            while(c < i && m_ptr[m] < FullBCRS2<T>::m_rowptr[m+1]-1)
            {
                m_ptr[m]++;
                c = FullBCRS2<T>::m_colind[m_ptr[m]];
            }
            if(c == i)
            {
//if(! (zero == m_blocks[ptr[m]]))
//{
//fe_fprintf(stderr, "fs m %d -- i %d <- T ptr[m] %d m_rowptr[m+1] %d\n", m, i, ptr[m], m_rowptr[m+1]);
//fprint(stderr, m_blocks[ptr[m]]);
//fe_fprintf(stderr, "\n");
//}
                //fe::multiply(tmp, m_blocks[ptr[m]], a_x[m]);
                //a_x[i] -= tmp;
                a_x[i] -= fe::transposeMultiply(FullBCRS2<T>::m_blocks[m_ptr[m]], a_x[m]);
            }
        }
//fe_fprintf(stderr, "fs m_rowptr[%d] %d col %d\n", i, m_rowptr[i], m_colind[m_rowptr[i]]);
        fe::transposeForeSub(FullBCRS2<T>::m_blocks[FullBCRS2<T>::m_rowptr[i]], a_x[i], a_x[i]);
//fprint(stderr, m_blocks[m_rowptr[i]]);
//fe_fprintf(stderr, "\n");
    }
}

template <typename T>
inline void UpperTriangularBCRS2<T>::backSolve(std::vector<t_v2> &a_x, const std::vector<t_v2> &a_b)
{
    //std::vector<t_v2> y(a_x.size());

    transposeForeSub(m_y, a_b);
    backSub(a_x, m_y);
}


template <typename T>
inline void UpperTriangularBCRS2<T>::incompleteSqrt(void)
{
    unsigned int d;
    unsigned int n = (unsigned int)(FullBCRS2<T>::m_rowptr.size()-1);
    t_matrix inv_z;
    t_matrix z;
    for(unsigned int k = 0; k < n-1; k++)
    {
        d = FullBCRS2<T>::m_rowptr[k];
        just_upper(FullBCRS2<T>::m_blocks[d]);
        //t_matrix tz = m_blocks[d] + transpose(m_blocks[d]);
        //fprint(stderr, m_blocks[d]); fe_fprintf(stderr, "\n");
        //nanCheck(m_blocks[d], "pre squareroot");
        squareroot(FullBCRS2<T>::m_blocks[d], FullBCRS2<T>::m_blocks[d]);
        //nanCheck(m_blocks[d], "post squareroot");
        z = FullBCRS2<T>::m_blocks[d];
        //fprint(stderr, m_blocks[d]); fe_fprintf(stderr, "\n");

        just_upper(z);

        //t_v2 sv(0.2,0.5,0.7);
        //t_v2 x, y;

        invert(inv_z, z);

        for(unsigned int i=d+1; i < FullBCRS2<T>::m_rowptr[k+1]; i++)
        {
            //fe::multiply(m_blocks[i], inv_z, m_blocks[i]);
            //fprint(stderr, m_blocks[i]); fe_fprintf(stderr, "\n");
            //fprint(stderr, inv_z); fe_fprintf(stderr, "\n");

            //nanCheck(m_blocks[i], "pre multiply A");
            //nanCheck(inv_z, "pre multiply inv_z");
            fe::multiply(z, FullBCRS2<T>::m_blocks[i], inv_z);
            setTranspose(z);

            //fe::multiply(m_blocks[i], inv_z, m_blocks[i]);
            //setTranspose(m_blocks[i]);

            //fprint(stderr, m_blocks[i]); fe_fprintf(stderr, "\n");
            //fprint(stderr, z); fe_fprintf(stderr, "\n");
            FullBCRS2<T>::m_blocks[i] = z;
        }

        for(unsigned int i=d+1; i < FullBCRS2<T>::m_rowptr[k+1]; i++)
        {
            z = transpose(FullBCRS2<T>::m_blocks[i]);
            unsigned int h = FullBCRS2<T>::m_colind[i];
            unsigned int g = i;

            t_matrix t;

            for(unsigned int j = FullBCRS2<T>::m_rowptr[h]; j < FullBCRS2<T>::m_rowptr[h+1]; j++)
            {
                for(; g < FullBCRS2<T>::m_rowptr[k+1] && FullBCRS2<T>::m_colind[g] <= FullBCRS2<T>::m_colind[j]; g++)
                {
                    if (FullBCRS2<T>::m_colind[g] == FullBCRS2<T>::m_colind[j])
                    {
#if 0
                        fe::multiply(t, z, FullBCRS2<T>::m_blocks[g]);
                        setTranspose(t);
                        subtract(FullBCRS2<T>::m_blocks[j], FullBCRS2<T>::m_blocks[j], t);
#endif

#if 1
                        const t_matrix &A = z;
                        const t_matrix &B = FullBCRS2<T>::m_blocks[g];
                        for(unsigned int ii = 0;ii < 2; ii++)
                        {
                            for(unsigned int jj = 0;jj < 2; jj++)
                            {
                                t_real sum=0.0;
                                sum += A(ii,0)*B(0,jj);
                                sum += A(ii,1)*B(1,jj);
                                t(ii,jj)=sum;
                            }
                        }

                        t_matrix &R = FullBCRS2<T>::m_blocks[j];
                        R(0,0) = R(0,0)-t(0,0);
                        R(0,1) = R(0,1)-t(1,0);
                        R(1,0) = R(1,0)-t(0,1);
                        R(1,1) = R(1,1)-t(1,1);
#endif
                    }
                }
            }
        }
    }
    d = FullBCRS2<T>::m_rowptr[n-1];
    t_matrix tmp;
    just_upper(FullBCRS2<T>::m_blocks[d]);
    squareroot(tmp, FullBCRS2<T>::m_blocks[d]);
    FullBCRS2<T>::m_blocks[d] = tmp;
}



} /* namespace */
} /* namespace */

#endif /* __SparseMatrix2x2_h__ */