Dear Eigen Developers,

We use Eigen in FreeCAD (freecadweb.org) for the underlying mathematics to solve Sketchs (Geometric Constraint Systems). A while ago (around Eigen 3.0.3), there was a patch submitted in Eigen, this:
https://bitbucket.org/eigen/eigen/src/e ... ew-default

that made the DogLeg algorithm we use to fail to converge in some cases. Part wishing for a future further fix by Eigen Developers, part discouraged by the tedious way to extract this "failure to converge" that we see in FreeCAD into a simple code that we could share with you to show where we find the issue, we opted at that time for providing our own implementation of that method in FreeCAD (that is your Eigen FullPivLU::Compute() method before the patch, what I will call legacy code in the following) instead of using the "unstable" patched-code that has been in Eigen since after. Please, note that I here talk about stability in layman terminology as I am not a mathematician myself. By this I simply mean that in some cases (corner cases), the patched version of the function made DogLeg fail to converge (or converge extremely slowly as to cause our users to actually kill the application and lose any not-yet-saved work).

With the advent of Eigen 3.3, and the change in the declaration of FullPivLU::Compute(), we have been doing some work to support Eigen 3.3 and finally I have decided it is about time to provide you with the information we have regarding this issue, mostly, so that if there is a bug, I may provide you with a corner case showing it. If it is a bug, you may decide to solve it, if not I would personally would like to know what is happening. Somehow, I would like to note that this "issue" does not affect FreeCAD in the sense that we can still continue applying the same solution as before (provide our own version of the function in FreeCAD, which as I said, is your version before that patch). I provide you this information for the sake of cooperation between open source projects and with the aim of improving both. So, enough writting, to the topic:

Git Repository

I have created this Git repository:
https://github.com/abdullahtahiriyo/Eig ... onvergence

It provides a simple project to investigate the Eigen issue in the corner case. This repository uses partly code from the FreeCAD project, but without requiring other dependency than Eigen. It takes under 30 seconds to compile.

What is this exactly?

This example demonstrates how a given system (aka blegh example) will cause our DogLeg algorithm to fail to converge for Eigen versions after this commit (around Eigen 3.0.3):
https://bitbucket.org/eigen/eigen/src/e ... ew-default

This failure to converge was reported here:
http://forum.freecadweb.org/viewtopic.p ... 1&start=40

How to use this?

This test was performed against the Eigen 3.3 of 5/11/2015 (7832:ef9fac3eb614), but should be reproducible with any older version (most probably also newer).

1. Clone the repo
2. Open main.cpp and configure the macro EIGEN_3_3_OR_BEYOND:
#define EIGEN_3_3_OR_BEYOND // if you are using eigen 3.3 or higher this must be set, if you are using eigen 3.2 this must be unset.
3. Comment the macro OLD_FULLPIVLU:
//#define OLD_FULLPIVLU // If enabled, this enforces old legacy FullPivLU::compute method
4. Compile the project:
cmake ..
make
5. Run the executable:
./myproject

You will get this output:

Code: Select all

    This is a demontration of an Eigen algorithm relevant to FreeCAD
   
    Populating values for the relevant Subsystem
    Creating the Subsystem
    Solving using DogLeg
    DL: tolg: 1e-80, tolx: 1e-80, tolf: 1e-10, convergence: 1e-10, xsize: 172, csize: 171, maxIter: 100
    DL, Iteration: 0, fx_inf(tolf): 0.000120285, g_inf(tolg): 0.000210914, delta(f(tolx)): 0.3, err(divergingLim): 5.01991e-08
    DL, Iteration: 1, fx_inf(tolf): 0.000120285, g_inf(tolg): 0.000210914, delta(f(tolx)): 0.15, err(divergingLim): 5.01991e-08
    DL, Iteration: 2, fx_inf(tolf): 0.000114152, g_inf(tolg): 0.000100878, delta(f(tolx)): 0.15, err(divergingLim): 3.37895e-08
    DL, Iteration: 3, fx_inf(tolf): 0.000112493, g_inf(tolg): 8.36503e-05, delta(f(tolx)): 0.15, err(divergingLim): 2.71141e-08
    ...
    DL, Iteration: 99, fx_inf(tolf): 5.01622e-05, g_inf(tolg): 1.73602e-05, delta(f(tolx)): 0.0421875, err(divergingLim): 9.25661e-09
    DL: stopcode: 4, Failed


