diff --git a/src/Mod/Sketcher/App/planegcs/GCS.cpp b/src/Mod/Sketcher/App/planegcs/GCS.cpp
index 4507019dcb..4836d0fb7e 100644
--- a/src/Mod/Sketcher/App/planegcs/GCS.cpp
+++ b/src/Mod/Sketcher/App/planegcs/GCS.cpp
@@ -108,184 +108,279 @@
 
 typedef Eigen::FullPivHouseholderQR<Eigen::MatrixXd>::IntDiagSizeVectorType MatrixIndexType;
 
-#ifdef _GCS_DEBUG
-void LogMatrix(std::string str, Eigen::MatrixXd matrix )
-{
-#ifdef _DEBUG_TO_FILE
-    std::ofstream stream;
-    stream.open("GCS_debug.txt", std::ofstream::out | std::ofstream::app);
-#else
-    // Debug code starts
-    std::stringstream stream;
-#endif
-
-    stream << '\n' << " " << str << " =" << '\n';
-    stream << "[" << '\n';
-    stream << matrix << '\n' ;
-    stream << "]" << '\n';
-
-#ifdef _DEBUG_TO_FILE
-    stream.flush();
-    stream.close();
-#else
-    const std::string tmp = stream.str();
-
-    Base::Console().Log(tmp.c_str());
-#endif
-}
-
-void LogMatrix(std::string str, MatrixIndexType matrix )
-{
-    #ifdef _DEBUG_TO_FILE
-    std::ofstream stream;
-    stream.open("GCS_debug.txt", std::ofstream::out | std::ofstream::app);
-    #else
-    // Debug code starts
-    std::stringstream stream;
-    #endif
-
-    stream << '\n' << " " << str << " =" << '\n';
-    stream << "[" << '\n';
-    stream << matrix << '\n' ;
-    stream << "]" << '\n';
-
-    #ifdef _DEBUG_TO_FILE
-    stream.flush();
-    stream.close();
-    #else
-    const std::string tmp = stream.str();
-
-    Base::Console().Log(tmp.c_str());
-    #endif
-}
-#endif
-
-void LogString(std::string str)
-{
-    #ifdef _DEBUG_TO_FILE
-    std::ofstream stream;
-    stream.open("GCS_debug.txt", std::ofstream::out | std::ofstream::app);
-    #else
-    // Debug code starts
-    std::stringstream stream;
-    #endif
-
-    stream << str << std::endl;
-
-    #ifdef _DEBUG_TO_FILE
-    stream.flush();
-    stream.close();
-    #else
-    const std::string tmp = stream.str();
-
-    Base::Console().Log(tmp.c_str());
-    #endif
-}
-
 #ifndef EIGEN_STOCK_FULLPIVLU_COMPUTE
 namespace Eigen {
 
-typedef Matrix<double,-1,-1,0,-1,-1> MatrixdType;
+    typedef Matrix<double,-1,-1,0,-1,-1> MatrixdType;
 
-template<>
-FullPivLU<MatrixdType>& FullPivLU<MatrixdType>::compute(const MatrixdType& matrix)
-{
-  m_isInitialized = true;
-  m_lu = matrix;
-
-  const Index size = matrix.diagonalSize();
-  const Index rows = matrix.rows();
-  const Index cols = matrix.cols();
-
-  // will store the transpositions, before we accumulate them at the end.
-  // can't accumulate on-the-fly because that will be done in reverse order for the rows.
-  m_rowsTranspositions.resize(matrix.rows());
-  m_colsTranspositions.resize(matrix.cols());
-  Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
-
-  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
-  m_maxpivot = RealScalar(0);
-  RealScalar cutoff(0);
-
-  for(Index k = 0; k < size; ++k)
-  {
-    // First, we need to find the pivot.
-
-    // biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
-    Index row_of_biggest_in_corner, col_of_biggest_in_corner;
-    RealScalar biggest_in_corner;
-    biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
-                        .cwiseAbs()
-                        .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
-    row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
-    col_of_biggest_in_corner += k; // need to add k to them.
-
-    // when k==0, biggest_in_corner is the biggest coeff absolute value in the original matrix
-    if(k == 0) cutoff = biggest_in_corner * NumTraits<Scalar>::epsilon();
-
-    // if the pivot (hence the corner) is "zero", terminate to avoid generating nan/inf values.
-    // Notice that using an exact comparison (biggest_in_corner==0) here, as Golub-van Loan do in
-    // their pseudo-code, results in numerical instability! The cutoff here has been validated
-    // by running the unit test 'lu' with many repetitions.
-    if(biggest_in_corner < cutoff)
+    template<>
+    FullPivLU<MatrixdType>& FullPivLU<MatrixdType>::compute(const MatrixdType& matrix)
     {
-      // before exiting, make sure to initialize the still uninitialized transpositions
-      // in a sane state without destroying what we already have.
-      m_nonzero_pivots = k;
-      for(Index i = k; i < size; ++i)
-      {
-        m_rowsTranspositions.coeffRef(i) = i;
-        m_colsTranspositions.coeffRef(i) = i;
-      }
-      break;
+        m_isInitialized = true;
+        m_lu = matrix;
+
+        const Index size = matrix.diagonalSize();
+        const Index rows = matrix.rows();
+        const Index cols = matrix.cols();
+
+        // will store the transpositions, before we accumulate them at the end.
+        // can't accumulate on-the-fly because that will be done in reverse order for the rows.
+        m_rowsTranspositions.resize(matrix.rows());
+        m_colsTranspositions.resize(matrix.cols());
+        Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
+
+        m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
+        m_maxpivot = RealScalar(0);
+        RealScalar cutoff(0);
+
+        for(Index k = 0; k < size; ++k)
+        {
+            // First, we need to find the pivot.
+
+            // biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
+            Index row_of_biggest_in_corner, col_of_biggest_in_corner;
+            RealScalar biggest_in_corner;
+            biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
+            .cwiseAbs()
+            .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
+            row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
+            col_of_biggest_in_corner += k; // need to add k to them.
+
+            // when k==0, biggest_in_corner is the biggest coeff absolute value in the original matrix
+            if(k == 0) cutoff = biggest_in_corner * NumTraits<Scalar>::epsilon();
+
+            // if the pivot (hence the corner) is "zero", terminate to avoid generating nan/inf values.
+            // Notice that using an exact comparison (biggest_in_corner==0) here, as Golub-van Loan do in
+            // their pseudo-code, results in numerical instability! The cutoff here has been validated
+            // by running the unit test 'lu' with many repetitions.
+            if(biggest_in_corner < cutoff)
+            {
+                // before exiting, make sure to initialize the still uninitialized transpositions
+                // in a sane state without destroying what we already have.
+                m_nonzero_pivots = k;
+                for(Index i = k; i < size; ++i)
+                {
+                    m_rowsTranspositions.coeffRef(i) = i;
+                    m_colsTranspositions.coeffRef(i) = i;
+                }
+                break;
+            }
+
+            if(biggest_in_corner > m_maxpivot) m_maxpivot = biggest_in_corner;
+
+            // Now that we've found the pivot, we need to apply the row/col swaps to
+            // bring it to the location (k,k).
+
+            m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
+            m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
+            if(k != row_of_biggest_in_corner) {
+                m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
+                ++number_of_transpositions;
+            }
+            if(k != col_of_biggest_in_corner) {
+                m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
+                ++number_of_transpositions;
+            }
+
+            // Now that the pivot is at the right location, we update the remaining
+            // bottom-right corner by Gaussian elimination.
+
+            if(k<rows-1)
+                m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
+            if(k<size-1)
+                m_lu.block(k+1,k+1,rows-k-1,cols-k-1).noalias() -= m_lu.col(k).tail(rows-k-1) * m_lu.row(k).tail(cols-k-1);
+        }
+
+        // the main loop is over, we still have to accumulate the transpositions to find the
+        // permutations P and Q
+
+        m_p.setIdentity(rows);
+        for(Index k = size-1; k >= 0; --k)
+            m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
+
+        m_q.setIdentity(cols);
+        for(Index k = 0; k < size; ++k)
+            m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
+
+        m_det_pq = (number_of_transpositions%2) ? -1 : 1;
+        return *this;
     }
 
-    if(biggest_in_corner > m_maxpivot) m_maxpivot = biggest_in_corner;
-
-    // Now that we've found the pivot, we need to apply the row/col swaps to
-    // bring it to the location (k,k).
-
-    m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
-    m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
-    if(k != row_of_biggest_in_corner) {
-      m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
-      ++number_of_transpositions;
-    }
-    if(k != col_of_biggest_in_corner) {
-      m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
-      ++number_of_transpositions;
-    }
-
-    // Now that the pivot is at the right location, we update the remaining
-    // bottom-right corner by Gaussian elimination.
-
-    if(k<rows-1)
-      m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
-    if(k<size-1)
-      m_lu.block(k+1,k+1,rows-k-1,cols-k-1).noalias() -= m_lu.col(k).tail(rows-k-1) * m_lu.row(k).tail(cols-k-1);
-  }
-
-  // the main loop is over, we still have to accumulate the transpositions to find the
-  // permutations P and Q
-
-  m_p.setIdentity(rows);
-  for(Index k = size-1; k >= 0; --k)
-    m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
-
-  m_q.setIdentity(cols);
-  for(Index k = 0; k < size; ++k)
-    m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
-
-  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
-  return *this;
-}
-
 } // Eigen
 #endif
 
 namespace GCS
 {
 
+class SolverReportingManager
+{
+public:
+    SolverReportingManager(SolverReportingManager const&) = delete;
+    SolverReportingManager(SolverReportingManager&&) = delete;
+    SolverReportingManager& operator=(SolverReportingManager const&) = delete;
+    SolverReportingManager& operator=(SolverReportingManager &&) = delete;
+
+    static SolverReportingManager& Manager();
+
+    inline void LogString(const std::string& str);
+
+    inline void LogToConsole(const std::string& str);
+
+    inline void LogToFile(const std::string& str);
+
+    void LogQRSystemInformation(const System &system, int paramsNum = 0, int constrNum = 0, int rank = 0);
+
+    void LogMatrix(std::string str, Eigen::MatrixXd matrix);
+    void LogMatrix(std::string str, MatrixIndexType matrix );
+
+private:
+    SolverReportingManager();
+    ~SolverReportingManager();
+
+    inline void initStream();
+    inline void flushStream();
+
+private:
+    #ifdef _DEBUG_TO_FILE
+    std::ofstream stream;
+    #endif
+};
+
+SolverReportingManager::SolverReportingManager()
+{
+    initStream();
+}
+
+SolverReportingManager::~SolverReportingManager()
+{
+    #ifdef _DEBUG_TO_FILE
+    stream.flush();
+    stream.close();
+    #endif
+}
+
+void SolverReportingManager::initStream()
+{
+    #ifdef _DEBUG_TO_FILE
+    if(!stream.is_open()) {
+        stream.open("GCS_debug.txt", std::ofstream::out | std::ofstream::app);
+    }
+    #endif
+}
+
+void SolverReportingManager::flushStream()
+{
+    // Akwardly in some systems flushing does not force the write to the file, requiring a close
+    #ifdef _DEBUG_TO_FILE
+    stream.flush();
+    stream.close();
+    #endif
+}
+
+SolverReportingManager& SolverReportingManager::Manager()
+{
+    static SolverReportingManager theInstance;
+
+    return theInstance;
+}
+
+void SolverReportingManager::LogToConsole(const std::string& str)
+{
+    if(str.size() < Base::Console().BufferSize)
+        Base::Console().Log(str.c_str());
+    else
+        Base::Console().Log("SolverReportingManager - Too long string suppressed");
+}
+
+void SolverReportingManager::LogToFile(const std::string& str)
+{
+    #ifdef _DEBUG_TO_FILE
+    initStream();
+
+    stream << str << std::endl;
+
+    flushStream();
+    #else
+    (void)(str); // silence unused parameter
+    LogToConsole("Debugging to file not enabled!");
+    #endif
+}
+
+void SolverReportingManager::LogString(const std::string& str)
+{
+    LogToConsole(str);
+
+    #ifdef _DEBUG_TO_FILE
+    LogToFile(str);
+    #endif
+}
+
+void SolverReportingManager::LogQRSystemInformation(const System &system, int paramsNum, int constrNum, int rank)
+{
+
+    std::stringstream tempstream;
+
+    tempstream  << (system.qrAlgorithm==EigenSparseQR?"EigenSparseQR":(system.qrAlgorithm==EigenDenseQR?"DenseQR":""));
+
+    if (paramsNum > 0) {
+        tempstream
+        #ifdef EIGEN_SPARSEQR_COMPATIBLE
+        << ", Threads: " << Eigen::nbThreads()
+        #endif
+        #ifdef EIGEN_VECTORIZE
+        << ", Vectorization: On"
+        #endif
+        << ", Pivot Threshold: " << system.qrpivotThreshold
+        << ", Params: " << paramsNum
+        << ", Constr: " << constrNum
+        << ", Rank: "   << rank
+        << std::endl;
+    }
+    else {
+        tempstream
+        #ifdef EIGEN_SPARSEQR_COMPATIBLE
+        << ", Threads: " << Eigen::nbThreads()
+        #endif
+        #ifdef EIGEN_VECTORIZE
+        << ", Vectorization: On"
+        #endif
+        << ", Empty Sketch, nothing to solve"
+        << std::endl;
+    }
+
+    LogString(tempstream.str());
+
+}
+
+
+#ifdef _GCS_DEBUG
+void SolverReportingManager::LogMatrix(std::string str, Eigen::MatrixXd matrix )
+{
+    std::stringstream tempstream;
+
+    tempstream << '\n' << " " << str << " =" << '\n';
+    tempstream << "[" << '\n';
+    tempstream << matrix << '\n' ;
+    tempstream << "]" << '\n';
+
+    LogString(tempstream.str());
+
+}
+
+void SolverReportingManager::LogMatrix(std::string str, MatrixIndexType matrix )
+{
+    std::stringstream tempstream;
+
+    stream << '\n' << " " << str << " =" << '\n';
+    stream << "[" << '\n';
+    stream << matrix << '\n' ;
+    stream << "]" << '\n';
+
+    LogString(tempstream.str());
+}
+#endif
+
+
 typedef boost::adjacency_list <boost::vecS, boost::vecS, boost::undirectedS> Graph;
 
 ///////////////////////////////////////
@@ -3664,6 +3759,45 @@ void System::undoSolution()
     resetToReference();
 }
 
