Skip to content

Commit

Permalink
GCS: Asynchronous full geometry parameter identification for Dense QR
Browse files Browse the repository at this point in the history
  • Loading branch information
@abdullahtahiriyo
abdullahtahiriyo committed Dec 19, 2020
1 parent 5e2a2cb commit b2fd91c
Show file tree
Hide file tree
Showing 2 changed files with 97 additions and 37 deletions.
120 changes: 84 additions & 36 deletions src/Mod/Sketcher/App/planegcs/GCS.cpp
View file
Expand Up @@ -3922,20 +3922,27 @@ SolverReportingManager::Manager().LogToFile("GCS::System::diagnose()\n");
#ifdef PROFILE_DIAGNOSE
Base::TimeInfo DenseQR_start_time;
#endif
int rank = 0; // rank is not cheap to retrieve from qrJT in DenseQR
Eigen::MatrixXd R;
Eigen::FullPivHouseholderQR<Eigen::MatrixXd> qrJT;
if (J.rows() > 0) {
int rank = 0; // rank is not cheap to retrieve from qrJT in DenseQR
Eigen::MatrixXd R;
Eigen::FullPivHouseholderQR<Eigen::MatrixXd> qrJT;
// Here we enforce calculating the two QR decompositions in parallel
// Care to wait() for the future before any prospective detection of conflicting/redundant, because the redundant solve
// modifies pdiagnoselist and it would not be thread-safe
//
// identifyDependentParametersSparseQR(J, jacobianconstraintmap, pdiagnoselist, true)
auto fut = std::async(std::launch::async,&System::identifyDependentParametersDenseQR,this,J,jacobianconstraintmap, pdiagnoselist, true);

makeDenseQRDecomposition( J, jacobianconstraintmap, qrJT, rank, R);
makeDenseQRDecomposition( J, jacobianconstraintmap, qrJT, rank, R);

int paramsNum = qrJT.rows();
int constrNum = qrJT.cols();
int paramsNum = qrJT.rows();
int constrNum = qrJT.cols();

if (J.rows() > 0) {
identifyDependentGeometryParametersInTransposedJacobianDenseQRDecomposition(
qrJT,
pdiagnoselist,
paramsNum, rank);
// This function is legacy code that was used to obtain partial geometry dependency information from a SINGLE Dense QR
// decomposition. I am reluctant to remove it from here until everything new is well tested.
//identifyDependentGeometryParametersInTransposedJacobianDenseQRDecomposition( qrJT, pdiagnoselist, paramsNum, rank);

fut.wait(); // wait for the execution of identifyDependentParametersSparseQR to finish

dofs = paramsNum - rank; // unless overconstraint, which will be overridden below

Expand Down Expand Up @@ -4014,48 +4021,63 @@ SolverReportingManager::Manager().LogToFile("GCS::System::diagnose()\n");
void System::makeDenseQRDecomposition( const Eigen::MatrixXd &J,
const std::map<int,int> &jacobianconstraintmap,
Eigen::FullPivHouseholderQR<Eigen::MatrixXd>& qrJT,
int &rank, Eigen::MatrixXd & R)
int &rank, Eigen::MatrixXd & R, bool transposeJ, bool silent)
{

#ifdef _GCS_DEBUG
SolverReportingManager::Manager().LogMatrix("J",J);
if(!silent)
SolverReportingManager::Manager().LogMatrix("J",J);
#endif

#ifdef _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX
Eigen::MatrixXd Q; // Obtaining the Q matrix with Sparse QR is buggy, see comments below
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;
// For a transposed J SJG rows are paramsNum and cols are constrNum
// For a non-transposed J SJG rows are constrNum and cols are paramsNum
int rowsNum = 0;
int colsNum = 0;

if (J.rows() > 0) {
qrJT.compute(J.topRows(jacobianconstraintmap.size()).transpose());
//Eigen::MatrixXd Q = qrJT.matrixQ ();
Eigen::MatrixXd JG;
if(transposeJ)
JG = J.topRows(jacobianconstraintmap.size()).transpose();
else
JG = J.topRows(jacobianconstraintmap.size());

paramsNum = qrJT.rows();
constrNum = qrJT.cols();
qrJT.setThreshold(qrpivotThreshold);
rank = qrJT.rank();
if (JG.rows() > 0 && JG.cols() > 0) {

if (constrNum >= paramsNum)
R = qrJT.matrixQR().triangularView<Eigen::Upper>();
else
R = qrJT.matrixQR().topRows(constrNum)
.triangularView<Eigen::Upper>();
qrJT.compute(JG);

rowsNum = qrJT.rows();
colsNum = qrJT.cols();
qrJT.setThreshold(qrpivotThreshold);
rank = qrJT.rank();

if (colsNum >= rowsNum)
R = qrJT.matrixQR().triangularView<Eigen::Upper>();
else
R = qrJT.matrixQR().topRows(colsNum)
.triangularView<Eigen::Upper>();
}
else {
rowsNum = JG.rows();
colsNum = JG.cols();
}

#ifdef _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX
R2 = qrJT.matrixQR();
Q = qrJT.matrixQ();
#endif
}

if(debugMode==IterationLevel) {
SolverReportingManager::Manager().LogQRSystemInformation(*this, paramsNum, constrNum, rank);
if(debugMode==IterationLevel && !silent) {
SolverReportingManager::Manager().LogQRSystemInformation(*this, rowsNum, colsNum, rank);
}

#ifdef _GCS_DEBUG_SOLVER_JACOBIAN_QR_DECOMPOSITION_TRIANGULAR_MATRIX
if (J.rows() > 0) {
if (J.rows() > 0 && !silent) {
SolverReportingManager::Manager().LogMatrix("R", R);

SolverReportingManager::Manager().LogMatrix("R2", R2);
Expand Down Expand Up @@ -4151,45 +4173,71 @@ void System::makeSparseQRDecomposition( const Eigen::MatrixXd &J,
#endif // EIGEN_SPARSEQR_COMPATIBLE
}

void System::identifyDependentParametersSparseQR( const Eigen::MatrixXd &J,
void System::identifyDependentParametersDenseQR( const Eigen::MatrixXd &J,
const std::map<int,int> &jacobianconstraintmap,
const GCS::VEC_pD &pdiagnoselist,
bool silent)
{
Eigen::FullPivHouseholderQR<Eigen::MatrixXd> qrJ;
Eigen::MatrixXd Rparams;

int rank;

makeDenseQRDecomposition( J, jacobianconstraintmap, qrJ, rank, Rparams, false, true);

identifyDependentParameters(qrJ, Rparams, rank, pdiagnoselist, silent);
}

void System::identifyDependentParametersSparseQR( const Eigen::MatrixXd &J,
const std::map<int,int> &jacobianconstraintmap,
const GCS::VEC_pD &pdiagnoselist,
bool silent)
{
Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int> > SqrJ;
Eigen::MatrixXd Rparams;

int nontransprank;

makeSparseQRDecomposition( J, jacobianconstraintmap, SqrJ, nontransprank, Rparams, false, true); // do not transpose allow to diagnose parameters

identifyDependentParameters(SqrJ, Rparams, nontransprank, pdiagnoselist, silent);
}

template <typename T>
void System::identifyDependentParameters( T & qrJ,
Eigen::MatrixXd &Rparams,
int rank,
const GCS::VEC_pD &pdiagnoselist,
bool silent)
{


//int constrNum = SqrJ.rows(); // this is the other way around than for the transposed J
//int paramsNum = SqrJ.cols();

eliminateNonZerosOverPivotInUpperTriangularMatrix(Rparams, nontransprank);
eliminateNonZerosOverPivotInUpperTriangularMatrix(Rparams, rank);

#ifdef _GCS_DEBUG
if(!silent)
SolverReportingManager::Manager().LogMatrix("Rparams", Rparams);
#endif
//std::vector< std::vector<double *> > dependencyGroups(SqrJ.cols()-rank);
for (int j=nontransprank; j < SqrJ.cols(); j++) {
for (int row=0; row < nontransprank; row++) {
for (int j=rank; j < qrJ.cols(); j++) {
for (int row=0; row < rank; row++) {
if (fabs(Rparams(row,j)) > 1e-10) {
int origCol = SqrJ.colsPermutation().indices()[row];
int origCol = qrJ.colsPermutation().indices()[row];

//dependencyGroups[j-rank].push_back(pdiagnoselist[origCol]);
pdependentparameters.push_back(pdiagnoselist[origCol]);
}
}
int origCol = SqrJ.colsPermutation().indices()[j];
int origCol = qrJ.colsPermutation().indices()[j];

//dependencyGroups[j-rank].push_back(pdiagnoselist[origCol]);
pdependentparameters.push_back(pdiagnoselist[origCol]);
}
}


void System::identifyDependentGeometryParametersInTransposedJacobianDenseQRDecomposition(
const Eigen::FullPivHouseholderQR<Eigen::MatrixXd>& qrJT,
const GCS::VEC_pD &pdiagnoselist,
Expand Down
14 changes: 13 additions & 1 deletion src/Mod/Sketcher/App/planegcs/GCS.h
View file
Expand Up @@ -136,7 +136,7 @@ namespace GCS
void makeDenseQRDecomposition( const Eigen::MatrixXd &J,
const std::map<int,int> &jacobianconstraintmap,
Eigen::FullPivHouseholderQR<Eigen::MatrixXd>& qrJT,
int &rank, Eigen::MatrixXd &R);
int &rank, Eigen::MatrixXd &R, bool transposeJ = true, bool silent = false);

void makeSparseQRDecomposition( const Eigen::MatrixXd &J,
const std::map<int,int> &jacobianconstraintmap,
Expand Down Expand Up @@ -170,6 +170,18 @@ namespace GCS
const GCS::VEC_pD &pdiagnoselist,
bool silent=true);

void identifyDependentParametersDenseQR( const Eigen::MatrixXd &J,
const std::map<int,int> &jacobianconstraintmap,
const GCS::VEC_pD &pdiagnoselist,
bool silent=true);

template <typename T>
void identifyDependentParameters( T & qrJ,
Eigen::MatrixXd &Rparams,
int rank,
const GCS::VEC_pD &pdiagnoselist,
bool silent=true);

#ifdef _GCS_EXTRACT_SOLVER_SUBSYSTEM_
void extractSubsystem(SubSystem *subsys, bool isRedundantsolving);
#endif
Expand Down

0 comments on commit b2fd91c

Please sign in to comment.