BCGThreaded.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_BCGThreaded_h__
#define __solve_BCGThreaded_h__

#define BCGT_DEBUG          FALSE
#define BCGT_TRACE          FALSE
#define BCGT_THREAD_DEBUG   FALSE
#define BCGT_VERIFY         (FE_CODEGEN<=FE_DEBUG)

namespace fe
{
namespace ext
{

/**************************************************************************//**
    @brief solve Ax=b for x

    @ingroup solve

    Uses Biconjugate-Gradient.  The matrix must be positive-definite,
    but not necessarily symmetric.  For symmetric matrices, a regular
    Conjugate-Gradient can be about twice as fast.

    The arguments are templated, so any argument types should work,
    given that they have the appropriate methods and operators.
*//***************************************************************************/
template <typename MATRIX, typename VECTOR>
class BCGThreaded
{
    public:

    class PointerCounted: public Counted
    {
        public:
                            PointerCounted(BCGThreaded* pointer)
                            {   m_pointer=pointer; }
            BCGThreaded*    m_pointer;
    };

    class Worker
    {
        public:
                            Worker(U32 id,sp< JobQueue<I32> > spJobQueue):
                                m_id(id),
                                m_spJobQueue(spJobQueue)                    {}

            void            operator()(void);

            void            setObject(sp<Counted> spObject)
                            {
                                m_spPointerCounted=spObject;
                                FEASSERT(m_spPointerCounted.isValid());
                                m_pBCGThreaded=m_spPointerCounted->m_pointer;
                            }

        private:
            U32                 m_id;
            sp< JobQueue<I32> > m_spJobQueue;
            sp<PointerCounted>  m_spPointerCounted;
            BCGThreaded*        m_pBCGThreaded;
    };

                BCGThreaded(void):
                        m_A(NULL),
                        m_x(NULL),
                        m_b(NULL),
                        m_threshold(1e-6f)
                {
                    sp<PointerCounted> spPointerCounted(
                            new PointerCounted(this));
                    m_spGang=new Gang<Worker,I32>();
                    m_spGang->start(2,spPointerCounted);
                }
virtual         ~BCGThreaded(void)
                {
                    m_spGang->post(-1);
                    m_spGang->post(-1);
                    m_spGang->finish();
                }

        void    solve(VECTOR& x, const MATRIX& A, const VECTOR& b);
        void    setThreshold(F64 threshold) { m_threshold=threshold; }

    private:
        void    solve(U32 thread);
        void    solve(U32 thread, VECTOR& x, const MATRIX& A,const VECTOR& b);

        sp< Gang<Worker,I32> >  m_spGang;

        const MATRIX*   m_A;
        VECTOR* m_x;
        const VECTOR*   m_b;

        VECTOR  r,r_1,r_2;      //* residual    (at k, k-1, k-2)
        VECTOR  rb,rb_1,rb_2;   //* second residual (r bar)
        VECTOR  p,p_1;          //* direction
        VECTOR  pb,pb_1;        //* second direction

        VECTOR  temp[2];        //* persistent temporary

