@@ -56,10 +56,9 @@ namespace roboptim
5656
5757 // compute number of non zeros elements in jacobian constraint.
5858 nnz_jac_g = 0 ;
59- typedef std::vector<differentiableFunctionPtr_t>::const_iterator
60- citer_t ;
61- for (citer_t it = differentiableConstraintFunctions_.begin ();
62- it != differentiableConstraintFunctions_.end (); ++it)
59+ typedef differentiableConstraints_t::const_iterator citer_t ;
60+ for (citer_t it = differentiableConstraints_.begin ();
61+ it != differentiableConstraints_.end (); ++it)
6362 {
6463 nnz_jac_g += (*it)->jacobian (x).nonZeros ();
6564 }
@@ -84,12 +83,12 @@ namespace roboptim
8483 assert (costFunction_->inputSize () == n_);
8584 assert (constraintsOutputSize () == m);
8685
87- if (!jacobian_ )
86+ if (!jacobianBuf_ )
8887 {
89- jacobian_ = differentiableFunction_t::jacobian_t
88+ jacobianBuf_ = differentiableFunction_t::jacobian_t
9089 (static_cast <differentiableFunction_t::matrix_t ::Index> (constraintsOutputSize ()),
9190 costFunction_->inputSize ());
92- jacobian_ ->reserve (nele_jac);
91+ jacobianBuf_ ->reserve (nele_jac);
9392 }
9493
9594 if (!values)
@@ -98,22 +97,24 @@ namespace roboptim
9897 (logger, " Looking for non-zeros elements." );
9998 LOG4CXX_TRACE (logger, " nele_jac = " << nele_jac);
10099
100+ // Clear just in case
101+ constraintJacobians_.clear ();
102+
101103 // Emptying iRow/jCol arrays.
102- memset (iRow, 0 , static_cast <std::size_t > (nele_jac) * sizeof (Index));
103- memset (jCol, 0 , static_cast <std::size_t > (nele_jac) * sizeof (Index));
104+ std:: memset (iRow, 0 , static_cast <std::size_t > (nele_jac) * sizeof (Index));
105+ std:: memset (jCol, 0 , static_cast <std::size_t > (nele_jac) * sizeof (Index));
104106
105107 // First evaluate the constraints in zero to build the
106108 // constraints jacobian.
107109 int idx = 0 ;
108- typedef std::vector<differentiableFunctionPtr_t>::const_iterator
109- citer_t ;
110+ typedef differentiableConstraints_t::const_iterator citer_t ;
110111 unsigned constraintId = 0 ;
111112
112113 typedef Eigen::Triplet<double > triplet_t ;
113114 std::vector<triplet_t > coefficients;
114115
115- for (citer_t it = differentiableConstraintFunctions_ .begin ();
116- it != differentiableConstraintFunctions_ .end ();
116+ for (citer_t it = differentiableConstraints_ .begin ();
117+ it != differentiableConstraints_ .end ();
117118 ++it, ++constraintId)
118119 {
119120 LOG4CXX_TRACE
@@ -154,10 +155,13 @@ namespace roboptim
154155 else // other use initial guess.
155156 x = *(solver_.problem ().startingPoint ());
156157
157- differentiableFunction_t::jacobian_t jacobian = (*it)->jacobian (x);
158- for (int k = 0 ; k < jacobian.outerSize (); ++k)
158+ constraintJacobians_.push_back ((*it)->jacobian (x));
159+ differentiableFunction_t::jacobian_t & tmp_jac = constraintJacobians_.back ();
160+ tmp_jac.makeCompressed ();
161+
162+ for (int k = 0 ; k < tmp_jac.outerSize (); ++k)
159163 for (differentiableFunction_t::jacobian_t ::InnerIterator
160- it (jacobian , k); it; ++it)
164+ it (tmp_jac , k); it; ++it)
161165 {
162166 const int row = static_cast <int > (idx + it.row ());
163167 const int col = static_cast <int > (it.col ());
@@ -167,18 +171,18 @@ namespace roboptim
167171 idx += (*it)->outputSize ();
168172 }
169173
170- jacobian_ ->setFromTriplets
174+ jacobianBuf_ ->setFromTriplets
171175 (coefficients.begin (), coefficients.end ());
172176
173177 LOG4CXX_TRACE
174- (logger, " full problem jacobian...\n " << *jacobian_ );
178+ (logger, " full problem jacobian...\n " << *jacobianBuf_ );
175179
176180 // Then look for non-zero values.
177181 LOG4CXX_TRACE (logger, " filling iRow and jCol..." );
178182 idx = 0 ;
179183
180- for (int k = 0 ; k < jacobian_ ->outerSize (); ++k)
181- for (differentiableFunction_t::jacobian_t ::InnerIterator it (*jacobian_ , k);
184+ for (int k = 0 ; k < jacobianBuf_ ->outerSize (); ++k)
185+ for (differentiableFunction_t::jacobian_t ::InnerIterator it (*jacobianBuf_ , k);
182186 it; ++it)
183187 {
184188 iRow[idx] = it.row (), jCol[idx] = it.col ();
@@ -196,34 +200,53 @@ namespace roboptim
196200
197201 Eigen::Map<const function_t ::vector_t > x_ (x, n);
198202
199- typedef std::vector<differentiableFunctionPtr_t>::const_iterator
200- citer_t ;
203+ typedef differentiableConstraints_t::const_iterator citer_t ;
201204
202- int idx = 0 ;
203205 int constraintId = 0 ;
204- for (citer_t it = differentiableConstraintFunctions_ .begin ();
205- it != differentiableConstraintFunctions_ .end (); ++it)
206+ for (citer_t it = differentiableConstraints_ .begin ();
207+ it != differentiableConstraints_ .end (); ++it, constraintId++ )
206208 {
207- // TODO: use middleRows once Eigen is fixed
208- // TODO: avoid allocation here (may be solved with
209- // http://eigen.tuxfamily.org/bz/show_bug.cgi?id=910)
210- copySparseBlock (*jacobian_, (*it)->jacobian (x_), idx, 0 );
211- idx += (*it)->outputSize ();
209+ typename differentiableFunction_t::matrix_t & jac = constraintJacobians_[constraintId];
210+ // Set the Jacobian to 0 while keeping its structure
211+ jac *= 0 .;
212+ (*it)->jacobian (jac, x_);
212213
213214 IpoptCheckGradient
214215 (*(*it), 0 , x_,
215- constraintId++ , solver_);
216+ static_cast < int > ( constraintId) , solver_);
216217 }
217218
218- // Copy jacobian values from internal sparse matrix.
219- idx = 0 ;
220- for (int k = 0 ; k < jacobian_->outerSize (); ++k)
221- for (differentiableFunction_t::jacobian_t ::InnerIterator it (*jacobian_, k);
222- it; ++it)
223- {
224- assert (idx < nele_jac);
225- values[idx++] = it.value ();
226- }
219+ // Copy jacobian values from internal sparse matrices.
220+ int idx = 0 ;
221+
222+ if (StorageOrder == Eigen::ColMajor)
223+ {
224+ for (int k = 0 ; k < solver_.problem ().function ().inputSize (); ++k)
225+ for (constraintJacobians_t::const_iterator
226+ g = constraintJacobians_.begin ();
227+ g != constraintJacobians_.end (); ++g)
228+ {
229+ for (differentiableFunction_t::jacobian_t ::InnerIterator it (*g, k);
230+ it; ++it)
231+ {
232+ assert (idx < nele_jac);
233+ values[idx++] = it.value ();
234+ }
235+ }
236+ }
237+ else
238+ {
239+ for (constraintJacobians_t::const_iterator
240+ g = constraintJacobians_.begin ();
241+ g != constraintJacobians_.end (); ++g)
242+ for (int k = 0 ; k < g->outerSize (); ++k)
243+ for (differentiableFunction_t::jacobian_t ::InnerIterator it (*g, k);
244+ it; ++it)
245+ {
246+ assert (idx < nele_jac);
247+ values[idx++] = it.value ();
248+ }
249+ }
227250
228251 return true ;
229252 }
0 commit comments