Note that DogLeg is limited to 100 iterations, and that in 100 iterations It did not converge.

5. Compile the project with OLD_FULLPIVLU defined (uncommented) and run the executable. The output is:

Code: Select all

    This is a demontration of an Eigen algorithm relevant to FreeCAD
   
    Populating values for the relevant Subsystem
    Creating the Subsystem
    Solving using DogLeg
    DL: tolg: 1e-80, tolx: 1e-80, tolf: 1e-10, convergence: 1e-10, xsize: 172, csize: 171, maxIter: 100
    DL, Iteration: 0, fx_inf(tolf): 9.96423e-06, g_inf(tolg): 9.96426e-06, delta(f(tolx)): 0.1, err(divergingLim): 4.96434e-11
    DL, Iteration: 1, fx_inf(tolf): 2.49106e-06, g_inf(tolg): 2.49106e-06, delta(f(tolx)): 0.1, err(divergingLim): 3.10271e-12
    DL, Iteration: 2, fx_inf(tolf): 6.22764e-07, g_inf(tolg): 6.22765e-07, delta(f(tolx)): 0.1, err(divergingLim): 1.93919e-13
    DL, Iteration: 3, fx_inf(tolf): 1.55691e-07, g_inf(tolg): 1.55691e-07, delta(f(tolx)): 0.1, err(divergingLim): 1.21199e-14
    DL, Iteration: 4, fx_inf(tolf): 3.89227e-08, g_inf(tolg): 3.89228e-08, delta(f(tolx)): 0.1, err(divergingLim): 7.57496e-16
    DL, Iteration: 5, fx_inf(tolf): 9.73068e-09, g_inf(tolg): 9.73069e-09, delta(f(tolx)): 0.1, err(divergingLim): 4.73436e-17
    DL, Iteration: 6, fx_inf(tolf): 2.43267e-09, g_inf(tolg): 2.43268e-09, delta(f(tolx)): 0.1, err(divergingLim): 2.95903e-18
    DL, Iteration: 7, fx_inf(tolf): 6.08168e-10, g_inf(tolg): 6.08168e-10, delta(f(tolx)): 0.1, err(divergingLim): 1.8498e-19
    DL, Iteration: 8, fx_inf(tolf): 1.52042e-10, g_inf(tolg): 1.52049e-10, delta(f(tolx)): 0.1, err(divergingLim): 1.16563e-20
    DL, Iteration: 9, fx_inf(tolf): 3.80105e-11, g_inf(tolg): 3.80105e-11, delta(f(tolx)): 0.1, err(divergingLim): 7.88338e-22
    DL: stopcode: 1, Success


The algorithm converged in just 10 iterations.

What to look for ?

Most of the code in the repository is just to build this specific system of equations, as the original was a FreeCAD sketch. Looking at main(), the only relevant bit to the convergence is the last but one line:

solve_DL(mysub, false);

This solve_DL() function is defined in main.cpp for convenience and applies the DogLeg algorithm. The only difference between applying the legacy code or the after-Eigen-patch code is in:

h_gn = Jx.fullPivLu().solve(-fx);

This function internally (within Eigen) calls the compute() method, either the legacy one or the patched one, depending on whether the OLD_FULLPIVLU macro is set or not.

Further information

I have a rather limited knowledge of numerical methods. I make you aware of this comment in the event that what is said there may be relevant for solving the issue:
http://forum.freecadweb.org/viewtopic.p ... =40#p36761

My last word is in the form of thanks for this amazing library.

Regards,
Abdullah

