sofa-framework
[SofaSparseSolver] Fix the assembling of the constraint matrix
It seems that there was an error in the code of the SparseLDLSolver in the assembling of the constraint matrix. I tried to see the execution of this function in the scene InextensiblePendulum.scn and it produced a wrong result at the first step (with the current master branch). After this fix, it starts with a coherent behaviour but then it still fails.
By submitting this pull request, I acknowledge that I have read, understand, and agree SOFA Developer Certificate of Origin (DCO).
Reviewers will merge this pull-request only if
it builds with SUCCESS for all platforms on the CI.
it does not generate new warnings.
it does not generate new unit test failures.
it does not generate new scene test failures.
it does not break API compatibility.
it is more than 1 week old (or has fast-merge label).
-
Indeed I changed the linear solver used in the scene. In the master branch, the code actually compute another matrix which is not the one we need. Hence the simulation fail at the first step. For this branch the smilation failed when the message "D(k,k) is zero" appear. I don't know if there is an issue with this scene or with the solver.
-
In the code, we exploit the LDL factorization in order to solve only one triangular system (and not two). We are allowed to do that by rewriting the matricial product we need to compute. It will then write as A *D^-1 * A^t. We need to apply D^-1 and it won't commute, that is why there is two matrices in memory (we could store only one with the Cholesky factorization).
@dessteph On my setup with the master branch, I manage to run the simulation several tens of iterations (see screenshot). So I would like to know what you mean by "the simulation fails".
After a certain number of iterations (on my computer), an error message appears in the terminal and the simulation becomes unstable (it explodes) and that doesn't append with LU or Cholesky.
@dessteph I suggest you start again your investigation but taking into account https://github.com/sofa-framework/sofa/pull/3072. In other words, in your CMake, enable the variable METIS-GKLIB_GKRAND. This way we should have comparable results across different architectures.
Rebase on #3050 and then create a unit test on SparseLDLSolver: doAddJMInvJtLocal
@dessteph I prepared the test in Sofa/Component/LinearSolver/Direct/tests/SparseLDLSolver_test.cpp. Could you add the check of the result at the end? You will need to use another software (matlab?) to compute the ground truth and compare it to the variable JMInvJt containing our result. Thanks
@dessteph What is the status of this pull request? Can you provide details why we should continue working on it or not?
I made an unitary test to compare the results of the LDL solver to those of the LU solver. The difference goes up to 9e-13. I would have like to compare LU and Cholesky to see if that difference is realy significative.
And what about the result of the LDL solver with and without your changes?
It seems it didn't need any fixing, there is a math error but compensated by a roundabout use of the CSC format.