diff --git a/src/Mod/Sketcher/App/planegcs/GCS.cpp b/src/Mod/Sketcher/App/planegcs/GCS.cpp index 12c6f3dcc5..c5093970b9 100644 --- a/src/Mod/Sketcher/App/planegcs/GCS.cpp +++ b/src/Mod/Sketcher/App/planegcs/GCS.cpp @@ -3961,40 +3961,7 @@ SolverReportingManager::Manager().LogToFile("GCS::System::diagnose()\n"); if (J.rows() > 0) { // Get dependent parameters - { - Eigen::MatrixXd Rparams; - Eigen::SparseQR, Eigen::COLAMDOrdering > SqrJ; - - int nontransprank; // will be the same as for the transpose, but better use a different name - - makeSparseQRDecomposition( J, jacobianconstraintmap, SqrJ, nontransprank, Rparams, false); // do not transpose allow to diagnose parameters - - //int constrNum = SqrJ.rows(); // this is the other way around than for the transposed J - //int paramsNum = SqrJ.cols(); - - eliminateNonZerosOverPivotInUpperTriangularMatrix(Rparams, rank); - -#ifdef _GCS_DEBUG - SolverReportingManager::Manager().LogMatrix("Rparams", Rparams); -#endif - //std::vector< std::vector > dependencyGroups(SqrJ.cols()-rank); - for (int j=nontransprank; j < SqrJ.cols(); j++) { - for (int row=0; row < rank; row++) { - if (fabs(Rparams(row,j)) > 1e-10) { - int origCol = SqrJ.colsPermutation().indices()[row]; - - //dependencyGroups[j-rank].push_back(pdiagnoselist[origCol]); - pdependentparameters.push_back(pdiagnoselist[origCol]); - } - } - int origCol = SqrJ.colsPermutation().indices()[j]; - - //dependencyGroups[j-rank].push_back(pdiagnoselist[origCol]); - pdependentparameters.push_back(pdiagnoselist[origCol]); - } - - } - + identifyDependentParametersSparseQR(J, jacobianconstraintmap, pdiagnoselist, true); // Detecting conflicting or redundant constraints if (constrNum > rank) { // conflicting or redundant constraints @@ -4023,7 +3990,8 @@ SolverReportingManager::Manager().LogToFile("GCS::System::diagnose()\n"); return dofs; } -void System::makeDenseQRDecomposition( const Eigen::MatrixXd &J, std::map &jacobianconstraintmap, +void System::makeDenseQRDecomposition( const Eigen::MatrixXd &J, + const std::map &jacobianconstraintmap, Eigen::FullPivHouseholderQR& qrJT, int &rank, Eigen::MatrixXd & R) { @@ -4077,9 +4045,10 @@ void System::makeDenseQRDecomposition( const Eigen::MatrixXd &J, std::map &jacobianconstraintmap, +void System::makeSparseQRDecomposition( const Eigen::MatrixXd &J, + const std::map &jacobianconstraintmap, Eigen::SparseQR, Eigen::COLAMDOrdering > &SqrJT, - int &rank, Eigen::MatrixXd & R, bool transposeJ) + int &rank, Eigen::MatrixXd & R, bool transposeJ, bool silent) { #ifdef EIGEN_SPARSEQR_COMPATIBLE Eigen::SparseMatrix SJ; @@ -4091,7 +4060,8 @@ void System::makeSparseQRDecomposition( const Eigen::MatrixXd &J, std::map 0) { + if (J.rows() > 0 && !silent) { SolverReportingManager::Manager().LogMatrix("R", R); @@ -4161,6 +4130,45 @@ void System::makeSparseQRDecomposition( const Eigen::MatrixXd &J, std::map &jacobianconstraintmap, + const GCS::VEC_pD &pdiagnoselist, + bool silent) +{ + Eigen::MatrixXd Rparams; + Eigen::SparseQR, Eigen::COLAMDOrdering > SqrJ; + + int nontransprank; + + makeSparseQRDecomposition( J, jacobianconstraintmap, SqrJ, nontransprank, Rparams, false, true); // do not transpose allow to diagnose parameters + + //int constrNum = SqrJ.rows(); // this is the other way around than for the transposed J + //int paramsNum = SqrJ.cols(); + + eliminateNonZerosOverPivotInUpperTriangularMatrix(Rparams, nontransprank); + +#ifdef _GCS_DEBUG + if(!silent) + SolverReportingManager::Manager().LogMatrix("Rparams", Rparams); +#endif + //std::vector< std::vector > dependencyGroups(SqrJ.cols()-rank); + for (int j=nontransprank; j < SqrJ.cols(); j++) { + for (int row=0; row < nontransprank; row++) { + if (fabs(Rparams(row,j)) > 1e-10) { + int origCol = SqrJ.colsPermutation().indices()[row]; + + //dependencyGroups[j-rank].push_back(pdiagnoselist[origCol]); + pdependentparameters.push_back(pdiagnoselist[origCol]); + } + } + int origCol = SqrJ.colsPermutation().indices()[j]; + + //dependencyGroups[j-rank].push_back(pdiagnoselist[origCol]); + pdependentparameters.push_back(pdiagnoselist[origCol]); + } +} + + void System::identifyDependentGeometryParametersInTransposedJacobianDenseQRDecomposition( const Eigen::FullPivHouseholderQR& qrJT, const GCS::VEC_pD &pdiagnoselist, diff --git a/src/Mod/Sketcher/App/planegcs/GCS.h b/src/Mod/Sketcher/App/planegcs/GCS.h index a771bd8728..3690899827 100644 --- a/src/Mod/Sketcher/App/planegcs/GCS.h +++ b/src/Mod/Sketcher/App/planegcs/GCS.h @@ -133,13 +133,15 @@ namespace GCS void makeReducedJacobian(Eigen::MatrixXd &J, std::map &jacobianconstraintmap, GCS::VEC_pD &pdiagnoselist, std::map< int , int> &tagmultiplicity); - void makeDenseQRDecomposition( const Eigen::MatrixXd &J, std::map &jacobianconstraintmap, + void makeDenseQRDecomposition( const Eigen::MatrixXd &J, + const std::map &jacobianconstraintmap, Eigen::FullPivHouseholderQR& qrJT, int &rank, Eigen::MatrixXd &R); - void makeSparseQRDecomposition( const Eigen::MatrixXd &J, std::map &jacobianconstraintmap, + void makeSparseQRDecomposition( const Eigen::MatrixXd &J, + const std::map &jacobianconstraintmap, Eigen::SparseQR, Eigen::COLAMDOrdering > &SqrJT, - int &rank, Eigen::MatrixXd &R, bool transposeJ = true); + int &rank, Eigen::MatrixXd &R, bool transposeJ = true, bool silent = false); // This function name is long for a reason: // - Only for DenseQR @@ -163,6 +165,11 @@ namespace GCS void eliminateNonZerosOverPivotInUpperTriangularMatrix(Eigen::MatrixXd &R, int rank); + void identifyDependentParametersSparseQR( const Eigen::MatrixXd &J, + const std::map &jacobianconstraintmap, + const GCS::VEC_pD &pdiagnoselist, + bool silent=true); + #ifdef _GCS_EXTRACT_SOLVER_SUBSYSTEM_ void extractSubsystem(SubSystem *subsys, bool isRedundantsolving); #endif