Eigen3 dynamic matrix with boost::multiprecision::mpfr_float

Question

I would like to make matrices and use them using the Eigen3 library, with my number type being Boost.Multiprecision's mpfr_float wrapper. I can make the matrices just fine, but working with them fails for everything I've tried except matrix addition. Merely multiplication of two identity matrices produces garbage results!

Here's a MWE:

#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/LU>
#include <boost/multiprecision/mpfr.hpp>
#include <iostream>

namespace Eigen{

using boost::multiprecision::mpfr_float;
    template<> struct NumTraits<boost::multiprecision::mpfr_float>
    {

    typedef boost::multiprecision::mpfr_float Real;
    typedef boost::multiprecision::mpfr_float NonInteger;
    typedef boost::multiprecision::mpfr_float Nested;
    enum {
        IsComplex = 0,
        IsInteger = 0,
        IsSigned = 1,
        RequireInitialization = 1,
        ReadCost = 20, //these have no impact
        AddCost = 30,
        MulCost = 40
    };


    inline static Real highest() { // these seem to have no impact
        return (mpfr_float(1) - epsilon()) * pow(mpfr_float(2),mpfr_get_emax()-1);
    }

    inline static Real lowest() {
        return -highest();
    }

    inline static Real dummy_precision(){
        return pow(mpfr_float(10),-int(mpfr_float::default_precision()-3));
    }

    inline static Real epsilon(){
        return pow(mpfr_float(10),-int( mpfr_float::default_precision()));
    }
    //http://www.manpagez.com/info/mpfr/mpfr-2.3.2/mpfr_31.php
};
} // namespace eigen


int main()
{
    int size = 10;

    typedef Eigen::Matrix<boost::multiprecision::mpfr_float, Eigen::Dynamic, Eigen::Dynamic> mp_matrix;
    mp_matrix A = mp_matrix::Identity(size, size);
    std::cout << A * A << std::endl; // produces nan's every other row!!!

    return 0;
}

It produces the identity matrix just fine, but on my machine, using the latest homebrew-distributed versions (and others) of the dependencies for this code (Boost 1.57, Eigen 3.2.4), my program produces NaN's every other row in the matrix:

  1   0   0   0   0   0   0   0   0   0
nan nan nan nan nan nan nan nan nan nan
  0   0   1   0   0   0   0   0   0   0
nan nan nan nan nan nan nan nan nan nan
  0   0   0   0   1   0   0   0   0   0
nan nan nan nan nan nan nan nan nan nan
  0   0   0   0   0   0   1   0   0   0
nan nan nan nan nan nan nan nan nan nan
  0   0   0   0   0   0   0   0   1   0
nan nan nan nan nan nan nan nan nan nan

Odd matrix sizes produce two rows of nan's at the bottom...

This does not seem to depend on the default precision, or the details of the NumTraits struct I define, or even whether I define one. I can inherit from GenericTraits<mpfr_float>, or not; I can say RequireInitialization = 1, or 0. I get NaN's. If I try to LU invert to solve a system, the returned matrix is entirely NaN. If the size of the matrix is 1x1, I even get a single NaN from matrix multiplication. Changing the various static functions has no impact either.

I feel like the strangest part is that if I define a custom complex class (not std::complex, for data loss reasons), with mpfr_float's as the underlying type for real and imaginary parts, I DO get functional matrices.

edit : Here is the complex type's NumTraits:

/**
 \brief this templated struct permits us to use the Float type in Eigen matrices.
 */
template<> struct NumTraits<mynamespace::complex> : NumTraits<boost::multiprecision::mpfr_float> // permits to get the epsilon, dummy_precision, lowest, highest functions
{
    typedef boost::multiprecision::mpfr_float Real;
    typedef boost::multiprecision::mpfr_float NonInteger;
    typedef mynamespace::complex Nested;
    enum {
        IsComplex = 1,
        IsInteger = 0,
        IsSigned = 1,
        RequireInitialization = 1, // yes, require initialization, otherwise get crashes
        ReadCost = 2 * NumTraits<Real>::ReadCost,
        AddCost = 2 * NumTraits<Real>::AddCost,
        MulCost = 4 * NumTraits<Real>::MulCost + 2 * NumTraits<Real>::AddCost
    };
};

Here is the complex class I wrote:

#include <boost/multiprecision/mpfr.hpp>
#include <boost/multiprecision/random.hpp>


#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
#include <boost/serialization/split_member.hpp>

#include <eigen3/Eigen/Core>

#include <assert.h>
namespace mynamespace {
    using boost::multiprecision::mpfr_float;

    class complex {

    private:

        mpfr_float real_, imag_; ///< the real and imaginary parts of the complex number


        // let the boost serialization library have access to the private members of this class.
        friend class boost::serialization::access;

