@@ -33,46 +33,69 @@ namespace dg = dynamicgraph;
33
33
/* --------------------------------------------------------------------- */
34
34
namespace Eigen {
35
35
36
- void pseudoInverse ( dg::Matrix& _inputMatrix,
37
- dg::Matrix& _inverseMatrix,
38
- const double threshold = 1e-6 )
39
- {
40
- JacobiSVD<dg::Matrix> svd (_inputMatrix, ComputeThinU | ComputeThinV);
41
- JacobiSVD<dg::Matrix>::SingularValuesType m_singularValues=svd.singularValues ();
42
- JacobiSVD<dg::Matrix>::SingularValuesType singularValues_inv;
43
- singularValues_inv.resizeLike (m_singularValues);
44
- for ( long i=0 ; i<m_singularValues.size (); ++i) {
45
- if ( m_singularValues (i) > threshold )
46
- singularValues_inv (i)=1.0 /m_singularValues (i);
47
- else singularValues_inv (i)=0 ;
48
- }
49
- _inverseMatrix = (svd.matrixV ()*singularValues_inv.asDiagonal ()*svd.matrixU ().transpose ());
50
- }
36
+ void pseudoInverse ( dg::Matrix& _inputMatrix,
37
+ dg::Matrix& _inverseMatrix,
38
+ const double threshold = 1e-6 ) {
39
+ JacobiSVD<dg::Matrix> svd (_inputMatrix, ComputeThinU | ComputeThinV);
40
+ JacobiSVD<dg::Matrix>::SingularValuesType m_singularValues=svd.singularValues ();
41
+ JacobiSVD<dg::Matrix>::SingularValuesType singularValues_inv;
42
+ singularValues_inv.resizeLike (m_singularValues);
43
+ for ( long i=0 ; i<m_singularValues.size (); ++i) {
44
+ if ( m_singularValues (i) > threshold )
45
+ singularValues_inv (i)=1.0 /m_singularValues (i);
46
+ else singularValues_inv (i)=0 ;
47
+ }
48
+ _inverseMatrix = (svd.matrixV ()*singularValues_inv.asDiagonal ()*svd.matrixU ().transpose ());
49
+ }
50
+
51
+ void dampedInverse ( dg::Matrix& _inputMatrix,
52
+ dg::Matrix& _inverseMatrix,
53
+ dg::Matrix& Uref,
54
+ dg::Vector& Sref,
55
+ dg::Matrix& Vref,
56
+ const double threshold = 1e-6 ) {
57
+ sotDEBUGIN (15 );
58
+ sotDEBUG (5 ) << " Input Matrix: " <<_inputMatrix<<std::endl;
59
+ JacobiSVD<dg::Matrix> svd (_inputMatrix, ComputeThinU | ComputeThinV);
60
+ JacobiSVD<dg::Matrix>::SingularValuesType m_singularValues=svd.singularValues ();
61
+ JacobiSVD<dg::Matrix>::SingularValuesType singularValues_inv;
62
+ singularValues_inv.resizeLike (m_singularValues);
63
+ for ( long i=0 ; i<m_singularValues.size (); ++i) {
64
+ if ( m_singularValues (i) > threshold )
65
+ singularValues_inv (i)=m_singularValues (i)/(m_singularValues (i)*m_singularValues (i)+threshold*threshold);
66
+ else singularValues_inv (i)=0 ;
67
+ }
68
+ dg::Matrix matrix_U (svd.matrixU ());
69
+ dg::Matrix matrix_V (svd.matrixV ());
70
+ _inverseMatrix = (matrix_V*singularValues_inv.asDiagonal ()*matrix_U.transpose ());
71
+ Uref = matrix_U; Vref = matrix_V; Sref = m_singularValues;
72
+
73
+ sotDEBUGOUT (15 );
74
+ }
75
+
76
+ void dampedInverse ( dg::Matrix& _inputMatrix,
77
+ dg::Matrix& _inverseMatrix,
78
+ const double threshold = 1e-6 ) {
79
+ sotDEBUGIN (15 );
80
+ sotDEBUG (5 ) << " Input Matrix: " <<_inputMatrix<<std::endl;
81
+ JacobiSVD<dg::Matrix> svd (_inputMatrix, ComputeThinU | ComputeThinV);
82
+ JacobiSVD<dg::Matrix>::SingularValuesType m_singularValues=svd.singularValues ();
83
+ JacobiSVD<dg::Matrix>::SingularValuesType singularValues_inv;
84
+ singularValues_inv.resizeLike (m_singularValues);
85
+ for ( long i=0 ; i<m_singularValues.size (); ++i) {
86
+ if ( m_singularValues (i) > threshold )
87
+ singularValues_inv (i)=m_singularValues (i)/(m_singularValues (i)*m_singularValues (i)+threshold*threshold);
88
+ else singularValues_inv (i)=0 ;
89
+ }
90
+ dg::Matrix Uref (svd.matrixU ());
91
+ dg::Matrix Vref (svd.matrixV ());
92
+
93
+ _inverseMatrix = (Vref*singularValues_inv.asDiagonal ()*Uref.transpose ());
94
+
95
+ sotDEBUGOUT (15 );
96
+ }
51
97
52
- void dampedInverse ( dg::Matrix& _inputMatrix,
53
- dg::Matrix& _inverseMatrix,
54
- const double threshold = 1e-6 ,
55
- dg::Matrix* Uref = NULL ,
56
- dg::Vector* Sref = NULL ,
57
- dg::Matrix* Vref = NULL ) {
58
- JacobiSVD<dg::Matrix> svd (_inputMatrix, ComputeThinU | ComputeThinV);
59
- JacobiSVD<dg::Matrix>::SingularValuesType m_singularValues=svd.singularValues ();
60
- JacobiSVD<dg::Matrix>::SingularValuesType singularValues_inv;
61
- singularValues_inv.resizeLike (m_singularValues);
62
- for ( long i=0 ; i<m_singularValues.size (); ++i) {
63
- if ( m_singularValues (i) > threshold )
64
- singularValues_inv (i)=m_singularValues (i)/(m_singularValues (i)*m_singularValues (i)+threshold*threshold);
65
- else singularValues_inv (i)=0 ;
66
- }
67
- MatrixXd svd_matrixV = svd.matrixV ();
68
- MatrixXd svd_matrixU = svd.matrixU ();
69
-
70
- _inverseMatrix = (svd_matrixV*singularValues_inv.asDiagonal ()*svd_matrixU.transpose ());
71
98
72
- if ( Uref ) Uref = &svd_matrixU;
73
- if ( Vref ) Vref = &svd_matrixV;
74
- if ( Sref ) Sref = &singularValues_inv;
75
- }
76
99
77
100
78
101
}
0 commit comments