        VECTOR  Ap;
        F64     m_threshold;
        F64     m_dot_r_1;
        F64     m_alpha;
        F64     m_beta;
        U32     m_N;
        BWORD   m_break;
};

template <typename MATRIX, typename VECTOR>
inline void BCGThreaded<MATRIX,VECTOR>::Worker::operator()(void)
{
    while(TRUE)
    {
        I32 job;
#if BCGT_THREAD_DEBUG
        feLog("BCGThreaded<>::Worker::operator %p thread %d wait\n",
                m_spJobQueue.raw(),m_id);
#endif
        m_spJobQueue->waitForJob(job);
        if(job<0)
        {
#if BCGT_THREAD_DEBUG
            feLog("BCGThreaded<>::Worker::operator %p thread %d break\n",
                m_spJobQueue.raw(),m_id);
#endif
            break;
        }
#if BCGT_THREAD_DEBUG
        feLog("BCGThreaded<>::Worker::operator %p thread %d solve %d\n",
                m_spJobQueue.raw(),m_id,job);
#endif
        m_pBCGThreaded->solve(job);
#if BCGT_THREAD_DEBUG
        feLog("BCGThreaded<>::Worker::operator %p thread %d deliver %d\n",
                m_spJobQueue.raw(),m_id,job);
#endif
        m_spJobQueue->deliver(job);
    }
}

template <typename MATRIX, typename VECTOR>
inline void BCGThreaded<MATRIX,VECTOR>::solve(VECTOR& x,
        const MATRIX& A, const VECTOR& b)
{
    m_A=&A;
    m_x=&x;
    m_b=&b;

#if BCGT_THREAD_DEBUG
    feLog("BCGThreaded<>::solve %p post\n",m_spGang.raw());
#endif
    m_spGang->post(0,1);
    U32 jobs=2;
    I32 job;
    while(jobs)
    {
        if(m_spGang->acceptDelivery(job))
        {
#if BCGT_THREAD_DEBUG
            feLog("BCGThreaded<>::solve %p acceptDelivery %d\n",
                    m_spGang.raw(),job);
#endif
            jobs--;
        }
    }

    m_A=NULL;
    m_x=NULL;
    m_b=NULL;

#if BCGT_THREAD_DEBUG
    feLog("BCGThreaded<>::solve %p complete\n",m_spGang.raw());
#endif
}

template <typename MATRIX, typename VECTOR>
inline void BCGThreaded<MATRIX,VECTOR>::solve(U32 thread)
{
    solve(thread,*m_x,*m_A,*m_b);
}

template <typename MATRIX, typename VECTOR>
inline void BCGThreaded<MATRIX,VECTOR>::solve(U32 thread,VECTOR& x,
        const MATRIX& A, const VECTOR& b)
{
    if(!thread)
    {
        m_N=size(b);
        if(size(x)!=m_N)
        {
            x=b;        // adopt size
        }
        if(size(Ap)!=m_N)
        {
            Ap=b;       // adopt size
        }
        set(x);
        temp[1]=x;
        m_break=FALSE;

#if BCGT_DEBUG
        feLog("\nA\n%s\nb=<%s>\n",c_print(A),c_print(b));
#endif

        if(magnitudeSquared(b)<m_threshold)
        {
#if BCGT_DEBUG
            feLog("BCGThreaded::solve has trivial solution\n");
#endif
            m_break=TRUE;
        }
    }

    m_spGang->synchronize();
    if(m_break)
    {
        return;
    }

    for(U32 k=1;k<=m_N;k++)
    {
        if(!thread)
        {
            if(k==1)
            {
                p=b;
                p_1=p;

                set(r);
                r_1=p;
                r_2=r_1;

                pb=b;
                rb_1=pb;
                set(rb_2);
                set(pb_1);

                m_dot_r_1=dot(rb_1,r_1);
                m_beta=0.0f;
            }
            else
            {
                r_2=r_1;
                r_1=r;
                p_1=p;

                rb_2=rb_1;
                rb_1=rb;
                pb_1=pb;

                m_dot_r_1=dot(rb_1,r_1);
                m_beta=m_dot_r_1/dot(rb_2,r_2);

//              p=r_1+m_beta*p_1;
                temp[0]=p_1;
                temp[0]*=m_beta;
                p=r_1;
                p+=temp[0];

//              pb=rb_1+m_beta*pb_1;
                temp[0]=pb_1;
                temp[0]*=m_beta;
                pb=rb_1;
                pb+=temp[0];
            }
        }

        m_spGang->synchronize();
        if(!thread)
        {
            transformVector(A,p,Ap);
            m_alpha=m_dot_r_1/dot(pb,Ap);

            BWORD zero_mag=(magnitudeSquared(p)==0.0f);

#if BCGT_TRACE==FALSE
            if(zero_mag)
#endif
            {
                feLog("\n%d m_alpha=%.6G beta=%.6G\n",k,m_alpha,m_beta);
                feLog("x<%s>\n",c_print(x));
                feLog("r<%s>\n",c_print(r));
                feLog("r_1<%s>\n",c_print(r_1));
                feLog("r_2<%s>\n",c_print(r_2));
                feLog("rb_1<%s>\n",c_print(rb_1));
                feLog("rb_2<%s>\n",c_print(rb_2));
                feLog("p<%s>\n",c_print(p));
                feLog("p_1<%s>\n",c_print(p_1));
                feLog("A*p<%s>\n",c_print(A*p));
                feLog("pb<%s>\n",c_print(pb));
                feLog("pb_1<%s>\n",c_print(pb_1));
            }

            if(zero_mag)
            {
                feX("BCGThreaded::solve","direction lost its magnitude");
            }

//          x+=m_alpha*p;
            temp[0]=p;
            temp[0]*=m_alpha;
            x+=temp[0];

//          r=r_1-m_alpha*(Ap);
            temp[0]=Ap;
            temp[0]*=m_alpha;
            r=r_1;
            r-=temp[0];

            if(magnitudeSquared(r)<m_threshold)
            {
#if BCGT_DEBUG
                feLog("BCGThreaded::solve early solve %d/%d\n",k,m_N);
#endif
                m_break=TRUE;
                m_spGang->synchronize();
                break;
            }

            if(k==m_N)
            {
                feLog("BCGThreaded::solve ran %d/%d\n",k,m_N);
            }

            m_spGang->synchronize();

            temp[1]*=m_alpha;
            rb=rb_1;
            rb-=temp[1];

#if BCGT_TRACE
            feLog("temp[1]<%s>\n",c_print(temp[1]));
            feLog("rb<%s>\n",c_print(rb));
#endif
        }
        else
        {
//          rb=rb_1-m_alpha*transposeMultiply(A,pb);
            transposeTransformVector(A,pb,temp[1]);

            m_spGang->synchronize();

            if(m_break)
            {
                break;
            }
        }
    }

    if(!thread)
    {
#if BCGT_DEBUG
        feLog("\nx=<%s>\nA*x=<%s>\n",c_print(x),c_print(A*x));
#endif

#if BCGT_VERIFY
        BWORD invalid=FALSE;
        for(U32 k=0;k<m_N;k++)
        {
            if(FE_INVALID_SCALAR(x[k]))
            {
                invalid=TRUE;
            }
        }
        VECTOR Ax=A*x;
        F64 distance=magnitude(Ax-b);
        if(invalid || distance>1.0f)
        {
            feLog("BCGThreaded::solve failed to converge (dist=%.6G)\n",
                    distance);
            if(size(x)<100)
            {
                feLog("  collecting state ...\n");
                feLog("A=\n%s\nx=<%s>\nA*x=<%s>\nb=<%s>\n",
                        c_print(A),c_print(x),
                        c_print(Ax),c_print(b));
            }
    //      feX("BCGThreaded::solve","failed to converge");
        }
#endif
    }
}

} /* namespace ext */
} /* namespace fe */

#endif /* __solve_BCGThreaded_h__ */