        template<class Archive>
        void save(Archive & ar, const unsigned int version) const {
            // note, version is always the latest when saving
            ar & real_;
            ar & imag_;
        }


        /**
         \brief load method for archiving a bertini::complex
         */
        template<class Archive>
        void load(Archive & ar, const unsigned int version) {
            ar & real_;
            ar & imag_;
        }

        BOOST_SERIALIZATION_SPLIT_MEMBER()


    public:

        complex():real_(), imag_(){}

        complex(double re) : real_(re), imag_("0.0"){}

        complex(const mpfr_float & re) : real_(re), imag_("0.0"){}

        complex(const std::string & re) : real_(re), imag_("0.0"){}

        complex(const mpfr_float & re, const mpfr_float & im) : real_(re), imag_(im) {}

        complex(double re, double im) : real_(re), imag_(im) {}

        complex(const std::string & re, const std::string & im) : real_(re), imag_(im) {}

        complex(const mpfr_float & re, const std::string & im) : real_(re), imag_(im) {}

        complex(const std::string & re, const mpfr_float & im) : real_(re), imag_(im) {}

        complex(complex&& other) : complex() {
            swap(*this, other);
        }

        complex(const complex & other) : real_(other.real_), imag_(other.imag_) {}

        friend void swap(complex& first, complex& second)  {
            using std::swap;
            swap(first.real_,second.real_);
            swap(first.imag_,second.imag_);
        }

        complex& operator=(complex other) {
            swap(*this, other);
            return *this;
        }

        mpfr_float real() const {return real_;}

        mpfr_float imag() const {return imag_;}

        void real(const mpfr_float & new_real){real_ = new_real;}

        void imag(const mpfr_float & new_imag){imag_ = new_imag;}

        void real(const std::string & new_real){real_ = mpfr_float(new_real);}

        void imag(const std::string & new_imag){imag_ = mpfr_float(new_imag);}


        complex& operator+=(const complex & rhs) {
            real_+=rhs.real_;
            imag_+=rhs.imag_;
            return *this;
        }

        complex& operator-=(const complex & rhs) {
            real_-=rhs.real_;
            imag_-=rhs.imag_;
            return *this;
        }

        complex& operator*=(const complex & rhs) {
            mpfr_float a = real_*rhs.real_ - imag_*rhs.imag_; // cache the real part of the result
            imag_ = real_*rhs.imag_ + imag_*rhs.real_;
            real_ = a;
            return *this;
        }

        complex& operator/=(const complex & rhs) {
            mpfr_float d = rhs.abs2();
            mpfr_float a = real_*rhs.real_ + imag_*rhs.imag_; // cache the numerator of the real part of the result
            imag_ = imag_*rhs.real_ - real_*rhs.imag_/d;
            real_ = a/d;

            return *this;
        }

        complex operator-() const 
            return complex(-real(), -imag());
        }

        mpfr_float abs2() const {
            return pow(real(),2)+pow(imag(),2);
        }

        mpfr_float abs() const {
            return sqrt(abs2());
        }

        mpfr_float arg() const {
            return boost::multiprecision::atan2(imag(),real());
        }

        mpfr_float norm() const {
            return abs2();
        }

        complex conj() const {
            return complex(real(), -imag());
        }

        void precision(unsigned int prec) {
            real_.precision(prec);
            imag_.precision(prec);
        }

        unsigned int precision() const {
            assert(real_.precision()==imag_.precision());
            return real_.precision();
        }

        friend std::ostream& operator<<(std::ostream& out, const complex & z) {
            out << "(" << z.real() << "," << z.imag() << ")";
            return out;
        }

