svn merge ^/trunk/blender -r41226:41227 .
[blender.git] / intern / itasc / WSDLSSolver.cpp
1 /** \file itasc/WSDLSSolver.cpp
2  *  \ingroup itasc
3  */
4  * WDLSSolver.hpp.cpp
5  *
6  *  Created on: Jan 8, 2009
7  *      Author: rubensmits
8  */
9
10 #include "WSDLSSolver.hpp"
11 #include "kdl/utilities/svd_eigen_HH.hpp"
12 #include <cstdio>
13
14 namespace iTaSC {
15
16 WSDLSSolver::WSDLSSolver() :
17         m_ns(0), m_nc(0), m_nq(0)
18
19 {
20         // default maximum speed: 50 rad/s
21         m_qmax = 50.0;
22 }
23
24 WSDLSSolver::~WSDLSSolver() {
25 }
26
27 bool WSDLSSolver::init(unsigned int _nq, unsigned int _nc, const std::vector<bool>& gc)
28 {
29         if (_nc == 0 || _nq == 0 || gc.size() != _nc)
30                 return false;
31         m_nc = _nc;
32         m_nq = _nq;
33         m_ns = std::min(m_nc,m_nq);
34     m_AWq = e_zero_matrix(m_nc,m_nq);
35     m_WyAWq = e_zero_matrix(m_nc,m_nq);
36     m_WyAWqt = e_zero_matrix(m_nq,m_nc);
37         m_S = e_zero_vector(std::max(m_nc,m_nq));
38         m_Wy_ydot = e_zero_vector(m_nc);
39         m_ytask = gc;
40         if (m_nq > m_nc) {
41                 m_transpose = true;
42             m_temp = e_zero_vector(m_nc);
43             m_U = e_zero_matrix(m_nc,m_nc);
44                 m_V = e_zero_matrix(m_nq,m_nc);
45             m_WqV = e_zero_matrix(m_nq,m_nc);
46         } else {
47                 m_transpose = false;
48             m_temp = e_zero_vector(m_nq);
49             m_U = e_zero_matrix(m_nc,m_nq);
50                 m_V = e_zero_matrix(m_nq,m_nq);
51             m_WqV = e_zero_matrix(m_nq,m_nq);
52         }
53     return true;
54 }
55
56 bool WSDLSSolver::solve(const e_matrix& A, const e_vector& Wy, const e_vector& ydot, const e_matrix& Wq, e_vector& qdot, e_scalar& nlcoef)
57 {
58         unsigned int i, j, l;
59         e_scalar N, M;
60
61         // Create the Weighted jacobian
62     m_AWq = (A*Wq).lazy();
63         for (i=0; i<m_nc; i++)
64                 m_WyAWq.row(i) = Wy(i)*m_AWq.row(i);
65
66     // Compute the SVD of the weighted jacobian
67         int ret;
68         if (m_transpose) {
69                 m_WyAWqt = m_WyAWq.transpose();
70                 ret = KDL::svd_eigen_HH(m_WyAWqt,m_V,m_S,m_U,m_temp);
71         } else {
72                 ret = KDL::svd_eigen_HH(m_WyAWq,m_U,m_S,m_V,m_temp);
73         }
74     if(ret<0)
75         return false;
76
77         m_Wy_ydot = Wy.cwise() * ydot;
78     m_WqV = (Wq*m_V).lazy();
79         qdot.setZero();
80         e_scalar maxDeltaS = e_scalar(0.0);
81         e_scalar prevS = e_scalar(0.0);
82         e_scalar maxS = e_scalar(1.0);
83         for(i=0;i<m_ns;++i) {
84                 e_scalar norm, mag, alpha, _qmax, Sinv, vmax, damp;
85                 e_scalar S = m_S(i);
86                 bool prev;
87                 if (S < KDL::epsilon)
88                         break;
89                 Sinv = e_scalar(1.)/S;
90                 if (i > 0) {
91                         if ((prevS-S) > maxDeltaS) {
92                                 maxDeltaS = (prevS-S);
93                                 maxS = prevS;
94                         }
95                 }
96                 N = M = e_scalar(0.);
97                 for (l=0, prev=m_ytask[0], norm=e_scalar(0.); l<m_nc; l++) {
98                         if (prev == m_ytask[l]) {
99                                 norm += m_U(l,i)*m_U(l,i);
100                         } else {
101                                 N += std::sqrt(norm);
102                                 norm = m_U(l,i)*m_U(l,i);
103                         }
104                         prev = m_ytask[l];
105                 }
106                 N += std::sqrt(norm);
107                 for (j=0; j<m_nq; j++) {
108                         for (l=0, prev=m_ytask[0], norm=e_scalar(0.), mag=e_scalar(0.); l<m_nc; l++) {
109                                 if (prev == m_ytask[l]) {
110                                         norm += m_WyAWq(l,j)*m_WyAWq(l,j);
111                                 } else {
112                                         mag += std::sqrt(norm);
113                                         norm = m_WyAWq(l,j)*m_WyAWq(l,j);
114                                 }
115                                 prev = m_ytask[l];
116                         }
117                         mag += std::sqrt(norm);
118                         M += fabs(m_V(j,i))*mag;
119                 }
120                 M *= Sinv;
121                 alpha = m_U.col(i).dot(m_Wy_ydot);
122                 _qmax = (N < M) ? m_qmax*N/M : m_qmax;
123                 vmax = m_WqV.col(i).cwise().abs().maxCoeff();
124                 norm = fabs(Sinv*alpha*vmax);
125                 if (norm > _qmax) {
126                         damp = Sinv*alpha*_qmax/norm;
127                 } else {
128                         damp = Sinv*alpha;
129                 }
130                 qdot += m_WqV.col(i)*damp;
131                 prevS = S;
132         }
133         if (maxDeltaS == e_scalar(0.0))
134                 nlcoef = e_scalar(KDL::epsilon);
135         else
136                 nlcoef = (maxS-maxDeltaS)/maxS;
137     return true;
138 }
139
140 }