Ceres: Update to the latest actual version
[blender.git] / extern / ceres / internal / ceres / reorder_program.cc
1 // Ceres Solver - A fast non-linear least squares minimizer
2 // Copyright 2015 Google Inc. All rights reserved.
3 // http://ceres-solver.org/
4 //
5 // Redistribution and use in source and binary forms, with or without
6 // modification, are permitted provided that the following conditions are met:
7 //
8 // * Redistributions of source code must retain the above copyright notice,
9 //   this list of conditions and the following disclaimer.
10 // * Redistributions in binary form must reproduce the above copyright notice,
11 //   this list of conditions and the following disclaimer in the documentation
12 //   and/or other materials provided with the distribution.
13 // * Neither the name of Google Inc. nor the names of its contributors may be
14 //   used to endorse or promote products derived from this software without
15 //   specific prior written permission.
16 //
17 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
21 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 // POSSIBILITY OF SUCH DAMAGE.
28 //
29 // Author: sameeragarwal@google.com (Sameer Agarwal)
30
31 #include "ceres/reorder_program.h"
32
33 #include <algorithm>
34 #include <numeric>
35 #include <vector>
36
37 #include "ceres/cxsparse.h"
38 #include "ceres/internal/port.h"
39 #include "ceres/ordered_groups.h"
40 #include "ceres/parameter_block.h"
41 #include "ceres/parameter_block_ordering.h"
42 #include "ceres/problem_impl.h"
43 #include "ceres/program.h"
44 #include "ceres/residual_block.h"
45 #include "ceres/solver.h"
46 #include "ceres/suitesparse.h"
47 #include "ceres/triplet_sparse_matrix.h"
48 #include "ceres/types.h"
49 #include "Eigen/SparseCore"
50
51 #ifdef CERES_USE_EIGEN_SPARSE
52 #include "Eigen/OrderingMethods"
53 #endif
54
55 #include "glog/logging.h"
56
57 namespace ceres {
58 namespace internal {
59
60 using std::map;
61 using std::set;
62 using std::string;
63 using std::vector;
64
65 namespace {
66
67 // Find the minimum index of any parameter block to the given
68 // residual.  Parameter blocks that have indices greater than
69 // size_of_first_elimination_group are considered to have an index
70 // equal to size_of_first_elimination_group.
71 static int MinParameterBlock(const ResidualBlock* residual_block,
72                              int size_of_first_elimination_group) {
73   int min_parameter_block_position = size_of_first_elimination_group;
74   for (int i = 0; i < residual_block->NumParameterBlocks(); ++i) {
75     ParameterBlock* parameter_block = residual_block->parameter_blocks()[i];
76     if (!parameter_block->IsConstant()) {
77       CHECK_NE(parameter_block->index(), -1)
78           << "Did you forget to call Program::SetParameterOffsetsAndIndex()? "
79           << "This is a Ceres bug; please contact the developers!";
80       min_parameter_block_position = std::min(parameter_block->index(),
81                                               min_parameter_block_position);
82     }
83   }
84   return min_parameter_block_position;
85 }
86
87 #if EIGEN_VERSION_AT_LEAST(3, 2, 2) && defined(CERES_USE_EIGEN_SPARSE)
88 Eigen::SparseMatrix<int> CreateBlockJacobian(
89     const TripletSparseMatrix& block_jacobian_transpose) {
90   typedef Eigen::SparseMatrix<int> SparseMatrix;
91   typedef Eigen::Triplet<int> Triplet;
92
93   const int* rows = block_jacobian_transpose.rows();
94   const int* cols = block_jacobian_transpose.cols();
95   int num_nonzeros = block_jacobian_transpose.num_nonzeros();
96   vector<Triplet> triplets;
97   triplets.reserve(num_nonzeros);
98   for (int i = 0; i < num_nonzeros; ++i) {
99     triplets.push_back(Triplet(cols[i], rows[i], 1));
100   }
101
102   SparseMatrix block_jacobian(block_jacobian_transpose.num_cols(),
103                               block_jacobian_transpose.num_rows());
104   block_jacobian.setFromTriplets(triplets.begin(), triplets.end());
105   return block_jacobian;
106 }
107 #endif
108
109 void OrderingForSparseNormalCholeskyUsingSuiteSparse(
110     const TripletSparseMatrix& tsm_block_jacobian_transpose,
111     const vector<ParameterBlock*>& parameter_blocks,
112     const ParameterBlockOrdering& parameter_block_ordering,
113     int* ordering) {
114 #ifdef CERES_NO_SUITESPARSE
115   LOG(FATAL) << "Congratulations, you found a Ceres bug! "
116              << "Please report this error to the developers.";
117 #else
118   SuiteSparse ss;
119   cholmod_sparse* block_jacobian_transpose =
120       ss.CreateSparseMatrix(
121           const_cast<TripletSparseMatrix*>(&tsm_block_jacobian_transpose));
122
123   // No CAMD or the user did not supply a useful ordering, then just
124   // use regular AMD.
125   if (parameter_block_ordering.NumGroups() <= 1 ||
126       !SuiteSparse::IsConstrainedApproximateMinimumDegreeOrderingAvailable()) {
127     ss.ApproximateMinimumDegreeOrdering(block_jacobian_transpose, &ordering[0]);
128   } else {
129     vector<int> constraints;
130     for (int i = 0; i < parameter_blocks.size(); ++i) {
131       constraints.push_back(
132           parameter_block_ordering.GroupId(
133               parameter_blocks[i]->mutable_user_state()));
134     }
135
136     // Renumber the entries of constraints to be contiguous integers
137     // as CAMD requires that the group ids be in the range [0,
138     // parameter_blocks.size() - 1].
139     MapValuesToContiguousRange(constraints.size(), &constraints[0]);
140     ss.ConstrainedApproximateMinimumDegreeOrdering(block_jacobian_transpose,
141                                                    &constraints[0],
142                                                    ordering);
143   }
144
145   VLOG(2) << "Block ordering stats: "
146           << " flops: " << ss.mutable_cc()->fl
147           << " lnz  : " << ss.mutable_cc()->lnz
148           << " anz  : " << ss.mutable_cc()->anz;
149
150   ss.Free(block_jacobian_transpose);
151 #endif  // CERES_NO_SUITESPARSE
152 }
153
154 void OrderingForSparseNormalCholeskyUsingCXSparse(
155     const TripletSparseMatrix& tsm_block_jacobian_transpose,
156     int* ordering) {
157 #ifdef CERES_NO_CXSPARSE
158   LOG(FATAL) << "Congratulations, you found a Ceres bug! "
159              << "Please report this error to the developers.";
160 #else  // CERES_NO_CXSPARSE
161   // CXSparse works with J'J instead of J'. So compute the block
162   // sparsity for J'J and compute an approximate minimum degree
163   // ordering.
164   CXSparse cxsparse;
165   cs_di* block_jacobian_transpose;
166   block_jacobian_transpose =
167       cxsparse.CreateSparseMatrix(
168             const_cast<TripletSparseMatrix*>(&tsm_block_jacobian_transpose));
169   cs_di* block_jacobian = cxsparse.TransposeMatrix(block_jacobian_transpose);
170   cs_di* block_hessian =
171       cxsparse.MatrixMatrixMultiply(block_jacobian_transpose, block_jacobian);
172   cxsparse.Free(block_jacobian);
173   cxsparse.Free(block_jacobian_transpose);
174
175   cxsparse.ApproximateMinimumDegreeOrdering(block_hessian, ordering);
176   cxsparse.Free(block_hessian);
177 #endif  // CERES_NO_CXSPARSE
178 }
179
180
181 #if EIGEN_VERSION_AT_LEAST(3, 2, 2)
182 void OrderingForSparseNormalCholeskyUsingEigenSparse(
183     const TripletSparseMatrix& tsm_block_jacobian_transpose,
184     int* ordering) {
185 #ifndef CERES_USE_EIGEN_SPARSE
186   LOG(FATAL) <<
187       "SPARSE_NORMAL_CHOLESKY cannot be used with EIGEN_SPARSE "
188       "because Ceres was not built with support for "
189       "Eigen's SimplicialLDLT decomposition. "
190       "This requires enabling building with -DEIGENSPARSE=ON.";
191 #else
192
193   // This conversion from a TripletSparseMatrix to a Eigen::Triplet
194   // matrix is unfortunate, but unavoidable for now. It is not a
195   // significant performance penalty in the grand scheme of
196   // things. The right thing to do here would be to get a compressed
197   // row sparse matrix representation of the jacobian and go from
198   // there. But that is a project for another day.
199   typedef Eigen::SparseMatrix<int> SparseMatrix;
200
201   const SparseMatrix block_jacobian =
202       CreateBlockJacobian(tsm_block_jacobian_transpose);
203   const SparseMatrix block_hessian =
204       block_jacobian.transpose() * block_jacobian;
205
206   Eigen::AMDOrdering<int> amd_ordering;
207   Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic, int> perm;
208   amd_ordering(block_hessian, perm);
209   for (int i = 0; i < block_hessian.rows(); ++i) {
210     ordering[i] = perm.indices()[i];
211   }
212 #endif  // CERES_USE_EIGEN_SPARSE
213 }
214 #endif
215
216 }  // namespace
217
218 bool ApplyOrdering(const ProblemImpl::ParameterMap& parameter_map,
219                    const ParameterBlockOrdering& ordering,
220                    Program* program,
221                    string* error) {
222   const int num_parameter_blocks =  program->NumParameterBlocks();
223   if (ordering.NumElements() != num_parameter_blocks) {
224     *error = StringPrintf("User specified ordering does not have the same "
225                           "number of parameters as the problem. The problem"
226                           "has %d blocks while the ordering has %d blocks.",
227                           num_parameter_blocks,
228                           ordering.NumElements());
229     return false;
230   }
231
232   vector<ParameterBlock*>* parameter_blocks =
233       program->mutable_parameter_blocks();
234   parameter_blocks->clear();
235
236   const map<int, set<double*> >& groups = ordering.group_to_elements();
237   for (map<int, set<double*> >::const_iterator group_it = groups.begin();
238        group_it != groups.end();
239        ++group_it) {
240     const set<double*>& group = group_it->second;
241     for (set<double*>::const_iterator parameter_block_ptr_it = group.begin();
242          parameter_block_ptr_it != group.end();
243          ++parameter_block_ptr_it) {
244       ProblemImpl::ParameterMap::const_iterator parameter_block_it =
245           parameter_map.find(*parameter_block_ptr_it);
246       if (parameter_block_it == parameter_map.end()) {
247         *error = StringPrintf("User specified ordering contains a pointer "
248                               "to a double that is not a parameter block in "
249                               "the problem. The invalid double is in group: %d",
250                               group_it->first);
251         return false;
252       }
253       parameter_blocks->push_back(parameter_block_it->second);
254     }
255   }
256   return true;
257 }
258
259 bool LexicographicallyOrderResidualBlocks(
260     const int size_of_first_elimination_group,
261     Program* program,
262     string* error) {
263   CHECK_GE(size_of_first_elimination_group, 1)
264       << "Congratulations, you found a Ceres bug! Please report this error "
265       << "to the developers.";
266
267   // Create a histogram of the number of residuals for each E block. There is an
268   // extra bucket at the end to catch all non-eliminated F blocks.
269   vector<int> residual_blocks_per_e_block(size_of_first_elimination_group + 1);
270   vector<ResidualBlock*>* residual_blocks = program->mutable_residual_blocks();
271   vector<int> min_position_per_residual(residual_blocks->size());
272   for (int i = 0; i < residual_blocks->size(); ++i) {
273     ResidualBlock* residual_block = (*residual_blocks)[i];
274     int position = MinParameterBlock(residual_block,
275                                      size_of_first_elimination_group);
276     min_position_per_residual[i] = position;
277     DCHECK_LE(position, size_of_first_elimination_group);
278     residual_blocks_per_e_block[position]++;
279   }
280
281   // Run a cumulative sum on the histogram, to obtain offsets to the start of
282   // each histogram bucket (where each bucket is for the residuals for that
283   // E-block).
284   vector<int> offsets(size_of_first_elimination_group + 1);
285   std::partial_sum(residual_blocks_per_e_block.begin(),
286                    residual_blocks_per_e_block.end(),
287                    offsets.begin());
288   CHECK_EQ(offsets.back(), residual_blocks->size())
289       << "Congratulations, you found a Ceres bug! Please report this error "
290       << "to the developers.";
291
292   CHECK(find(residual_blocks_per_e_block.begin(),
293              residual_blocks_per_e_block.end() - 1, 0) !=
294         residual_blocks_per_e_block.end())
295       << "Congratulations, you found a Ceres bug! Please report this error "
296       << "to the developers.";
297
298   // Fill in each bucket with the residual blocks for its corresponding E block.
299   // Each bucket is individually filled from the back of the bucket to the front
300   // of the bucket. The filling order among the buckets is dictated by the
301   // residual blocks. This loop uses the offsets as counters; subtracting one
302   // from each offset as a residual block is placed in the bucket. When the
303   // filling is finished, the offset pointerts should have shifted down one
304   // entry (this is verified below).
305   vector<ResidualBlock*> reordered_residual_blocks(
306       (*residual_blocks).size(), static_cast<ResidualBlock*>(NULL));
307   for (int i = 0; i < residual_blocks->size(); ++i) {
308     int bucket = min_position_per_residual[i];
309
310     // Decrement the cursor, which should now point at the next empty position.
311     offsets[bucket]--;
312
313     // Sanity.
314     CHECK(reordered_residual_blocks[offsets[bucket]] == NULL)
315         << "Congratulations, you found a Ceres bug! Please report this error "
316         << "to the developers.";
317
318     reordered_residual_blocks[offsets[bucket]] = (*residual_blocks)[i];
319   }
320
321   // Sanity check #1: The difference in bucket offsets should match the
322   // histogram sizes.
323   for (int i = 0; i < size_of_first_elimination_group; ++i) {
324     CHECK_EQ(residual_blocks_per_e_block[i], offsets[i + 1] - offsets[i])
325         << "Congratulations, you found a Ceres bug! Please report this error "
326         << "to the developers.";
327   }
328   // Sanity check #2: No NULL's left behind.
329   for (int i = 0; i < reordered_residual_blocks.size(); ++i) {
330     CHECK(reordered_residual_blocks[i] != NULL)
331         << "Congratulations, you found a Ceres bug! Please report this error "
332         << "to the developers.";
333   }
334
335   // Now that the residuals are collected by E block, swap them in place.
336   swap(*program->mutable_residual_blocks(), reordered_residual_blocks);
337   return true;
338 }
339
340 // Pre-order the columns corresponding to the schur complement if
341 // possible.
342 void MaybeReorderSchurComplementColumnsUsingSuiteSparse(
343     const ParameterBlockOrdering& parameter_block_ordering,
344     Program* program) {
345 #ifndef CERES_NO_SUITESPARSE
346   SuiteSparse ss;
347   if (!SuiteSparse::IsConstrainedApproximateMinimumDegreeOrderingAvailable()) {
348     return;
349   }
350
351   vector<int> constraints;
352   vector<ParameterBlock*>& parameter_blocks =
353       *(program->mutable_parameter_blocks());
354
355   for (int i = 0; i < parameter_blocks.size(); ++i) {
356     constraints.push_back(
357         parameter_block_ordering.GroupId(
358             parameter_blocks[i]->mutable_user_state()));
359   }
360
361   // Renumber the entries of constraints to be contiguous integers as
362   // CAMD requires that the group ids be in the range [0,
363   // parameter_blocks.size() - 1].
364   MapValuesToContiguousRange(constraints.size(), &constraints[0]);
365
366   // Compute a block sparse presentation of J'.
367   scoped_ptr<TripletSparseMatrix> tsm_block_jacobian_transpose(
368       program->CreateJacobianBlockSparsityTranspose());
369
370   cholmod_sparse* block_jacobian_transpose =
371       ss.CreateSparseMatrix(tsm_block_jacobian_transpose.get());
372
373   vector<int> ordering(parameter_blocks.size(), 0);
374   ss.ConstrainedApproximateMinimumDegreeOrdering(block_jacobian_transpose,
375                                                  &constraints[0],
376                                                  &ordering[0]);
377   ss.Free(block_jacobian_transpose);
378
379   const vector<ParameterBlock*> parameter_blocks_copy(parameter_blocks);
380   for (int i = 0; i < program->NumParameterBlocks(); ++i) {
381     parameter_blocks[i] = parameter_blocks_copy[ordering[i]];
382   }
383
384   program->SetParameterOffsetsAndIndex();
385 #endif
386 }
387
388 void MaybeReorderSchurComplementColumnsUsingEigen(
389     const int size_of_first_elimination_group,
390     const ProblemImpl::ParameterMap& parameter_map,
391     Program* program) {
392 #if !EIGEN_VERSION_AT_LEAST(3, 2, 2) || !defined(CERES_USE_EIGEN_SPARSE)
393   return;
394 #else
395
396   scoped_ptr<TripletSparseMatrix> tsm_block_jacobian_transpose(
397       program->CreateJacobianBlockSparsityTranspose());
398
399   typedef Eigen::SparseMatrix<int> SparseMatrix;
400   const SparseMatrix block_jacobian =
401       CreateBlockJacobian(*tsm_block_jacobian_transpose);
402   const int num_rows = block_jacobian.rows();
403   const int num_cols = block_jacobian.cols();
404
405   // Vertically partition the jacobian in parameter blocks of type E
406   // and F.
407   const SparseMatrix E =
408       block_jacobian.block(0,
409                            0,
410                            num_rows,
411                            size_of_first_elimination_group);
412   const SparseMatrix F =
413       block_jacobian.block(0,
414                            size_of_first_elimination_group,
415                            num_rows,
416                            num_cols - size_of_first_elimination_group);
417
418   // Block sparsity pattern of the schur complement.
419   const SparseMatrix block_schur_complement =
420       F.transpose() * F - F.transpose() * E * E.transpose() * F;
421
422   Eigen::AMDOrdering<int> amd_ordering;
423   Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic, int> perm;
424   amd_ordering(block_schur_complement, perm);
425
426   const vector<ParameterBlock*>& parameter_blocks = program->parameter_blocks();
427   vector<ParameterBlock*> ordering(num_cols);
428
429   // The ordering of the first size_of_first_elimination_group does
430   // not matter, so we preserve the existing ordering.
431   for (int i = 0; i < size_of_first_elimination_group; ++i) {
432     ordering[i] = parameter_blocks[i];
433   }
434
435   // For the rest of the blocks, use the ordering computed using AMD.
436   for (int i = 0; i < block_schur_complement.cols(); ++i) {
437     ordering[size_of_first_elimination_group + i] =
438         parameter_blocks[size_of_first_elimination_group + perm.indices()[i]];
439   }
440
441   swap(*program->mutable_parameter_blocks(), ordering);
442   program->SetParameterOffsetsAndIndex();
443 #endif
444 }
445
446 bool ReorderProgramForSchurTypeLinearSolver(
447     const LinearSolverType linear_solver_type,
448     const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
449     const ProblemImpl::ParameterMap& parameter_map,
450     ParameterBlockOrdering* parameter_block_ordering,
451     Program* program,
452     string* error) {
453   if (parameter_block_ordering->NumElements() !=
454       program->NumParameterBlocks()) {
455     *error = StringPrintf(
456         "The program has %d parameter blocks, but the parameter block "
457         "ordering has %d parameter blocks.",
458         program->NumParameterBlocks(),
459         parameter_block_ordering->NumElements());
460     return false;
461   }
462
463   if (parameter_block_ordering->NumGroups() == 1) {
464     // If the user supplied an parameter_block_ordering with just one
465     // group, it is equivalent to the user supplying NULL as an
466     // parameter_block_ordering. Ceres is completely free to choose the
467     // parameter block ordering as it sees fit. For Schur type solvers,
468     // this means that the user wishes for Ceres to identify the
469     // e_blocks, which we do by computing a maximal independent set.
470     vector<ParameterBlock*> schur_ordering;
471     const int size_of_first_elimination_group =
472         ComputeStableSchurOrdering(*program, &schur_ordering);
473
474     CHECK_EQ(schur_ordering.size(), program->NumParameterBlocks())
475         << "Congratulations, you found a Ceres bug! Please report this error "
476         << "to the developers.";
477
478     // Update the parameter_block_ordering object.
479     for (int i = 0; i < schur_ordering.size(); ++i) {
480       double* parameter_block = schur_ordering[i]->mutable_user_state();
481       const int group_id = (i < size_of_first_elimination_group) ? 0 : 1;
482       parameter_block_ordering->AddElementToGroup(parameter_block, group_id);
483     }
484
485     // We could call ApplyOrdering but this is cheaper and
486     // simpler.
487     swap(*program->mutable_parameter_blocks(), schur_ordering);
488   } else {
489     // The user provided an ordering with more than one elimination
490     // group.
491
492     // Verify that the first elimination group is an independent set.
493     const set<double*>& first_elimination_group =
494         parameter_block_ordering
495         ->group_to_elements()
496         .begin()
497         ->second;
498     if (!program->IsParameterBlockSetIndependent(first_elimination_group)) {
499       *error =
500           StringPrintf("The first elimination group in the parameter block "
501                        "ordering of size %zd is not an independent set",
502                        first_elimination_group.size());
503       return false;
504     }
505
506     if (!ApplyOrdering(parameter_map,
507                        *parameter_block_ordering,
508                        program,
509                        error)) {
510       return false;
511     }
512   }
513
514   program->SetParameterOffsetsAndIndex();
515
516   const int size_of_first_elimination_group =
517       parameter_block_ordering->group_to_elements().begin()->second.size();
518
519   if (linear_solver_type == SPARSE_SCHUR) {
520     if (sparse_linear_algebra_library_type == SUITE_SPARSE) {
521       MaybeReorderSchurComplementColumnsUsingSuiteSparse(
522           *parameter_block_ordering,
523           program);
524     } else if (sparse_linear_algebra_library_type == EIGEN_SPARSE) {
525       MaybeReorderSchurComplementColumnsUsingEigen(
526           size_of_first_elimination_group,
527           parameter_map,
528           program);
529     }
530   }
531
532   // Schur type solvers also require that their residual blocks be
533   // lexicographically ordered.
534   if (!LexicographicallyOrderResidualBlocks(size_of_first_elimination_group,
535                                             program,
536                                             error)) {
537     return false;
538   }
539
540   return true;
541 }
542
543 bool ReorderProgramForSparseNormalCholesky(
544     const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
545     const ParameterBlockOrdering& parameter_block_ordering,
546     Program* program,
547     string* error) {
548   if (parameter_block_ordering.NumElements() != program->NumParameterBlocks()) {
549     *error = StringPrintf(
550         "The program has %d parameter blocks, but the parameter block "
551         "ordering has %d parameter blocks.",
552         program->NumParameterBlocks(),
553         parameter_block_ordering.NumElements());
554     return false;
555   }
556
557   // Compute a block sparse presentation of J'.
558   scoped_ptr<TripletSparseMatrix> tsm_block_jacobian_transpose(
559       program->CreateJacobianBlockSparsityTranspose());
560
561   vector<int> ordering(program->NumParameterBlocks(), 0);
562   vector<ParameterBlock*>& parameter_blocks =
563       *(program->mutable_parameter_blocks());
564
565   if (sparse_linear_algebra_library_type == SUITE_SPARSE) {
566     OrderingForSparseNormalCholeskyUsingSuiteSparse(
567         *tsm_block_jacobian_transpose,
568         parameter_blocks,
569         parameter_block_ordering,
570         &ordering[0]);
571   } else if (sparse_linear_algebra_library_type == CX_SPARSE) {
572     OrderingForSparseNormalCholeskyUsingCXSparse(
573         *tsm_block_jacobian_transpose,
574         &ordering[0]);
575   } else if (sparse_linear_algebra_library_type == EIGEN_SPARSE) {
576 #if EIGEN_VERSION_AT_LEAST(3, 2, 2)
577        OrderingForSparseNormalCholeskyUsingEigenSparse(
578         *tsm_block_jacobian_transpose,
579         &ordering[0]);
580 #else
581     // For Eigen versions less than 3.2.2, there is nothing to do as
582     // older versions of Eigen do not expose a method for doing
583     // symbolic analysis on pre-ordered matrices, so a block
584     // pre-ordering is a bit pointless.
585
586     return true;
587 #endif
588   }
589
590   // Apply ordering.
591   const vector<ParameterBlock*> parameter_blocks_copy(parameter_blocks);
592   for (int i = 0; i < program->NumParameterBlocks(); ++i) {
593     parameter_blocks[i] = parameter_blocks_copy[ordering[i]];
594   }
595
596   program->SetParameterOffsetsAndIndex();
597   return true;
598 }
599
600 }  // namespace internal
601 }  // namespace ceres