I am trying to get the eigen3 package, particularly the JacobiSVD() routines to work, but running into problems.
Here's a sample program, it crashes with a float error (probably divide by zero)

I have followed the information given by the eigen3 tux group, on implementing a new number type and the num traits seems right.

My goal is very high accuracy (100 to 2000 digits) SVD factoring. I redid the source code again here, so that you can see the expected results.

Does anyone know how to fix the float error?

- American Gentile

Code: Select all

#include <gmpxx.h>
#include <iostream>
#include <eigen3/Eigen/Dense>

using namespace std;
using namespace Eigen;

// add in Eigen num traits

namespace Eigen {
 
template<> struct NumTraits<mpq_class> : GenericNumTraits<mpq_class>
{
   typedef mpf_class Real;
   typedef mpf_class NonInteger;
   typedef mpf_class Nested;
   typedef mpf_class Literal;
   static inline Real epsilon() { return 0; }
   static inline Real dummy_precision() { return 0; }
   static inline int digits10() { return 0; }
 
  enum {
    IsComplex = 0,
    IsInteger = 0,
    IsSigned = 1,
    RequireInitialization = 1,
    ReadCost = 1,
    AddCost = 3,
    MulCost = 30
  };
};
// add new typedef for matrix_type
typedef Matrix<mpf_class,Dynamic,Dynamic> matrix_type;
}

int main() {

   mpf_set_default_prec(128); // 38 bits decimal precision

   // these are 8 3-d coordinate points

   mpf_class pt1_x("-0.56656949032236077386013072551434190014");
   mpf_class pt1_y("-0.36279379728181290638975617243847270109");
   mpf_class pt1_z("0.32664219511729088023603923691018758812");

   mpf_class pt2_x("-0.56656949032236077386013072551434190014");
   mpf_class pt2_y("0.36279379728181290638975617243847270109");
   mpf_class pt2_z("0.32664219511729088023603923691018758812");

   mpf_class pt3_x("0.0039438609082687298612700637868822416726");
   mpf_class pt3_y("0.57614383882976695865220682816995171747");
   mpf_class pt3_z("-0.0022737394199980766529093877946320654819");

   mpf_class pt4_x("0.24983837281500445810831904053722344881");
   mpf_class pt4_y("0.38833714139730269358036754420691108744");
   mpf_class pt4_z("-0.14403838525507757655483361372390547624");

   mpf_class pt5_x("0.31278725659908758589054162119023620967");
   mpf_class pt5_y("0.28023506772424802131449095951480123048");
   mpf_class pt5_z("-0.18033007044221522702829623539165004638");

   mpf_class pt6_x("0.31278725659908758589054162119023620967");
   mpf_class pt6_y("-0.28023506772424802131449095951480123048");
   mpf_class pt6_z("-0.18033007044221522702829623539165004638");

   mpf_class pt7_x("0.24983837281500445810831904053722344881");
   mpf_class pt7_y("-0.38833714139730269358036754420691108744");
   mpf_class pt7_z("-0.14403838525507757655483361372390547624");

   mpf_class pt8_x("0.0039438609082687298612700637868822416726");
   mpf_class pt8_y("-0.57614383882976695865220682816995171747");
   mpf_class pt8_z("-0.0022737394199980766529093877946320654819");

   matrix_type M(8,3); // constructor which sets dynamic size = 8 x 3

   // load up the entries in the matrix

   M(0,0) = pt1_x; M(0,1) = pt1_y; M(0,2) = pt1_z;
   M(1,0) = pt2_x; M(1,1) = pt2_y; M(1,2) = pt2_z;
   M(2,0) = pt3_x; M(2,1) = pt3_y; M(2,2) = pt3_z;
   M(3,0) = pt4_x; M(3,1) = pt4_y; M(3,2) = pt4_z;
   M(4,0) = pt5_x; M(4,1) = pt5_y; M(4,2) = pt5_z;
   M(5,0) = pt6_x; M(5,1) = pt6_y; M(5,2) = pt6_z;
   M(6,0) = pt7_x; M(6,1) = pt7_y; M(6,2) = pt7_z;
   M(7,0) = pt8_x; M(7,1) = pt8_y; M(7,2) = pt8_z;

   //cout << "M = " << endl << M << endl;

   // attempt to run SVD decomposition

   // all we want is 1 column V but very accurate values

   cout << "-------------------------------------------" << endl;
   cout << "Output of this program:" << endl;
   cout << "Floating point exception (core dumped)" << endl;
   cout << "" << endl;
   cout << "This program crashes with floating point exception (divide by zero)" << endl;
   cout << "" << endl;
   cout << "running gdb debugger shows:" << endl;
   cout << " at /usr/local/include/eigen3/Eigen/src/misc/RealSvd2x2.h:41" << endl;
   cout << "41       RealScalar u = t / d;" << endl;
   cout << "(gdb) print t" << endl;
   cout << "$3 = {mp = {{_mp_prec = 2, _mp_size = -2, _mp_exp = 1, _mp_d = 0x6bb700}}}" << endl;
   cout << "(gdb) print d" << endl;
   cout << "$4 = {mp = {{_mp_prec = 2, _mp_size = 0, _mp_exp = 0, _mp_d = 0x6ba660}}}" << endl;
   cout << "it is inferred that d = 0, above" << endl;
   cout << "-------------------------------------------" << endl;

   cout << "  V matrix of the SVD decomposition" << endl;
   cout << "[0.49946451054082379982192394858171300293  0.86633434810713543089673327400005492051 0]" << endl;
   cout << "" << endl;
   cout << "[0                                         0  1.0000000000000000000000000000000000000]" << endl;
   cout << "" << endl;
   cout << "[0.86633434810713543089673327400005492051 -0.49946451054082379982192394858171300293 0]" << endl;
   cout << "" << endl;
   cout << "First column of V:" << endl;
   cout << "[0.49946451054082379982192394858171300293, 0, 0.86633434810713543089673327400005492051]" << endl;
   cout << "Magnitude of First Col of V = 1" << endl;

   cout << "" << endl;
   cout << "compile line:" << endl;
   cout << "              g++ -o test-eigen source/test-eigen.cpp -Wall -lgmpxx -lgmp -ggdb" << endl;

   JacobiSVD<matrix_type> svd(M, ComputeThinV);

   cout << svd.matrixV() << endl;

   return(0);
}


I was able to resolve my problem with the C++ code by adding 2 includes on the header file,
#include <eigen3/Eigen/Core>
#include <eigen3/Eigen/SVD>
which then solved my main problem.

As to the code above, I cannot get it to run yet, even though I have changed the includes to match the above.

UPDATE: The compile problem is fixed, but I still cannot get my source code, nor the program above to run, they both get float exception error (divide by zero) which seems to be coming from unitialized values for RealScalar, as apparently they are NOT being updated and converted to mpf_class gmp C++ float class variables.

Any ideas? I tried to carefully follow the mpq_class implementation given by the example on the tuk page, but the code continues to exception crash.

- American code programmer struggling to implement eigen3 with mpf_class floats.