+void System::makeReducedJacobian(Eigen::MatrixXd &J,
+                                 std::map<int,int> &jacobianconstraintmap,
+                                 GCS::VEC_pD &pdiagnoselist,
+                                 std::map< int , int> &tagmultiplicity)
+{
+    // construct specific parameter list for diagonose ignoring driven constraint parameters
+    for (int j=0; j < int(plist.size()); j++) {
+        auto result1 = std::find(std::begin(pdrivenlist), std::end(pdrivenlist), plist[j]);
+
+        if (result1 == std::end(pdrivenlist)) {
+            pdiagnoselist.push_back(plist[j]);
+        }
+    }
+
+
+    J = Eigen::MatrixXd::Zero(clist.size(), pdiagnoselist.size());
+
+    int jacobianconstraintcount=0;
+    int allcount=0;
+    for (std::vector<Constraint *>::iterator constr=clist.begin(); constr != clist.end(); ++constr) {
+        (*constr)->revertParams();
+        ++allcount;
+        if ((*constr)->getTag() >= 0 && (*constr)->isDriving()) {
+            jacobianconstraintcount++;
+            for (int j=0; j < int(pdiagnoselist.size()); j++) {
+                J(jacobianconstraintcount-1,j) = (*constr)->grad(pdiagnoselist[j]);
+            }
+
+            // parallel processing: create tag multiplicity map
+            if(tagmultiplicity.find((*constr)->getTag()) == tagmultiplicity.end())
+                tagmultiplicity[(*constr)->getTag()] = 0;
+            else
+                tagmultiplicity[(*constr)->getTag()]++;
+
+            jacobianconstraintmap[jacobianconstraintcount-1] = allcount-1;
+        }
+    }
+}
+
 int System::diagnose(Algorithm alg)
 {
     // Analyses the constrainess grad of the system and provides feedback
@@ -3681,16 +3815,15 @@ int System::diagnose(Algorithm alg)
         dofs = -1;
         return dofs;
     }
