MatrixPower.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 __geometry_MatrixPower_h__
#define __geometry_MatrixPower_h__

#define MRP_DEBUG       FALSE
#define MRP_VALIDATE    (FE_CODEGEN<=FE_DEBUG)

namespace fe
{
namespace ext
{


// NOTE fraction: 23 bits for single and 52 for double
/**************************************************************************//**
    @brief solve B = A^^power, where A is a matrix

    @ingroup geometry

    The power can be any arbitrary real number.

    Execution time is roughly proportional to the number of set bits in
    the integer portion of the floating point power and a fixed number
    of iterations for the fractional part.

    The number of iterations used to compute of the fractional portion
    of the power can be changed.  The maximum error after each iteration
    is half of the previous iteration, starting with one half.  The entire
    integer portion of the power is always computed.
*//***************************************************************************/
template <typename MATRIX>
class MatrixPower
{
    public:
                MatrixPower(void):
                    m_iterations(16)                                        {}

                template <typename T>
        void    solve(MATRIX& B, const MATRIX& A, T a_power) const;

        void    setIterations(U32 iterations)   { m_iterations=iterations; }

    private:
        MatrixSqrt<MATRIX>  m_matrixSqrt;
        U32                 m_iterations;
};

template <typename MATRIX>
template <typename T>
inline void MatrixPower<MATRIX>::solve(MATRIX& B, const MATRIX& A,
        T a_power) const
{
    T absolute=a_power;

#if MRP_DEBUG
    feLog("\nA\n%s\npower=%.6G\n",print(A).c_str(),absolute);
#endif

    const BWORD inverted=(absolute<0.0);
    if(inverted)
    {
        absolute= -absolute;
    }

    U32 whole=U32(absolute);
    F32 fraction=absolute-whole;

#if MRP_DEBUG
    feLog("\nwhole=%d\nfraction=%.6G\n",whole,fraction);
#endif

    MATRIX R;
    setIdentity(R);

    MATRIX partial=A;
    F32 contribution=1.0;
    U32 iteration;
    for(iteration=0;iteration<m_iterations;iteration++)
    {
        m_matrixSqrt.solve(partial,partial);
        contribution*=0.5;
#if MRP_DEBUG
    feLog("\ncontribution=%.6G\nfraction=%.6G\n",contribution,fraction);
#endif
        if(fraction>=contribution)
        {
            R*=partial;
            fraction-=contribution;
        }
    }

    partial=A;
    while(whole)
    {
#if MRP_DEBUG
    feLog("\nwhole=%d\n",whole);
#endif
        if(whole&1)
        {
            R*=partial;
        }
        whole>>=1;
        if(whole)
        {
            partial*=partial;
        }
    }

#if MRP_VALIDATE
    BWORD invalid=FALSE;
    for(U32 m=0;m<width(R);m++)
    {
        for(U32 n=0;n<height(R);n++)
        {
            if(FE_INVALID_SCALAR(R(m,n)))
            {
                invalid=TRUE;
            }
        }
    }
    if(invalid)
    {
        feLog("MatrixPower< %s >::solve invalid results power=%.6G\n",
                FE_TYPESTRING(MATRIX).c_str(),a_power);
        feLog("\nA\n%s\n",print(A).c_str());
        feLog("\nB\n%s\n",print(R).c_str());

        feX("MatrixPower<>::solve","invalid result");
    }
#endif

    if(inverted)
    {
        invert(B,R);
    }
    else
    {
        B=R;
    }
}

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

#endif /* __geometry_MatrixPower_h__ */