        friend std::istream& operator>>(std::istream& in, complex & z) {
            std::string gotten;
            in >> gotten;

            if (gotten[0]=='(') {
                if (*(gotten.end()-1)!=')') {
                    in.setstate(std::ios::failbit);
                    z.real("NaN");
                    z.imag("NaN");
                    return in;
                }
                else{
                    // try to find a comma in the string.
                    size_t comma_pos = gotten.find(",");

                    // if the second character, have no numbers in the real part.
                    // if the second to last character, have no numbers in the imag part.

                    if (comma_pos!=std::string::npos){
                        if (comma_pos==1 || comma_pos==gotten.size()-2) {
                            in.setstate(std::ios::failbit);
                            z.real("NaN");
                            z.imag("NaN");
                            return in;
                        }
                        else{
                            z.real(gotten.substr(1, comma_pos-1));
                            z.imag(gotten.substr(comma_pos+1, gotten.size()-2 - (comma_pos)));
                            return in;
                        }
                    }
                    // did not find a comma
                    else{
                        z.real(gotten.substr(1,gotten.size()-2));
                        z.imag("0.0");
                        return in;
                    }

                }
            }
            else{
                z.real(gotten);
                z.imag("0.0");
                return in;
            }
        }
    }; // end declaration of the mynamespace::complex number class

    inline complex operator+(complex lhs, const complex & rhs){
        lhs += rhs;
        return lhs;
    }

    inline complex operator+(complex lhs, const mpfr_float & rhs) {
        lhs.real(lhs.real()+rhs);
        return lhs;
    }

    inline complex operator+(const mpfr_float & lhs, complex rhs) {
        return rhs+lhs;
    }

    inline complex operator-(complex lhs, const complex & rhs){
        lhs -= rhs;
        return lhs;
    }

    inline complex operator-(complex lhs, const mpfr_float & rhs) {
        lhs.real(lhs.real()-rhs);
        return lhs;
    }

    inline complex operator-(const mpfr_float & lhs, complex rhs) {
        rhs.real(lhs - rhs.real());
        return rhs;
    }

    inline complex operator*(complex lhs, const complex & rhs){
        lhs *= rhs;
        return lhs;
    }

    inline complex operator*(complex lhs, const mpfr_float & rhs) {
        lhs.real(lhs.real()*rhs);
        lhs.imag(lhs.imag()*rhs);
        return lhs;
    }

    inline complex operator*(const mpfr_float & lhs, complex rhs) {
        return rhs*lhs; // it commutes!
    }

    inline complex operator/(complex lhs, const complex & rhs){
        lhs /= rhs;
        return lhs;
    }

    inline complex operator/(complex lhs, const mpfr_float & rhs) {
        lhs.real(lhs.real()/rhs);
        lhs.imag(lhs.imag()/rhs);
        return lhs;
    }

    inline complex operator/(const mpfr_float & lhs, const complex & rhs) {
        mpfr_float d = rhs.abs2();
        return complex(lhs*rhs.real()/d, -lhs*rhs.imag()/d);
    }

    inline mpfr_float real(const complex & z) {
        return z.real();
    }

    inline mpfr_float imag(const complex & z) {
        return z.imag();
    }

    inline complex conj(const complex & z) {
        return z.conj();
    }

    inline mpfr_float abs2(const complex & z) {
        return z.abs2();
    }

    inline mpfr_float abs(const complex & z) {
        return boost::multiprecision::sqrt(abs2(z));
    }

    inline mpfr_float arg(const complex & z) {
        return boost::multiprecision::atan2(z.imag(),z.real());
    }

    inline complex inverse(const complex & z) {
        mpfr_float d = z.abs2();

        return complex(z.real()/d, -z.imag()/d);
    }

    inline complex square(const complex & z) {
        return complex(z.real()*z.real() - z.imag()*z.imag(), mpfr_float("2.0")*z.real()*z.imag());
    }

    inline complex pow(const complex & z, int power) {
        if (power < 0) {
            return pow(inverse(z), -power);
        }
        else if (power==0)
            return complex("1.0","0.0");
        else if(power==1)
            return z;
        else if(power==2)
            return z*z;
        else if(power==3)
            return z*z*z;
        else {
            unsigned int p(power);
            complex result("1.0","0.0"), z_to_the_current_power_of_two = z;
            // have copy of p in memory, can freely modify it.
            do {
                if ( (p & 1) == 1 ) { // get the lowest bit of the number
                    result *= z_to_the_current_power_of_two;
                }
                z_to_the_current_power_of_two *= z_to_the_current_power_of_two; // square z_to_the_current_power_of_two
            } while (p  >>= 1);

            return result;
        }
    }

    inline complex polar(const mpfr_float & rho, const mpfr_float & theta) {
        return complex(rho*cos(theta), rho*sin(theta));
    }

    inline complex sqrt(const complex & z) {
        return polar(sqrt(abs(z)), arg(z)/2);
    }
} // re: namespace

What am I doing wrong? What can I do to Eigen / NumTraits / etc to get matrix operations to work correctly?

Answer 1

This strange nan behavior was caused be a bug in the boost/multiprecision/mpfr.hpp file, appearing in Boost versions 1.56, 1.57 and presumably earlier. The assignment operator did not test for self-assignment -- so the number was set to nan.

Why this happened every other row, and on the last row if the matrix size was odd, is still unclear to me. However, adding a test for self-assignment, à la the commit on GitHub here, fixes the problem. An official patch from the maintainer is forthcoming, but will/did not make it into 1.58. I give thanks to John on the Boost user's mailing list for finding the bug.

If you want to fix this bug in your installation of Boost, in boost/multiprecision/mpfr.hpp, simply wrap the contents (not including return) of the reference assignment operator in if(this != &o){...}.

One lingering question I have regards the Eigen implementation -- what is causing the self-assignment? Is it a problem in general?