-#ifdef _GCS_DEBUG
+
 #ifdef _DEBUG_TO_FILE
-    std::ofstream stream;
-    stream.open("GCS_debug.txt", std::ofstream::out | std::ofstream::app);
-    stream << "GCS::System::diagnose()" << std::endl;
-    stream.flush();
-    stream.close();
-#endif
+SolverReportingManager::Manager().LogToFile("GCS::System::diagnose()\n");
 #endif
 
+    // Input parameters' lists:
+    // plist            =>  list of all the parameters of the system, e.g. each coordinate of a point
+    // pdrivenlist      =>  list of the parameters that are driven by other parameters (e.g. value of driven constraints)
+
     // When adding an external geometry or a constraint on an external geometry the array 'plist' is empty.
     // So, we must abort here because otherwise we would create an invalid matrix and make the application
     // eventually crash. This fixes issues #0002372/#0002373.
@@ -3704,45 +3837,30 @@ int System::diagnose(Algorithm alg)
     conflictingTags.clear();
     redundantTags.clear();
 
-    // construct specific parameter list for diagonose ignoring driven constraint parameters
-    GCS::VEC_pD pdiagnoselist;
-    for (int j=0; j < int(plist.size()); j++) {
-        auto result1 = std::find(std::begin(pdrivenlist), std::end(pdrivenlist), plist[j]);
+    // This QR diagnosis uses a reduced Jacobian matrix to calculate the rank of the system and identify
+    // conflicting and redundant constraints.
+    //
+    // reduced Jacobian matrix
+    // The Jacobian has been reduced to:
+    // 1. only contain driving constraints, but keep a full size (zero padded).
+    // 2. remove the parameters of the values of driven constraints.
+    Eigen::MatrixXd J;
 
-        if (result1 == std::end(pdrivenlist)) {
-            pdiagnoselist.push_back(plist[j]);
-        }
-    }
-
-    // map tag to a tag multiplicity (the number of solver constraints associated with the same tag)
-    std::map< int , int> tagmultiplicity;
-
-    Eigen::MatrixXd J = Eigen::MatrixXd::Zero(clist.size(), pdiagnoselist.size());
-
-    // The jacobian has been reduced to only contain driving constraints. Identification
-    // of constraint indices from this reduced jacobian requires a mapping.
+    // maps the index of the rows of the reduced jacobian matrix (solver constraints) to
+    // the index those constraints would have in a full size Jacobian matrix
     std::map<int,int> jacobianconstraintmap;
 