Hi,
it would have been better to report the issue right at the time you found it! Anyway, I tried your self-contained example, and it works fine using either 3.1, 3.2 or 3.3, with or without OLD_FULLPIVLU. I always get the same output:

Code: Select all

Creating the Subsystem
Solving using DogLeg
DL: tolg: 1e-80, tolx: 1e-80, tolf: 1e-10, convergence: 1e-10, xsize: 172, csize: 171, maxIter: 100
DL, Iteration: 0, fx_inf(tolf): 9.96427e-06, g_inf(tolg): 9.96431e-06, delta(f(tolx)): 0.1, err(divergingLim): 4.96438e-11
DL, Iteration: 1, fx_inf(tolf): 2.49107e-06, g_inf(tolg): 2.49107e-06, delta(f(tolx)): 0.1, err(divergingLim): 3.10274e-12
DL, Iteration: 2, fx_inf(tolf): 6.22767e-07, g_inf(tolg): 6.22768e-07, delta(f(tolx)): 0.1, err(divergingLim): 1.93921e-13
DL, Iteration: 3, fx_inf(tolf): 1.55692e-07, g_inf(tolg): 1.55692e-07, delta(f(tolx)): 0.1, err(divergingLim): 1.21201e-14
DL, Iteration: 4, fx_inf(tolf): 3.89229e-08, g_inf(tolg): 3.8923e-08, delta(f(tolx)): 0.1, err(divergingLim): 7.57503e-16
DL, Iteration: 5, fx_inf(tolf): 9.73073e-09, g_inf(tolg): 9.73074e-09, delta(f(tolx)): 0.1, err(divergingLim): 4.7344e-17
DL, Iteration: 6, fx_inf(tolf): 2.43268e-09, g_inf(tolg): 2.43269e-09, delta(f(tolx)): 0.1, err(divergingLim): 2.95906e-18
DL, Iteration: 7, fx_inf(tolf): 6.08171e-10, g_inf(tolg): 6.08171e-10, delta(f(tolx)): 0.1, err(divergingLim): 1.84981e-19
DL, Iteration: 8, fx_inf(tolf): 1.52042e-10, g_inf(tolg): 1.52045e-10, delta(f(tolx)): 0.1, err(divergingLim): 1.16211e-20
DL, Iteration: 9, fx_inf(tolf): 3.80105e-11, g_inf(tolg): 3.80105e-11, delta(f(tolx)): 0.1, err(divergingLim): 7.65869e-22
DL: stopcode: 1, Success