-    int jacobianconstraintcount=0;
-    int allcount=0;
-    for (std::vector<Constraint *>::iterator constr=clist.begin(); constr != clist.end(); ++constr) {
-        (*constr)->revertParams();
-        ++allcount;
-        if ((*constr)->getTag() >= 0 && (*constr)->isDriving()) {
-            jacobianconstraintcount++;
-            for (int j=0; j < int(pdiagnoselist.size()); j++) {
-                    J(jacobianconstraintcount-1,j) = (*constr)->grad(pdiagnoselist[j]);
-            }
+    // list of parameters to be diagnosed in this routine (removes value parameters from driven constraints)
+    GCS::VEC_pD pdiagnoselist;
 
-            // parallel processing: create tag multiplicity map
-            if(tagmultiplicity.find((*constr)->getTag()) == tagmultiplicity.end())
-                tagmultiplicity[(*constr)->getTag()] = 0;
-            else
-                tagmultiplicity[(*constr)->getTag()]++;
+    // tag multiplicity gives the number of solver constraints associated with the same tag
+    // A tag generally corresponds to the Sketcher constraint index - There are special tag values, like 0 and -1.
+    std::map< int , int> tagmultiplicity;
 
-            jacobianconstraintmap[jacobianconstraintcount-1] = allcount-1;
-        }
-    }
+
+    makeReducedJacobian(J, jacobianconstraintmap, pdiagnoselist, tagmultiplicity);
+
+    // QR decomposition method selection: SparseQR vs DenseQR
 
 #ifdef EIGEN_SPARSEQR_COMPATIBLE
     Eigen::SparseMatrix<double> SJ;
@@ -3763,8 +3881,10 @@ int System::diagnose(Algorithm alg)
     }
 #endif
 
+
+
 #ifdef _GCS_DEBUG
-    LogMatrix("J",J);
+    SolverReportingManager::Manager().LogMatrix("J",J);
 #endif
 
     Eigen::MatrixXd R;
@@ -3774,6 +3894,7 @@ int System::diagnose(Algorithm alg)
     Eigen::MatrixXd R2; // Intended for a trapezoidal matrix, where R is the top triangular matrix of the R2 trapezoidal matrix
 #endif
 
+
     int paramsNum = 0;
     int constrNum = 0;
     int rank = 0;
@@ -3781,7 +3902,7 @@ int System::diagnose(Algorithm alg)
 
     if(qrAlgorithm==EigenDenseQR){
         if (J.rows() > 0) {
-            qrJT.compute(J.topRows(jacobianconstraintcount).transpose());
+            qrJT.compute(J.topRows(jacobianconstraintmap.size()).transpose());
             //Eigen::MatrixXd Q = qrJT.matrixQ ();
 
             paramsNum = qrJT.rows();
@@ -3804,7 +3925,7 @@ int System::diagnose(Algorithm alg)
 #ifdef EIGEN_SPARSEQR_COMPATIBLE
     else if(qrAlgorithm==EigenSparseQR){
         if (SJ.rows() > 0) {
-            auto SJT = SJ.topRows(jacobianconstraintcount).transpose();
+            auto SJT = SJ.topRows(jacobianconstraintmap.size()).transpose();
             if (SJT.rows() > 0 && SJT.cols() > 0) {
                 SqrJT.compute(SJT);
                 // Do not ask for Q Matrix!!
@@ -3839,53 +3960,22 @@ int System::diagnose(Algorithm alg)
 #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;
-        }
-        else {
-            stream
-#ifdef EIGEN_SPARSEQR_COMPATIBLE
-                    << ", Threads: " << Eigen::nbThreads()
-#endif
-#ifdef EIGEN_VECTORIZE
-                    << ", Vectorization: On"
-#endif
-                    << ", Empty Sketch, nothing to solve";
-        }
-
-        const std::string tmp = stream.str();
-
-        LogString(tmp);
-
+        SolverReportingManager::Manager().LogQRSystemInformation(*this, paramsNum, constrNum, rank);
     }
 
     if (J.rows() > 0) {
 #ifdef _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX
-        LogMatrix("R", R);
+        SolverReportingManager::Manager().LogMatrix("R", R);
 
-        LogMatrix("R2", R2);
+        SolverReportingManager::Manager().LogMatrix("R2", R2);
 
         if(qrAlgorithm == EigenDenseQR){ // There is no rowsTranspositions in SparseQR. obtaining Q is buggy in Eigen for SparseQR
-            LogMatrix("Q", Q);
-            LogMatrix("RowTransp", qrJT.rowsTranspositions());
+            SolverReportingManager::Manager().LogMatrix("Q", Q);
+            SolverReportingManager::Manager().LogMatrix("RowTransp", qrJT.rowsTranspositions());
         }
 #ifdef SPARSE_Q_MATRIX
         else if(qrAlgorithm == EigenSparseQR) {
-            LogMatrix("Q", Q);
+            SolverReportingManager::Manager().LogMatrix("Q", Q);
         }
 #endif
 #endif
@@ -3945,7 +4035,7 @@ int System::diagnose(Algorithm alg)
 #ifdef _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX
         std::string tmp = stream.str();
 
-        LogString(tmp);
+        SolverReportingManager::Manager().LogString(tmp);
 #endif
 
         // If not independent, must be dependent
@@ -3993,7 +4083,7 @@ int System::diagnose(Algorithm alg)
         stream << "]" << std::endl;
 
         tmp = stream.str();
-        LogString(tmp);
+        SolverReportingManager::Manager().LogString(tmp);
 #endif
         for( auto param : depParamCols) {
             pdependentparameters.push_back(pdiagnoselist[param]);
@@ -4057,7 +4147,7 @@ int System::diagnose(Algorithm alg)
         stream << "]" << std::endl;
 
         tmp = stream.str();
-        LogString(tmp);
+        SolverReportingManager::Manager().LogString(tmp);
 #endif
 
             // try to remove the conflicting constraints and solve the
diff --git a/src/Mod/Sketcher/App/planegcs/GCS.h b/src/Mod/Sketcher/App/planegcs/GCS.h
index b6d227d683..1e8bc8cb43 100644
--- a/src/Mod/Sketcher/App/planegcs/GCS.h
+++ b/src/Mod/Sketcher/App/planegcs/GCS.h
@@ -51,18 +51,18 @@ namespace GCS
         LevenbergMarquardt = 1,
         DogLeg = 2
     };
-    
+
     enum DogLegGaussStep {
         FullPivLU = 0,
         LeastNormFullPivLU = 1,
         LeastNormLdlt = 2
     };
-    
+
     enum QRAlgorithm {
         EigenDenseQR = 0,
         EigenSparseQR = 1
     };
-    
+
     enum DebugMode {
         NoDebug = 0,
         Minimal = 1,
@@ -77,7 +77,7 @@ namespace GCS
         VEC_pD plist; // list of the unknown parameters
         VEC_pD pdrivenlist; // list of parameters of driven constraints
         MAP_pD_I pIndex;
-        
+
         VEC_pD pdependentparameters; // list of dependent parameters by the system
 
         std::vector<Constraint *> clist;
@@ -107,6 +107,8 @@ namespace GCS
         int solve_LM(SubSystem *subsys, bool isRedundantsolving=false);
         int solve_DL(SubSystem *subsys, bool isRedundantsolving=false);
 
+        void makeReducedJacobian(Eigen::MatrixXd &J, std::map<int,int> &jacobianconstraintmap, GCS::VEC_pD &pdiagnoselist, std::map< int , int> &tagmultiplicity);
+
         #ifdef _GCS_EXTRACT_SOLVER_SUBSYSTEM_
         void extractSubsystem(SubSystem *subsys, bool isRedundantsolving);
         #endif
@@ -122,18 +124,18 @@ namespace GCS
         double qrpivotThreshold;
         DebugMode debugMode;
         double LM_eps;
-        double LM_eps1;          
+        double LM_eps1;
         double LM_tau;
         double DL_tolg;
-        double DL_tolx;          
+        double DL_tolx;
         double DL_tolf;
         double LM_epsRedundant;
-        double LM_eps1Redundant;          
+        double LM_eps1Redundant;
         double LM_tauRedundant;
         double DL_tolgRedundant;
-        double DL_tolxRedundant;          
-        double DL_tolfRedundant;        
-    
+        double DL_tolxRedundant;
+        double DL_tolfRedundant;
+
     public:
         System();
         /*System(std::vector<Constraint *> clist_);*/
@@ -227,7 +229,7 @@ namespace GCS
                                    double* n1, double* n2,
                                    bool flipn1, bool flipn2,
                                    int tagId, bool driving = true);
-        
+
         // internal alignment constraints
         int addConstraintInternalAlignmentPoint2Ellipse(Ellipse &e, Point &p1, InternalAlignmentType alignmentType, int tagId=0, bool driving = true);
         int addConstraintInternalAlignmentEllipseMajorDiameter(Ellipse &e, Point &p1, Point &p2, int tagId=0, bool driving = true);
@@ -251,7 +253,7 @@ namespace GCS
         // If there's only one, a signed value is returned.
         // Effectively, it calculates the error of a UI constraint
         double calculateConstraintErrorByTag(int tagId);
-        
+
         void rescaleConstraint(int id, double coeff);
 
         void declareUnknowns(VEC_pD &params);