(I'm using clang on osx, 64bits, SSE2)

I also tried with the i387 FPU. Again, OLD_FULLPIVLU has no effect, and I get:

Code: Select all

Creating the Subsystem
Solving using DogLeg
DL: tolg: 1e-80, tolx: 1e-80, tolf: 1e-10, convergence: 1e-10, xsize: 172, csize: 171, maxIter: 100
DL, Iteration: 0, fx_inf(tolf): 1.81471e-09, g_inf(tolg): 1.70823e-09, delta(f(tolx)): 0.1, err(divergingLim): 1.83212e-18
DL, Iteration: 1, fx_inf(tolf): 4.79737e-10, g_inf(tolg): 4.79943e-10, delta(f(tolx)): 0.1, err(divergingLim): 1.18768e-19
DL, Iteration: 2, fx_inf(tolf): 1.19934e-10, g_inf(tolg): 1.19978e-10, delta(f(tolx)): 0.1, err(divergingLim): 7.48942e-21
DL, Iteration: 3, fx_inf(tolf): 2.99849e-11, g_inf(tolg): 2.99938e-11, delta(f(tolx)): 0.1, err(divergingLim): 5.05654e-22
DL: stopcode: 1, Success


ok, after looking a bit more are your matrices, they appear to define an under-constraint problem by construction (less equations than unknowns), and some given equations are even linearly dependent. In this case using FullPivLU::solve is not recommended because it will pick an arbitrary solution which might not be the one you are looking for. In this case, unless you have additional knowledge to insert to your problem, a default choice consist in computing the least-norm solution. You can compute it using either:

h_gn = Jx.bdcSvd(Eigen::ComputeFullU|Eigen::ComputeFullV).solve(-fx);

or x3 faster:

h_gn = Jx.adjoint()*(Jx*Jx.adjoint()).fullPivLu().solve(-fx);

or even x2 faster than the previous one:

h_gn = Jx.adjoint()*(Jx*Jx.adjoint()).ldlt().solve(-fx);

Hi Gael,

First I apologise for pursuing it so late. I will try to report as soon as possible in the future.

Second thank you very much for your trials and advice.

I have just re-tried the example against these:
eigen-eigen-b30b87236a1b-3.2.7
eigen-eigen-ca142d0540d3-3.1.0
eigen-eigen-ffa86ffb5570-3.2.0

in addition to Eigen 3.3 of 5/11/2015 (7832:ef9fac3eb614).

I systematically get the failure to converge with all of them using the patched code. I run Ubuntu 14.04 with this configuration:
http://forum.freecadweb.org/viewtopic.p ... 93#p106193

As you see I have asked other users of the FreeCAD project to do some testing to see if we can solve it without stealing more of your valuable time. I will report here back with our findings.

ok, after looking a bit more are your matrices, they appear to define an under-constraint problem by construction (less equations than unknowns), and some given equations are even linearly dependent. In this case using FullPivLU::solve is not recommended because it will pick an arbitrary solution which might not be the one you are looking for. In this case, unless you have additional knowledge to insert to your problem, a default choice consist in computing the least-norm solution

Yes, it is true. Some of the geometric constraints are partially redundant by definition (there is little we can do about this). The user may additionally introduce fully redundant geometric constraints (and even conflicting ones). For the second, we do a rank revealing QR decomposition to know if we have redundant/constraints (rank and constraint number comparison). If we do have them, we remove the constraint from the system and solve it. If the solver success, the constraint is marked as redundant and the user is notified that he has to remove one of the redundant constraints. If the solver fails, the constraint is marked as conflicting and the user is notified that he has to remove one of the conflicting constraints. The example I have self-contained, is one of this solvings (redundant solving) to detect whether redundant or conflicting. We do it like this:

Code: Select all

    if(qrAlgorithm==EigenDenseQR){
        if (J.rows() > 0) {
            qrJT=Eigen::FullPivHouseholderQR<Eigen::MatrixXd>(J.topRows(count).transpose());
            Eigen::MatrixXd Q = qrJT.matrixQ ();
           
            paramsNum = qrJT.rows();
            constrNum = qrJT.cols();
            qrJT.setThreshold(qrpivotThreshold);
            rank = qrJT.rank();

            if (constrNum >= paramsNum)
                R = (qrJT.matrixQR().triangularView<Eigen::Upper>());
            else
                R = qrJT.matrixQR().topRows(constrNum)
                                .triangularView<Eigen::Upper>();       
        }
    }
    #ifdef EIGEN_SPARSEQR_COMPATIBLE
    else if(qrAlgorithm==EigenSparseQR){
        if (SJ.rows() > 0) {
            SqrJT.compute(SJ.topRows(count).transpose());
           
            // Do not ask for Q Matrix!!
            // At Eigen 3.2 still has a bug that this only works for square matrices
            // if enabled it will crash
            //Eigen::SparseMatrix<double> Q = qrJT.matrixQ();
            //qrJT.matrixQ().evalTo(Q);
           
            paramsNum = SqrJT.rows();
            constrNum = SqrJT.cols();
            SqrJT.setPivotThreshold(qrpivotThreshold);
            rank = SqrJT.rank();
           
            if (constrNum >= paramsNum)
                R = SqrJT.matrixR().triangularView<Eigen::Upper>();
            else
                R = SqrJT.matrixR().topRows(constrNum)
                                    .triangularView<Eigen::Upper>();   
        }
    }
    #endif
   
    if(debugMode==IterationLevel) {
        std::stringstream stream;
        stream  << (qrAlgorithm==EigenSparseQR?"EigenSparseQR":(qrAlgorithm==EigenDenseQR?"DenseQR":""));       
       
        if (J.rows() > 0) {
            stream
            #ifdef EIGEN_SPARSEQR_COMPATIBLE
                    << ", Threads: " << Eigen::nbThreads()
            #endif
                    #ifdef EIGEN_VECTORIZE
                    << ", Vectorization: On"
                    #endif           
                    << ", Pivot Threshold: " << qrpivotThreshold                   
                    << ", Params: " << paramsNum
                    << ", Constr: " << constrNum
                    << ", Rank: "   << rank         << "\n";
        }
        else {
            stream 
            #ifdef EIGEN_SPARSEQR_COMPATIBLE
                    << ", Threads: " << Eigen::nbThreads()
            #endif
                    #ifdef EIGEN_VECTORIZE
                    << ", Vectorization: On"
                    #endif               
                    << ", Empty Sketch, nothing to solve" << "\n";           
        }
       
        const std::string tmp = stream.str();
        Base::Console().Log(tmp.c_str());       
    }
       
    if (J.rows() > 0) {
       
        #ifdef _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX
        // Debug code starts
        std::stringstream stream;
       
        stream << "[";
        stream << R ;
        stream << "]";
       
        const std::string tmp = stream.str();
       
        Base::Console().Log(tmp.c_str());
        // Debug code ends
        #endif

        if (constrNum > rank) { // conflicting or redundant constraints
            for (int i=1; i < rank; i++) {
                // eliminate non zeros above pivot
                assert(R(i,i) != 0);
                for (int row=0; row < i; row++) {
                    if (R(row,i) != 0) {
                        double coef=R(row,i)/R(i,i);
                        R.block(row,i+1,1,constrNum-i-1) -= coef * R.block(i,i+1,1,constrNum-i-1);
                        R(row,i) = 0;
                    }
                }
            }
            std::vector< std::vector<Constraint *> > conflictGroups(constrNum-rank);
            for (int j=rank; j < constrNum; j++) {
                for (int row=0; row < rank; row++) {
                    if (fabs(R(row,j)) > 1e-10) {
                        int origCol = 0;
                       
                        if(qrAlgorithm==EigenDenseQR)
                            origCol=qrJT.colsPermutation().indices()[row];
                        #ifdef EIGEN_SPARSEQR_COMPATIBLE
                        else if(qrAlgorithm==EigenSparseQR)
                            origCol=SqrJT.colsPermutation().indices()[row];
                        #endif
                       
                        conflictGroups[j-rank].push_back(clist[origCol]);
                    }
                }
                int origCol = 0;
                       
                if(qrAlgorithm==EigenDenseQR)
                    origCol=qrJT.colsPermutation().indices()[j];
               
                #ifdef EIGEN_SPARSEQR_COMPATIBLE
                else if(qrAlgorithm==EigenSparseQR)
                    origCol=SqrJT.colsPermutation().indices()[j];
                #endif
               
                conflictGroups[j-rank].push_back(clist[origCol]);
            }

            // try to remove the conflicting constraints and solve the
            // system in order to check if the removed constraints were
            // just redundant but not really conflicting
            std::set<Constraint *> skipped;
            SET_I satisfiedGroups;
            while (1) {
                std::map< Constraint *, SET_I > conflictingMap;
                for (std::size_t i=0; i < conflictGroups.size(); i++) {
                    if (satisfiedGroups.count(i) == 0) {
                        for (std::size_t j=0; j < conflictGroups[i].size(); j++) {
                            Constraint *constr = conflictGroups[i][j];
                            if (constr->getTag() != 0) // exclude constraints tagged with zero
                                conflictingMap[constr].insert(i);
                        }
                    }
                }
                if (conflictingMap.empty())
                    break;
           
                int maxPopularity = 0;
                Constraint *mostPopular = NULL;
                for (std::map< Constraint *, SET_I >::const_iterator it=conflictingMap.begin();
                     it != conflictingMap.end(); ++it) {
                    if (static_cast<int>(it->second.size()) > maxPopularity ||
                        (static_cast<int>(it->second.size()) == maxPopularity && mostPopular &&
                         it->first->getTag() > mostPopular->getTag())) {
                        mostPopular = it->first;
                        maxPopularity = it->second.size();
                    }
                }
                if (maxPopularity > 0) {
                    skipped.insert(mostPopular);
                    for (SET_I::const_iterator it=conflictingMap[mostPopular].begin();
                         it != conflictingMap[mostPopular].end(); ++it)
                        satisfiedGroups.insert(*it);
                }
            }

            std::vector<Constraint *> clistTmp;
            clistTmp.reserve(clist.size());
            for (std::vector<Constraint *>::iterator constr=clist.begin();
                 constr != clist.end(); ++constr)
                if (skipped.count(*constr) == 0)
                    clistTmp.push_back(*constr);

            SubSystem *subSysTmp = new SubSystem(clistTmp, plist);
            int res = solve(subSysTmp,true,alg,true);
           
            if(debugMode==Minimal || debugMode==IterationLevel) {
                std::string solvername;
                switch (alg) {
                    case 0:
                        solvername = "BFGS";
                        break;
                    case 1: // solving with the LevenbergMarquardt solver
                        solvername = "LevenbergMarquardt";
                        break;
                    case 2: // solving with the BFGS solver
                        solvername = "DogLeg";
                        break;
                }   
               
            Base::Console().Log("Sketcher::RedundantSolving-%s-\n",solvername.c_str());
                       
            }
           
            if (res == Success) {
                subSysTmp->applySolution();
                for (std::set<Constraint *>::const_iterator constr=skipped.begin();
                     constr != skipped.end(); ++constr) {
                    double err = (*constr)->error();
                    if (err * err < convergenceRedundant)
                        redundant.insert(*constr);
                }
                resetToReference();
               
                if(debugMode==Minimal || debugMode==IterationLevel) {                   
                    Base::Console().Log("Sketcher Redundant solving: %d redundants\n",redundant.size());         
                }

                std::vector< std::vector<Constraint *> > conflictGroupsOrig=conflictGroups;
                conflictGroups.clear();
                for (int i=conflictGroupsOrig.size()-1; i >= 0; i--) {
                    bool isRedundant = false;
                    for (std::size_t j=0; j < conflictGroupsOrig[i].size(); j++) {
                        if (redundant.count(conflictGroupsOrig[i][j]) > 0) {
                            isRedundant = true;
                            break;
                        }
                    }
                    if (!isRedundant)
                        conflictGroups.push_back(conflictGroupsOrig[i]);
                    else
                        constrNum--;
                }
            }
            delete subSysTmp;

            // simplified output of conflicting tags
            SET_I conflictingTagsSet;
            for (std::size_t i=0; i < conflictGroups.size(); i++) {
                for (std::size_t j=0; j < conflictGroups[i].size(); j++) {
                    conflictingTagsSet.insert(conflictGroups[i][j]->getTag());
                }
            }
            conflictingTagsSet.erase(0); // exclude constraints tagged with zero
            conflictingTags.resize(conflictingTagsSet.size());
            std::copy(conflictingTagsSet.begin(), conflictingTagsSet.end(),
                      conflictingTags.begin());

            // output of redundant tags
            SET_I redundantTagsSet;
            for (std::set<Constraint *>::iterator constr=redundant.begin();
                 constr != redundant.end(); ++constr)
                redundantTagsSet.insert((*constr)->getTag());
            // remove tags represented at least in one non-redundant constraint
            for (std::vector<Constraint *>::iterator constr=clist.begin();
                 constr != clist.end(); ++constr)
                if (redundant.count(*constr) == 0)
                    redundantTagsSet.erase((*constr)->getTag());
            redundantTags.resize(redundantTagsSet.size());
            std::copy(redundantTagsSet.begin(), redundantTagsSet.end(),
                      redundantTags.begin());

            if (paramsNum == rank && constrNum > rank) { // over-constrained
                hasDiagnosis = true;
                dofs = paramsNum - constrNum;
                return dofs;
            }
        }

        hasDiagnosis = true;
        dofs = paramsNum - rank;
        return dofs;
    }
    hasDiagnosis = true;
    dofs = plist.size();
    return dofs;


ggael wrote:You can compute it using either:

h_gn = Jx.bdcSvd(Eigen::ComputeFullU|Eigen::ComputeFullV).solve(-fx);

or x3 faster:

h_gn = Jx.adjoint()*(Jx*Jx.adjoint()).fullPivLu().solve(-fx);

or even x2 faster than the previous one:

h_gn = Jx.adjoint()*(Jx*Jx.adjoint()).ldlt().solve(-fx);

Oh thanks for this!! I am still digesting it (and will take me a while). The first option is unsupported Eigen code and I have not tried it yet, but the second converges in my self-contained example in 10 iterations and the third converges in 7 iterations.

Do not lose more time with this. I will come back to you when more people have confirmed that my results can be reproduced in their computers (or can not).

Thank you very much!!!!!
Abdullah

ok, actually I can reproduce on a Linux system (I don't why this make such a difference - does your test rely on random numbers?)

Anyway, there was an issue in solve() which did not used the rank() method but nonzeroPivots(). Fixed in devel branch: https://bitbucket.org/eigen/eigen/commits/616017adf663/

I'll wait a bit before backporting.

ggael wrote:ok, actually I can reproduce on a Linux system (I don't why this make such a difference - does your test rely on random numbers?)

Short answer: Not to my knowledge.

Long answer: I have written some code in FreeCAD to automatically generate the c++ code that would generate an identical system to the one being solved. I took a lot of care to make sure that all the information was transferred. The behaviour between FreeCAD and the self-contained example is the same. The self-contained example is tiny in size compared to FreeCAD, so a random operation that would potentially introduce a random number into the system (that does not exist to my knowledge) in the code of FreeCAD that is not embedded into the self-contained example, would nevertheless be deterministic for the self-contained example (same data is always loaded). I really cannot think of any part of the self-contained example that introduces random data. However, I would lie if I say I know every single line in the self-contained example (code borrowed from FreeCAD).

ggael wrote:Anyway, there was an issue in solve() which did not used the rank() method but nonzeroPivots(). Fixed in devel branch: https://bitbucket.org/eigen/eigen/commits/616017adf663/

I'll wait a bit before backporting.

I can confirm that when linked against Eigen (7907:616017adf663), FreeCAD converges in a few iterations in Ubuntu 14.04.3 LTS using gcc 4.8 with Eigens own implementation (i.e. with your fix).

Other FreeCAD developers and testers have found some problems when compiling and executing the example in different platforms:
http://forum.freecadweb.org/viewtopic.p ... 41#p106251

I would like to make sure we identify any further existing issue, as it seems that some parts, either of my project or eigen related parts might have compiling issues, for example when compiling with MSVC 64 bits (Windows), and also when executing the Windows builds.

I have also asked for further information and tests with the fixed version under other platforms. I will wait for their cooperation (I use neither Windows nor OSX) and report back here.

Thank you so much for taking the time to look into this.
Abdullah

I've backported the change in FullPivLU::solve, and older MSVC builds are fixed too.

Hi Gael,

The tests performed so far are positive. We will switch back to the mainstream Eigen implementation in the near future, when Debian updates their buggy Eigen-3.3-alpha1 package (currently it is just a distribution issue, not a bug anymore, or so we think). We will let you know if anything else arises.

Thank you very much as always for your prompt support and avid bug fixing

Regards,
Abdullah