diff options
author | Sergey Sharybin <sergey.vfx@gmail.com> | 2014-01-13 17:11:08 +0400 |
---|---|---|
committer | Sergey Sharybin <sergey.vfx@gmail.com> | 2014-01-13 17:17:28 +0400 |
commit | a6ceb4a498dd74e855bc8a5a47120d9a1b8f879d (patch) | |
tree | 1c8afe112af66da8b2f344798a9ba3bdfa76fb3a /extern/libmv/third_party | |
parent | 4c9a3a53bd3154954bc2bcf73684b87a79332a71 (diff) |
Re-bundle new Ceres library
- Fixes some harmless issues (in cases blender doesn't use Ceres yet)
- Fixes some compilation warnings
Diffstat (limited to 'extern/libmv/third_party')
18 files changed, 488 insertions, 339 deletions
diff --git a/extern/libmv/third_party/ceres/ChangeLog b/extern/libmv/third_party/ceres/ChangeLog index 854712fe228..e80015534d6 100644 --- a/extern/libmv/third_party/ceres/ChangeLog +++ b/extern/libmv/third_party/ceres/ChangeLog @@ -1,3 +1,136 @@ +commit 80a53eebfd28bfc032cedbf7852d5c56eb1d5af5 +Author: Sameer Agarwal <sameeragarwal@google.com> +Date: Thu Jan 9 12:40:54 2014 -0800 + + Faster LBFGS. + + 1. Use column major storage for the various matrices used by + LowRankInverseHessian. Since all the operations are on columns. + + 2. Use a circular buffer to keep track of history of the LBFGS updates + so that an update does not require copying the entire history. This + makes the updates O(1) rather than O(rank). + + The implementation has been checked against the denoising code + where it gives numerically identical results. The overhead of the + LBFGS code is now near negligible as compared to the gradient evaluation. + + On a sample problem + + before 1050ms after: 630ms + + Change-Id: I537ba506ac35fc4960b304c10d923a8dea2ae031 + +commit f55780063620e7a3dcfe7e018d6488bf6a5b29da +Author: Sameer Agarwal <sameeragarwal@google.com> +Date: Wed Jan 8 10:43:31 2014 -0800 + + Reduce logging verbosity. + + When user specifies Solver::Options::logging_type = SILENT, + ensure that the minimizer does not log anything. + + Change-Id: I94e34dae504881ab36d4a66e6adb7a19a227363e + +commit 85561eee951c91e578984c6d3eecf0073acabb64 +Author: Sameer Agarwal <sameeragarwal@google.com> +Date: Tue Jan 7 22:22:14 2014 -0800 + + Use int32 for parameter block sizes. + + CostFunction now uses int32 instead of int16 + to store the size of its parameter blocks. + + This is an API breaking change. + + Change-Id: I032ea583bc7ea4b3009be25d23a3be143749c73e + +commit a7fda3317b1a97702750bea96ac3ef3d1a2afb49 +Author: Alex Stewart <alexs.mac@gmail.com> +Date: Mon Jan 6 10:25:42 2014 +0000 + + Fix typos in error messages in line search config checks. + + Change-Id: I3ae2ae58328e996598e3e32c12869d2b10109ef7 + +commit f695322eb8c5ff118f0d27f68d46d557338e5db1 +Author: Sameer Agarwal <sameeragarwal@google.com> +Date: Sat Jan 4 14:28:23 2014 -0800 + + Remove a compilation warning on windows. + + Only define NOMINMAX if it is not already defined. + + Thanks to Pierre Moulon for this fix. + + Change-Id: Ia5dc0f5ff2afe10e4c7e97a57f54297d82052b21 + +commit b811041d78d80518db153ef3030bcbdbaf80df8d +Author: Sergey Sharybin <sergey.vfx@gmail.com> +Date: Thu Jan 2 15:19:17 2014 +0600 + + Code cleanup: fix no previous declaration warnings + + Moved some internally used functions into an anonymous namespace. + + Change-Id: Ie82df61b0608abac79ccc9f7b14e7f7e04ab733d + +commit f14f6bf9b7d3fbd2cab939cf4ad615b317e93c83 +Author: Sameer Agarwal <sameeragarwal@google.com> +Date: Thu Dec 26 09:50:45 2013 -0800 + + Speed up SPARSE_NORMAL_CHOLESKY when using CX_SPARSE. + + When using sparse cholesky factorization to solve the linear + least squares problem: + + Ax = b + + There are two sources of computational complexity. + + 1. Computing H = A'A + 2. Computing the sparse Cholesky factorization of H. + + Doing 1. using CX_SPARSE is particularly expensive, as it uses + a generic cs_multiply function which computes the structure of + the matrix H everytime, reallocates memory and does not take + advantage of the fact that the matrix being computed is a symmetric + outer product. + + This change adds a custom symmetric outer product algorithm for + CompressedRowSparseMatrix. + + It has a symbolic phase, where it computes the sparsity structure + of the output matrix and a "program" which allows the actual + multiplication routine to determine exactly which entry in the + values array each term in the product contributes to. + + With these two bits of information, the outer product H = A'A + can be computed extremely fast without any reasoning about + the structure of H. + + Further gains in efficiency are made by exploiting the block + structure of A. + + With this change, SPARSE_NORMAL_CHOLESKY with CX_SPARSE as the + backend results in > 300% speedup for some problems. + + The symbolic analysis phase of the solver is a bit more expensive + now but the increased cost is made up in 3-4 iterations. + + Change-Id: I5e4a72b4d03ba41b378a2634330bc22b299c0f12 + +commit d79f886eb87cb064e19eb12c1ad3d45bbed92198 +Author: Sameer Agarwal <sameeragarwal@google.com> +Date: Mon Dec 30 07:39:10 2013 -0800 + + Refactor line search error checking code. + + Move the error checking code into its own function + so that it can be used in upcoming changes. + + Change-Id: Icf348e5a8bbe8f8b663f04fb8cfc9a2149b12f22 + commit 2b16b0080b6e673eaaf9ed478c9e971d9fcd65de Author: Sameer Agarwal <sameeragarwal@google.com> Date: Fri Dec 20 15:22:26 2013 -0800 @@ -553,131 +686,3 @@ Date: Tue Nov 5 13:10:27 2013 +0000 information as to which sub-modules are missing. Change-Id: I6eed94af49263540b8f87917b75c41b8f49658a0 - -commit 9ba0b352a282f08b1b6368a5690434407d7c81af -Author: Sameer Agarwal <sameeragarwal@google.com> -Date: Tue Nov 5 13:04:56 2013 -0800 - - Lint and other cleanups from William Rucklidge - - Change-Id: I7fb23c2db85f0f121204560b79f1966f3d584431 - -commit 69bd65ff4368ce2841519f00ff48c5284c1743a3 -Author: Alex Stewart <alexs.mac@gmail.com> -Date: Mon Nov 4 23:01:14 2013 +0000 - - Downgrading warning messages when optional deps are not found. - - - Now when find_package() is called for a dependency without the - REQUIRED or QUIET qualifiers, we emit no priority (above STATUS, but - below WARNING) messages and continue. - - Change-Id: I8cdeda7a8f6c91d45fb7f24fb366244c6c9b66e1 - -commit b0a8731fcdde31e6c37a54e8c1e1c00f853c0d5c -Author: Alex Stewart <alexs.mac@gmail.com> -Date: Mon Nov 4 20:32:40 2013 +0000 - - Removing duplicate SuiteSparse found message. - - - Also flipping ordering of variables in - find_package_handle_standard_args() so that the automatically - generated message prints the include directories, not TRUE. - - Change-Id: I2bf62eacd5c96f27152e9542b9a74651243a584e - -commit 6fed9fe0de9d1737095c24e19ad8df9735b7e572 -Author: Alex Stewart <alexs.mac@gmail.com> -Date: Mon Nov 4 18:33:05 2013 +0000 - - Fix FindPackage scripts to emit warnings, not errors if not found. - - - Previously we used message priority: SEND_ERROR when a package was - not found and find_package() was called without QUIET or REQUIRED, - which emits an error message, and prevents generation, but continues - configuration. - - The fact SEND_ERROR induces an error message was confusing for users - as it implies that something bad has happened and they cannot - continue, when in fact we were disabling the option in question - and were thus able to continue, all they had to do was re-configure. - - - This commit also reorders the search lists for includes/libraries - so that we always search user installed locations (e.g. /usr/local) - before system installed locations. Thus we will now always prefer - a user install to a system install if both are available, which is - likely to be the users desired intention. - - Change-Id: Ide84919f27d3373f31282f70c685720cd77a6723 - -commit cada337149cbc4b9e6f2bae14593b87ecf8f1a5c -Author: Alex Stewart <alexs.mac@gmail.com> -Date: Mon Nov 4 18:08:24 2013 +0000 - - Fixing CXSparse include directories statement. - - - Reported as issue #135: - https://code.google.com/p/ceres-solver/issues/detail?id=135. - - CXSPARSE_INCLUDE was the legacy include directory variable, since - the buildsystem updates we now use the CMake standard: - CXSPARSE_INCLUDE_DIRS. - - Change-Id: Iab0c2de14d524bb9e9da230bc574b5e6f09e1f31 - -commit c71085ed326239dc2d318d848ded9a99e4e3c107 -Author: Sameer Agarwal <sameeragarwal@google.com> -Date: Thu Oct 31 13:56:38 2013 -0700 - - Update to 1.8.0rc1. - - Change-Id: Iaa10fd5a20be2ef84aca0119306c44669d87cc5d - -commit 88a703f44ff0d6d5d4601584fa77f5ce853025f4 -Author: Petter Strandmark <petter.strandmark@gmail.com> -Date: Thu Oct 31 21:13:48 2013 +0100 - - Fix compilation in Visual C++ 2013. - - I had to fix the following things to make Ceres compile in 2013: - * Not link to 'm' (GNU math library). - * Excplicitly convert an std::ostream to bool. - * Include <algorithm> for std::max. - - Change-Id: I3ff65413baf8711364360d46dd71fd553fa63e72 - -commit f06b9face5bfbbc2b338aa2460bee2298a3865c5 -Author: Sameer Agarwal <sameeragarwal@google.com> -Date: Sun Oct 27 21:38:13 2013 -0700 - - Add support for multiple visibility clustering algorithms. - - The original visibility based preconditioning paper and - implementation only used the canonical views algorithm. - - This algorithm for large dense graphs can be particularly - expensive. As its worst case complexity is cubic in size - of the graph. - - Further, for many uses the SCHUR_JACOBI preconditioner - was both effective enough while being cheap. It however - suffers from a fatal flaw. If the camera parameter blocks - are split between two or more parameter blocks, e.g, - extrinsics and intrinsics. The preconditioner because - it is block diagonal will not capture the interactions - between them. - - Using CLUSTER_JACOBI or CLUSTER_TRIDIAGONAL will fix - this problem but as mentioned above this can be quite - expensive depending on the problem. - - This change extends the visibility based preconditioner - to allow for multiple clustering algorithms. And adds - a simple thresholded single linkage clustering algorithm - which allows you to construct versions of CLUSTER_JACOBI - and CLUSTER_TRIDIAGONAL preconditioners that are cheap - to construct and are more effective than SCHUR_JACOBI. - - Currently the constants controlling the threshold above - which edges are considered in the single linkage algorithm - are not exposed. This would be done in a future change. - - Change-Id: I7ddc36790943f24b19c7f08b10694ae9a822f5c9 diff --git a/extern/libmv/third_party/ceres/include/ceres/cost_function.h b/extern/libmv/third_party/ceres/include/ceres/cost_function.h index 8013e962616..722ac7732ea 100644 --- a/extern/libmv/third_party/ceres/include/ceres/cost_function.h +++ b/extern/libmv/third_party/ceres/include/ceres/cost_function.h @@ -115,7 +115,7 @@ class CostFunction { double* residuals, double** jacobians) const = 0; - const vector<int16>& parameter_block_sizes() const { + const vector<int32>& parameter_block_sizes() const { return parameter_block_sizes_; } @@ -124,7 +124,7 @@ class CostFunction { } protected: - vector<int16>* mutable_parameter_block_sizes() { + vector<int32>* mutable_parameter_block_sizes() { return ¶meter_block_sizes_; } @@ -135,7 +135,7 @@ class CostFunction { private: // Cost function signature metadata: number of inputs & their sizes, // number of outputs (residuals). - vector<int16> parameter_block_sizes_; + vector<int32> parameter_block_sizes_; int num_residuals_; CERES_DISALLOW_COPY_AND_ASSIGN(CostFunction); }; diff --git a/extern/libmv/third_party/ceres/include/ceres/cost_function_to_functor.h b/extern/libmv/third_party/ceres/include/ceres/cost_function_to_functor.h index fa1012de741..0d01f772a3b 100644 --- a/extern/libmv/third_party/ceres/include/ceres/cost_function_to_functor.h +++ b/extern/libmv/third_party/ceres/include/ceres/cost_function_to_functor.h @@ -127,7 +127,7 @@ class CostFunctionToFunctor { << N3 << ", " << N4 << ", " << N5 << ", " << N6 << ", " << N7 << ", " << N8 << ", " << N9; - const vector<int16>& parameter_block_sizes = + const vector<int32>& parameter_block_sizes = cost_function->parameter_block_sizes(); const int num_parameter_blocks = (N0 > 0) + (N1 > 0) + (N2 > 0) + (N3 > 0) + (N4 > 0) + @@ -679,7 +679,7 @@ class CostFunctionToFunctor { template <typename JetT> bool EvaluateWithJets(const JetT** inputs, JetT* output) const { const int kNumParameters = N0 + N1 + N2 + N3 + N4 + N5 + N6 + N7 + N8 + N9; - const vector<int16>& parameter_block_sizes = + const vector<int32>& parameter_block_sizes = cost_function_->parameter_block_sizes(); const int num_parameter_blocks = parameter_block_sizes.size(); const int num_residuals = cost_function_->num_residuals(); @@ -732,7 +732,7 @@ class CostFunctionToFunctor { output[i].v.setZero(); for (int j = 0; j < num_parameter_blocks; ++j) { - const int16 block_size = parameter_block_sizes[j]; + const int32 block_size = parameter_block_sizes[j]; for (int k = 0; k < parameter_block_sizes[j]; ++k) { output[i].v += jacobian_blocks[j][i * block_size + k] * inputs[j][k].v; diff --git a/extern/libmv/third_party/ceres/include/ceres/dynamic_numeric_diff_cost_function.h b/extern/libmv/third_party/ceres/include/ceres/dynamic_numeric_diff_cost_function.h index c2bfb3223cb..2b6e8260286 100644 --- a/extern/libmv/third_party/ceres/include/ceres/dynamic_numeric_diff_cost_function.h +++ b/extern/libmv/third_party/ceres/include/ceres/dynamic_numeric_diff_cost_function.h @@ -104,7 +104,7 @@ class DynamicNumericDiffCostFunction : public CostFunction { << "You must call DynamicNumericDiffCostFunction::SetNumResiduals() " << "before DynamicNumericDiffCostFunction::Evaluate()."; - const vector<int16>& block_sizes = parameter_block_sizes(); + const vector<int32>& block_sizes = parameter_block_sizes(); CHECK(!block_sizes.empty()) << "You must call DynamicNumericDiffCostFunction::AddParameterBlock() " << "before DynamicNumericDiffCostFunction::Evaluate()."; diff --git a/extern/libmv/third_party/ceres/include/ceres/gradient_checker.h b/extern/libmv/third_party/ceres/include/ceres/gradient_checker.h index 3ec8056a02d..79ebae5f13e 100644 --- a/extern/libmv/third_party/ceres/include/ceres/gradient_checker.h +++ b/extern/libmv/third_party/ceres/include/ceres/gradient_checker.h @@ -119,7 +119,7 @@ class GradientChecker { // Do a consistency check between the term and the template parameters. CHECK_EQ(M, term->num_residuals()); const int num_residuals = M; - const vector<int16>& block_sizes = term->parameter_block_sizes(); + const vector<int32>& block_sizes = term->parameter_block_sizes(); const int num_blocks = block_sizes.size(); CHECK_LE(num_blocks, 5) << "Unable to test functions that take more " diff --git a/extern/libmv/third_party/ceres/internal/ceres/block_structure.h b/extern/libmv/third_party/ceres/internal/ceres/block_structure.h index f509067d216..2abb36a0af5 100644 --- a/extern/libmv/third_party/ceres/internal/ceres/block_structure.h +++ b/extern/libmv/third_party/ceres/internal/ceres/block_structure.h @@ -47,7 +47,7 @@ namespace internal { class BlockStructureProto; -typedef int16 BlockSize; +typedef int32 BlockSize; struct Block { Block() : size(-1), position(-1) {} diff --git a/extern/libmv/third_party/ceres/internal/ceres/compressed_row_sparse_matrix.cc b/extern/libmv/third_party/ceres/internal/ceres/compressed_row_sparse_matrix.cc index 439270ed3bc..34f31ab6c7f 100644 --- a/extern/libmv/third_party/ceres/internal/ceres/compressed_row_sparse_matrix.cc +++ b/extern/libmv/third_party/ceres/internal/ceres/compressed_row_sparse_matrix.cc @@ -125,7 +125,7 @@ CompressedRowSparseMatrix::CompressedRowSparseMatrix( // Find the cumulative sum of the row counts. for (int i = 1; i < num_rows_ + 1; ++i) { - rows_[i] += rows_[i-1]; + rows_[i] += rows_[i - 1]; } CHECK_EQ(num_nonzeros(), m.num_nonzeros()); @@ -217,8 +217,8 @@ void CompressedRowSparseMatrix::DeleteRows(int delta_rows) { num_rows_ -= delta_rows; rows_.resize(num_rows_ + 1); - // Walk the list of row blocks untill we reach the new number of - // rows and then drop the rest of the row blocks. + // Walk the list of row blocks until we reach the new number of rows + // and the drop the rest of the row blocks. int num_row_blocks = 0; int num_rows = 0; while (num_row_blocks < row_blocks_.size() && num_rows < num_rows_) { @@ -380,6 +380,155 @@ CompressedRowSparseMatrix* CompressedRowSparseMatrix::Transpose() const { return transpose; } +namespace { +// A ProductTerm is a term in the outer product of a matrix with +// itself. +struct ProductTerm { + ProductTerm(const int row, const int col, const int index) + : row(row), col(col), index(index) { + } + + bool operator<(const ProductTerm& right) const { + if (row == right.row) { + if (col == right.col) { + return index < right.index; + } + return col < right.col; + } + return row < right.row; + } + + int row; + int col; + int index; +}; + +CompressedRowSparseMatrix* +CompressAndFillProgram(const int num_rows, + const int num_cols, + const vector<ProductTerm>& product, + vector<int>* program) { + CHECK_GT(product.size(), 0); + + // Count the number of unique product term, which in turn is the + // number of non-zeros in the outer product. + int num_nonzeros = 1; + for (int i = 1; i < product.size(); ++i) { + if (product[i].row != product[i - 1].row || + product[i].col != product[i - 1].col) { + ++num_nonzeros; + } + } + + CompressedRowSparseMatrix* matrix = + new CompressedRowSparseMatrix(num_rows, num_cols, num_nonzeros); + + int* crsm_rows = matrix->mutable_rows(); + std::fill(crsm_rows, crsm_rows + num_rows + 1, 0); + int* crsm_cols = matrix->mutable_cols(); + std::fill(crsm_cols, crsm_cols + num_nonzeros, 0); + + CHECK_NOTNULL(program)->clear(); + program->resize(product.size()); + + // Iterate over the sorted product terms. This means each row is + // filled one at a time, and we are able to assign a position in the + // values array to each term. + // + // If terms repeat, i.e., they contribute to the same entry in the + // result matrix), then they do not affect the sparsity structure of + // the result matrix. + int nnz = 0; + crsm_cols[0] = product[0].col; + crsm_rows[product[0].row + 1]++; + (*program)[product[0].index] = nnz; + for (int i = 1; i < product.size(); ++i) { + const ProductTerm& previous = product[i - 1]; + const ProductTerm& current = product[i]; + + // Sparsity structure is updated only if the term is not a repeat. + if (previous.row != current.row || previous.col != current.col) { + crsm_cols[++nnz] = current.col; + crsm_rows[current.row + 1]++; + } + + // All terms get assigned the position in the values array where + // their value is accumulated. + (*program)[current.index] = nnz; + } + + for (int i = 1; i < num_rows + 1; ++i) { + crsm_rows[i] += crsm_rows[i - 1]; + } + + return matrix; +} + +} // namespace + +CompressedRowSparseMatrix* +CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram( + const CompressedRowSparseMatrix& m, + vector<int>* program) { + CHECK_NOTNULL(program)->clear(); + CHECK_GT(m.num_nonzeros(), 0) << "Congratulations, " + << "you found a bug in Ceres. Please report it."; + + vector<ProductTerm> product; + const vector<int>& row_blocks = m.row_blocks(); + int row_block_begin = 0; + // Iterate over row blocks + for (int row_block = 0; row_block < row_blocks.size(); ++row_block) { + const int row_block_end = row_block_begin + row_blocks[row_block]; + // Compute the outer product terms for just one row per row block. + const int r = row_block_begin; + // Compute the lower triangular part of the product. + for (int idx1 = m.rows()[r]; idx1 < m.rows()[r + 1]; ++idx1) { + for (int idx2 = m.rows()[r]; idx2 <= idx1; ++idx2) { + product.push_back(ProductTerm(m.cols()[idx1], m.cols()[idx2], product.size())); + } + } + row_block_begin = row_block_end; + } + CHECK_EQ(row_block_begin, m.num_rows()); + sort(product.begin(), product.end()); + return CompressAndFillProgram(m.num_cols(), m.num_cols(), product, program); +} + +void CompressedRowSparseMatrix::ComputeOuterProduct( + const CompressedRowSparseMatrix& m, + const vector<int>& program, + CompressedRowSparseMatrix* result) { + result->SetZero(); + double* values = result->mutable_values(); + const vector<int>& row_blocks = m.row_blocks(); + + int cursor = 0; + int row_block_begin = 0; + const double* m_values = m.values(); + const int* m_rows = m.rows(); + // Iterate over row blocks. + for (int row_block = 0; row_block < row_blocks.size(); ++row_block) { + const int row_block_end = row_block_begin + row_blocks[row_block]; + const int saved_cursor = cursor; + for (int r = row_block_begin; r < row_block_end; ++r) { + // Reuse the program segment for each row in this row block. + cursor = saved_cursor; + const int row_begin = m_rows[r]; + const int row_end = m_rows[r + 1]; + for (int idx1 = row_begin; idx1 < row_end; ++idx1) { + const double v1 = m_values[idx1]; + for (int idx2 = row_begin; idx2 <= idx1; ++idx2, ++cursor) { + values[program[cursor]] += v1 * m_values[idx2]; + } + } + } + row_block_begin = row_block_end; + } + + CHECK_EQ(row_block_begin, m.num_rows()); + CHECK_EQ(cursor, program.size()); +} } // namespace internal } // namespace ceres diff --git a/extern/libmv/third_party/ceres/internal/ceres/compressed_row_sparse_matrix.h b/extern/libmv/third_party/ceres/internal/ceres/compressed_row_sparse_matrix.h index c5721eb888a..06b86896f23 100644 --- a/extern/libmv/third_party/ceres/internal/ceres/compressed_row_sparse_matrix.h +++ b/extern/libmv/third_party/ceres/internal/ceres/compressed_row_sparse_matrix.h @@ -128,6 +128,32 @@ class CompressedRowSparseMatrix : public SparseMatrix { const double* diagonal, const vector<int>& blocks); + // Compute the sparsity structure of the product m.transpose() * m + // and create a CompressedRowSparseMatrix corresponding to it. + // + // Also compute a "program" vector, which for every term in the + // outer product points to the entry in the values array of the + // result matrix where it should be accumulated. + // + // This program is used by the ComputeOuterProduct function below to + // compute the outer product. + // + // Since the entries of the program are the same for rows with the + // same sparsity structure, the program only stores the result for + // one row per row block. The ComputeOuterProduct function reuses + // this information for each row in the row block. + static CompressedRowSparseMatrix* CreateOuterProductMatrixAndProgram( + const CompressedRowSparseMatrix& m, + vector<int>* program); + + // Compute the values array for the expression m.transpose() * m, + // where the matrix used to store the result and a program have been + // created using the CreateOuterProductMatrixAndProgram function + // above. + static void ComputeOuterProduct(const CompressedRowSparseMatrix& m, + const vector<int>& program, + CompressedRowSparseMatrix* result); + private: int num_rows_; int num_cols_; diff --git a/extern/libmv/third_party/ceres/internal/ceres/gradient_checking_cost_function.cc b/extern/libmv/third_party/ceres/internal/ceres/gradient_checking_cost_function.cc index 550301359ad..bca22e6de03 100644 --- a/extern/libmv/third_party/ceres/internal/ceres/gradient_checking_cost_function.cc +++ b/extern/libmv/third_party/ceres/internal/ceres/gradient_checking_cost_function.cc @@ -93,7 +93,7 @@ class GradientCheckingCostFunction : public CostFunction { DO_NOT_TAKE_OWNERSHIP, relative_step_size); - const vector<int16>& parameter_block_sizes = + const vector<int32>& parameter_block_sizes = function->parameter_block_sizes(); for (int i = 0; i < parameter_block_sizes.size(); ++i) { finite_diff_cost_function->AddParameterBlock(parameter_block_sizes[i]); @@ -117,7 +117,7 @@ class GradientCheckingCostFunction : public CostFunction { int num_residuals = function_->num_residuals(); // Make space for the jacobians of the two methods. - const vector<int16>& block_sizes = function_->parameter_block_sizes(); + const vector<int32>& block_sizes = function_->parameter_block_sizes(); vector<Matrix> term_jacobians(block_sizes.size()); vector<Matrix> finite_difference_jacobians(block_sizes.size()); vector<double*> term_jacobian_pointers(block_sizes.size()); diff --git a/extern/libmv/third_party/ceres/internal/ceres/low_rank_inverse_hessian.cc b/extern/libmv/third_party/ceres/internal/ceres/low_rank_inverse_hessian.cc index 9aeafecfb36..4816e3c20f0 100644 --- a/extern/libmv/third_party/ceres/internal/ceres/low_rank_inverse_hessian.cc +++ b/extern/libmv/third_party/ceres/internal/ceres/low_rank_inverse_hessian.cc @@ -28,6 +28,8 @@ // // Author: sameeragarwal@google.com (Sameer Agarwal) +#include <list> + #include "ceres/internal/eigen.h" #include "ceres/low_rank_inverse_hessian.h" #include "glog/logging.h" @@ -77,7 +79,6 @@ LowRankInverseHessian::LowRankInverseHessian( : num_parameters_(num_parameters), max_num_corrections_(max_num_corrections), use_approximate_eigenvalue_scaling_(use_approximate_eigenvalue_scaling), - num_corrections_(0), approximate_eigenvalue_scale_(1.0), delta_x_history_(num_parameters, max_num_corrections), delta_gradient_history_(num_parameters, max_num_corrections), @@ -96,29 +97,20 @@ bool LowRankInverseHessian::Update(const Vector& delta_x, return false; } - if (num_corrections_ == max_num_corrections_) { - // TODO(sameeragarwal): This can be done more efficiently using - // a circular buffer/indexing scheme, but for simplicity we will - // do the expensive copy for now. - delta_x_history_.block(0, 0, num_parameters_, max_num_corrections_ - 1) = - delta_x_history_ - .block(0, 1, num_parameters_, max_num_corrections_ - 1); - - delta_gradient_history_ - .block(0, 0, num_parameters_, max_num_corrections_ - 1) = - delta_gradient_history_ - .block(0, 1, num_parameters_, max_num_corrections_ - 1); - - delta_x_dot_delta_gradient_.head(num_corrections_ - 1) = - delta_x_dot_delta_gradient_.tail(num_corrections_ - 1); - } else { - ++num_corrections_; + + int next = indices_.size(); + // Once the size of the list reaches max_num_corrections_, simulate + // a circular buffer by removing the first element of the list and + // making it the next position where the LBFGS history is stored. + if (next == max_num_corrections_) { + next = indices_.front(); + indices_.pop_front(); } - delta_x_history_.col(num_corrections_ - 1) = delta_x; - delta_gradient_history_.col(num_corrections_ - 1) = delta_gradient; - delta_x_dot_delta_gradient_(num_corrections_ - 1) = - delta_x_dot_delta_gradient; + indices_.push_back(next); + delta_x_history_.col(next) = delta_x; + delta_gradient_history_.col(next) = delta_gradient; + delta_x_dot_delta_gradient_(next) = delta_x_dot_delta_gradient; approximate_eigenvalue_scale_ = delta_x_dot_delta_gradient / delta_gradient.squaredNorm(); return true; @@ -131,12 +123,16 @@ void LowRankInverseHessian::RightMultiply(const double* x_ptr, search_direction = gradient; - Vector alpha(num_corrections_); + const int num_corrections = indices_.size(); + Vector alpha(num_corrections); - for (int i = num_corrections_ - 1; i >= 0; --i) { - alpha(i) = delta_x_history_.col(i).dot(search_direction) / - delta_x_dot_delta_gradient_(i); - search_direction -= alpha(i) * delta_gradient_history_.col(i); + for (std::list<int>::const_reverse_iterator it = indices_.rbegin(); + it != indices_.rend(); + ++it) { + const double alpha_i = delta_x_history_.col(*it).dot(search_direction) / + delta_x_dot_delta_gradient_(*it); + search_direction -= alpha_i * delta_gradient_history_.col(*it); + alpha(*it) = alpha_i; } if (use_approximate_eigenvalue_scaling_) { @@ -177,10 +173,12 @@ void LowRankInverseHessian::RightMultiply(const double* x_ptr, << "approximation."; } - for (int i = 0; i < num_corrections_; ++i) { - const double beta = delta_gradient_history_.col(i).dot(search_direction) / - delta_x_dot_delta_gradient_(i); - search_direction += delta_x_history_.col(i) * (alpha(i) - beta); + for (std::list<int>::const_iterator it = indices_.begin(); + it != indices_.end(); + ++it) { + const double beta = delta_gradient_history_.col(*it).dot(search_direction) / + delta_x_dot_delta_gradient_(*it); + search_direction += delta_x_history_.col(*it) * (alpha(*it) - beta); } } diff --git a/extern/libmv/third_party/ceres/internal/ceres/low_rank_inverse_hessian.h b/extern/libmv/third_party/ceres/internal/ceres/low_rank_inverse_hessian.h index 7d293d09422..19ab7606aa0 100644 --- a/extern/libmv/third_party/ceres/internal/ceres/low_rank_inverse_hessian.h +++ b/extern/libmv/third_party/ceres/internal/ceres/low_rank_inverse_hessian.h @@ -34,6 +34,8 @@ #ifndef CERES_INTERNAL_LOW_RANK_INVERSE_HESSIAN_H_ #define CERES_INTERNAL_LOW_RANK_INVERSE_HESSIAN_H_ +#include <list> + #include "ceres/internal/eigen.h" #include "ceres/linear_operator.h" @@ -93,11 +95,11 @@ class LowRankInverseHessian : public LinearOperator { const int num_parameters_; const int max_num_corrections_; const bool use_approximate_eigenvalue_scaling_; - int num_corrections_; double approximate_eigenvalue_scale_; - Matrix delta_x_history_; - Matrix delta_gradient_history_; + ColMajorMatrix delta_x_history_; + ColMajorMatrix delta_gradient_history_; Vector delta_x_dot_delta_gradient_; + std::list<int> indices_; }; } // namespace internal diff --git a/extern/libmv/third_party/ceres/internal/ceres/minimizer.h b/extern/libmv/third_party/ceres/internal/ceres/minimizer.h index 3d9da997d24..ee77726d2b4 100644 --- a/extern/libmv/third_party/ceres/internal/ceres/minimizer.h +++ b/extern/libmv/third_party/ceres/internal/ceres/minimizer.h @@ -107,7 +107,7 @@ class Minimizer { options.line_search_sufficient_curvature_decrease; max_line_search_step_expansion = options.max_line_search_step_expansion; - is_silent = false; + is_silent = (options.logging_type == SILENT); evaluator = NULL; trust_region_strategy = NULL; jacobian = NULL; diff --git a/extern/libmv/third_party/ceres/internal/ceres/mutex.h b/extern/libmv/third_party/ceres/internal/ceres/mutex.h index 0c48ed352b5..564c39f1b12 100644 --- a/extern/libmv/third_party/ceres/internal/ceres/mutex.h +++ b/extern/libmv/third_party/ceres/internal/ceres/mutex.h @@ -112,7 +112,9 @@ // To avoid macro definition of ERROR. # define NOGDI // To avoid macro definition of min/max. -# define NOMINMAX +# ifndef NOMINMAX +# define NOMINMAX +# endif # include <windows.h> typedef CRITICAL_SECTION MutexType; #elif defined(CERES_HAVE_PTHREAD) && defined(CERES_HAVE_RWLOCK) diff --git a/extern/libmv/third_party/ceres/internal/ceres/polynomial.cc b/extern/libmv/third_party/ceres/internal/ceres/polynomial.cc index c4b0243d4f6..75f43de409a 100644 --- a/extern/libmv/third_party/ceres/internal/ceres/polynomial.cc +++ b/extern/libmv/third_party/ceres/internal/ceres/polynomial.cc @@ -120,7 +120,6 @@ Vector RemoveLeadingZeros(const Vector& polynomial_in) { } return polynomial_in.tail(polynomial_in.size() - i); } -} // namespace void FindLinearPolynomialRoots(const Vector& polynomial, Vector* real, @@ -178,6 +177,7 @@ void FindQuadraticPolynomialRoots(const Vector& polynomial, (*imaginary)(1) = -sqrt_D / (2.0 * a); } } +} // namespace bool FindPolynomialRoots(const Vector& polynomial_in, Vector* real, diff --git a/extern/libmv/third_party/ceres/internal/ceres/problem_impl.cc b/extern/libmv/third_party/ceres/internal/ceres/problem_impl.cc index 09eed9e614c..37cd351bf2a 100644 --- a/extern/libmv/third_party/ceres/internal/ceres/problem_impl.cc +++ b/extern/libmv/third_party/ceres/internal/ceres/problem_impl.cc @@ -224,7 +224,7 @@ ResidualBlock* ProblemImpl::AddResidualBlock( cost_function->parameter_block_sizes().size()); // Check the sizes match. - const vector<int16>& parameter_block_sizes = + const vector<int32>& parameter_block_sizes = cost_function->parameter_block_sizes(); if (!options_.disable_all_safety_checks) { diff --git a/extern/libmv/third_party/ceres/internal/ceres/solver_impl.cc b/extern/libmv/third_party/ceres/internal/ceres/solver_impl.cc index 6c579a6098d..cf0e2515ab6 100644 --- a/extern/libmv/third_party/ceres/internal/ceres/solver_impl.cc +++ b/extern/libmv/third_party/ceres/internal/ceres/solver_impl.cc @@ -484,7 +484,7 @@ void SolverImpl::TrustRegionSolve(const Solver::Options& original_options, "Terminating: Function tolerance reached. " "No non-constant parameter blocks found."; summary->termination_type = CONVERGENCE; - VLOG(1) << summary->message; + VLOG_IF(1, options.logging_type != SILENT) << summary->message; summary->initial_cost = summary->fixed_cost; summary->final_cost = summary->fixed_cost; @@ -620,107 +620,106 @@ void SolverImpl::TrustRegionSolve(const Solver::Options& original_options, #ifndef CERES_NO_LINE_SEARCH_MINIMIZER -void SolverImpl::LineSearchSolve(const Solver::Options& original_options, - ProblemImpl* original_problem_impl, - Solver::Summary* summary) { - double solver_start_time = WallTimeInSeconds(); - - Program* original_program = original_problem_impl->mutable_program(); - ProblemImpl* problem_impl = original_problem_impl; - - // Reset the summary object to its default values. - *CHECK_NOTNULL(summary) = Solver::Summary(); - - summary->minimizer_type = LINE_SEARCH; - SummarizeGivenProgram(*original_program, summary); - - summary->line_search_direction_type = - original_options.line_search_direction_type; - summary->max_lbfgs_rank = original_options.max_lbfgs_rank; - summary->line_search_type = original_options.line_search_type; - summary->line_search_interpolation_type = - original_options.line_search_interpolation_type; - summary->nonlinear_conjugate_gradient_type = - original_options.nonlinear_conjugate_gradient_type; - +bool LineSearchOptionsAreValid(const Solver::Options& options, + string* message) { // Validate values for configuration parameters supplied by user. - if ((original_options.line_search_direction_type == ceres::BFGS || - original_options.line_search_direction_type == ceres::LBFGS) && - original_options.line_search_type != ceres::WOLFE) { - summary->message = + if ((options.line_search_direction_type == ceres::BFGS || + options.line_search_direction_type == ceres::LBFGS) && + options.line_search_type != ceres::WOLFE) { + *message = string("Invalid configuration: require line_search_type == " "ceres::WOLFE when using (L)BFGS to ensure that underlying " "assumptions are guaranteed to be satisfied."); - LOG(ERROR) << summary->message; - return; + return false; } - if (original_options.max_lbfgs_rank <= 0) { - summary->message = + if (options.max_lbfgs_rank <= 0) { + *message = string("Invalid configuration: require max_lbfgs_rank > 0"); - LOG(ERROR) << summary->message; - return; + return false; } - if (original_options.min_line_search_step_size <= 0.0) { - summary->message = - "Invalid configuration: min_line_search_step_size <= 0.0."; - LOG(ERROR) << summary->message; - return; + if (options.min_line_search_step_size <= 0.0) { + *message = + "Invalid configuration: require min_line_search_step_size > 0.0."; + return false; } - if (original_options.line_search_sufficient_function_decrease <= 0.0) { - summary->message = + if (options.line_search_sufficient_function_decrease <= 0.0) { + *message = string("Invalid configuration: require ") + - string("line_search_sufficient_function_decrease <= 0.0."); - LOG(ERROR) << summary->message; - return; + string("line_search_sufficient_function_decrease > 0.0."); + return false; } - if (original_options.max_line_search_step_contraction <= 0.0 || - original_options.max_line_search_step_contraction >= 1.0) { - summary->message = string("Invalid configuration: require ") + + if (options.max_line_search_step_contraction <= 0.0 || + options.max_line_search_step_contraction >= 1.0) { + *message = string("Invalid configuration: require ") + string("0.0 < max_line_search_step_contraction < 1.0."); - LOG(ERROR) << summary->message; - return; + return false; } - if (original_options.min_line_search_step_contraction <= - original_options.max_line_search_step_contraction || - original_options.min_line_search_step_contraction > 1.0) { - summary->message = string("Invalid configuration: require ") + + if (options.min_line_search_step_contraction <= + options.max_line_search_step_contraction || + options.min_line_search_step_contraction > 1.0) { + *message = string("Invalid configuration: require ") + string("max_line_search_step_contraction < ") + string("min_line_search_step_contraction <= 1.0."); - LOG(ERROR) << summary->message; - return; + return false; } // Warn user if they have requested BISECTION interpolation, but constraints // on max/min step size change during line search prevent bisection scaling // from occurring. Warn only, as this is likely a user mistake, but one which // does not prevent us from continuing. LOG_IF(WARNING, - (original_options.line_search_interpolation_type == ceres::BISECTION && - (original_options.max_line_search_step_contraction > 0.5 || - original_options.min_line_search_step_contraction < 0.5))) + (options.line_search_interpolation_type == ceres::BISECTION && + (options.max_line_search_step_contraction > 0.5 || + options.min_line_search_step_contraction < 0.5))) << "Line search interpolation type is BISECTION, but specified " << "max_line_search_step_contraction: " - << original_options.max_line_search_step_contraction << ", and " + << options.max_line_search_step_contraction << ", and " << "min_line_search_step_contraction: " - << original_options.min_line_search_step_contraction + << options.min_line_search_step_contraction << ", prevent bisection (0.5) scaling, continuing with solve regardless."; - if (original_options.max_num_line_search_step_size_iterations <= 0) { - summary->message = string("Invalid configuration: require ") + + if (options.max_num_line_search_step_size_iterations <= 0) { + *message = string("Invalid configuration: require ") + string("max_num_line_search_step_size_iterations > 0."); - LOG(ERROR) << summary->message; - return; + return false; } - if (original_options.line_search_sufficient_curvature_decrease <= - original_options.line_search_sufficient_function_decrease || - original_options.line_search_sufficient_curvature_decrease > 1.0) { - summary->message = string("Invalid configuration: require ") + + if (options.line_search_sufficient_curvature_decrease <= + options.line_search_sufficient_function_decrease || + options.line_search_sufficient_curvature_decrease > 1.0) { + *message = string("Invalid configuration: require ") + string("line_search_sufficient_function_decrease < ") + string("line_search_sufficient_curvature_decrease < 1.0."); - LOG(ERROR) << summary->message; - return; + return false; } - if (original_options.max_line_search_step_expansion <= 1.0) { - summary->message = string("Invalid configuration: require ") + + if (options.max_line_search_step_expansion <= 1.0) { + *message = string("Invalid configuration: require ") + string("max_line_search_step_expansion > 1.0."); + return false; + } + return true; +} + +void SolverImpl::LineSearchSolve(const Solver::Options& original_options, + ProblemImpl* original_problem_impl, + Solver::Summary* summary) { + double solver_start_time = WallTimeInSeconds(); + + Program* original_program = original_problem_impl->mutable_program(); + ProblemImpl* problem_impl = original_problem_impl; + + // Reset the summary object to its default values. + *CHECK_NOTNULL(summary) = Solver::Summary(); + + SummarizeGivenProgram(*original_program, summary); + summary->minimizer_type = LINE_SEARCH; + summary->line_search_direction_type = + original_options.line_search_direction_type; + summary->max_lbfgs_rank = original_options.max_lbfgs_rank; + summary->line_search_type = original_options.line_search_type; + summary->line_search_interpolation_type = + original_options.line_search_interpolation_type; + summary->nonlinear_conjugate_gradient_type = + original_options.nonlinear_conjugate_gradient_type; + + if (!LineSearchOptionsAreValid(original_options, &summary->message)) { LOG(ERROR) << summary->message; return; } @@ -805,7 +804,7 @@ void SolverImpl::LineSearchSolve(const Solver::Options& original_options, "Terminating: Function tolerance reached. " "No non-constant parameter blocks found."; summary->termination_type = CONVERGENCE; - VLOG(1) << summary->message; + VLOG_IF(1, options.logging_type != SILENT) << summary->message; const double post_process_start_time = WallTimeInSeconds(); SetSummaryFinalCost(summary); diff --git a/extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.cc b/extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.cc index b2dba440e10..880adf7ebca 100644 --- a/extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.cc +++ b/extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.cc @@ -77,27 +77,50 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImpl( const double* b, const LinearSolver::PerSolveOptions& per_solve_options, double * x) { + + const int num_cols = A->num_cols(); + VectorRef(x, num_cols).setZero(); + A->LeftMultiply(b, x); + + if (per_solve_options.D != NULL) { + // Temporarily append a diagonal block to the A matrix, but undo + // it before returning the matrix to the user. + scoped_ptr<CompressedRowSparseMatrix> regularizer; + if (A->col_blocks().size() > 0) { + regularizer.reset(CompressedRowSparseMatrix::CreateBlockDiagonalMatrix( + per_solve_options.D, A->col_blocks())); + } else { + regularizer.reset(new CompressedRowSparseMatrix( + per_solve_options.D, num_cols)); + } + A->AppendRows(*regularizer); + } + + LinearSolver::Summary summary; switch (options_.sparse_linear_algebra_library_type) { case SUITE_SPARSE: - return SolveImplUsingSuiteSparse(A, b, per_solve_options, x); + summary = SolveImplUsingSuiteSparse(A, per_solve_options, x); + break; case CX_SPARSE: - return SolveImplUsingCXSparse(A, b, per_solve_options, x); + summary = SolveImplUsingCXSparse(A, per_solve_options, x); + break; default: LOG(FATAL) << "Unknown sparse linear algebra library : " << options_.sparse_linear_algebra_library_type; } - LOG(FATAL) << "Unknown sparse linear algebra library : " - << options_.sparse_linear_algebra_library_type; - return LinearSolver::Summary(); + if (per_solve_options.D != NULL) { + A->DeleteRows(num_cols); + } + + return summary; } #ifndef CERES_NO_CXSPARSE LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( CompressedRowSparseMatrix* A, - const double* b, const LinearSolver::PerSolveOptions& per_solve_options, - double * x) { + double * rhs_and_solution) { EventLogger event_logger("SparseNormalCholeskySolver::CXSparse::Solve"); LinearSolver::Summary summary; @@ -105,29 +128,6 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( summary.termination_type = LINEAR_SOLVER_SUCCESS; summary.message = "Success."; - const int num_cols = A->num_cols(); - Vector Atb = Vector::Zero(num_cols); - A->LeftMultiply(b, Atb.data()); - - if (per_solve_options.D != NULL) { - // Temporarily append a diagonal block to the A matrix, but undo - // it before returning the matrix to the user. - scoped_ptr<CompressedRowSparseMatrix> regularizer; - if (A->col_blocks().size() > 0) { - regularizer.reset(CompressedRowSparseMatrix::CreateBlockDiagonalMatrix( - per_solve_options.D, A->col_blocks())); - } else { - regularizer.reset(new CompressedRowSparseMatrix( - per_solve_options.D, num_cols)); - } - A->AppendRows(*regularizer); - } - - VectorRef(x, num_cols).setZero(); - - // Wrap the augmented Jacobian in a compressed sparse column matrix. - cs_di At = cxsparse_.CreateSparseMatrixTransposeView(A); - // Compute the normal equations. J'J delta = J'f and solve them // using a sparse Cholesky factorization. Notice that when compared // to SuiteSparse we have to explicitly compute the transpose of Jt, @@ -135,13 +135,18 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( // factorized. CHOLMOD/SuiteSparse on the other hand can just work // off of Jt to compute the Cholesky factorization of the normal // equations. - cs_di* A2 = cxsparse_.TransposeMatrix(&At); - cs_di* AtA = cxsparse_.MatrixMatrixMultiply(&At, A2); - - cxsparse_.Free(A2); - if (per_solve_options.D != NULL) { - A->DeleteRows(num_cols); + if (outer_product_.get() == NULL) { + outer_product_.reset( + CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram( + *A, &pattern_)); } + + CompressedRowSparseMatrix::ComputeOuterProduct( + *A, pattern_, outer_product_.get()); + cs_di AtA_view = + cxsparse_.CreateSparseMatrixTransposeView(outer_product_.get()); + cs_di* AtA = &AtA_view; + event_logger.AddEvent("Setup"); // Compute symbolic factorization if not available. @@ -160,23 +165,17 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( summary.termination_type = LINEAR_SOLVER_FATAL_ERROR; summary.message = "CXSparse failure. Unable to find symbolic factorization."; - } else if (cxsparse_.SolveCholesky(AtA, cxsparse_factor_, Atb.data())) { - VectorRef(x, Atb.rows()) = Atb; - } else { + } else if (!cxsparse_.SolveCholesky(AtA, cxsparse_factor_, rhs_and_solution)) { summary.termination_type = LINEAR_SOLVER_FAILURE; } event_logger.AddEvent("Solve"); - - cxsparse_.Free(AtA); - event_logger.AddEvent("Teardown"); return summary; } #else LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( CompressedRowSparseMatrix* A, - const double* b, const LinearSolver::PerSolveOptions& per_solve_options, - double * x) { + double * rhs_and_solution) { LOG(FATAL) << "No CXSparse support in Ceres."; // Unreachable but MSVC does not know this. @@ -187,9 +186,8 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( #ifndef CERES_NO_SUITESPARSE LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( CompressedRowSparseMatrix* A, - const double* b, const LinearSolver::PerSolveOptions& per_solve_options, - double * x) { + double * rhs_and_solution) { EventLogger event_logger("SparseNormalCholeskySolver::SuiteSparse::Solve"); LinearSolver::Summary summary; summary.termination_type = LINEAR_SOLVER_SUCCESS; @@ -197,24 +195,6 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( summary.message = "Success."; const int num_cols = A->num_cols(); - Vector Atb = Vector::Zero(num_cols); - A->LeftMultiply(b, Atb.data()); - - if (per_solve_options.D != NULL) { - // Temporarily append a diagonal block to the A matrix, but undo - // it before returning the matrix to the user. - scoped_ptr<CompressedRowSparseMatrix> regularizer; - if (A->col_blocks().size() > 0) { - regularizer.reset(CompressedRowSparseMatrix::CreateBlockDiagonalMatrix( - per_solve_options.D, A->col_blocks())); - } else { - regularizer.reset(new CompressedRowSparseMatrix( - per_solve_options.D, num_cols)); - } - A->AppendRows(*regularizer); - } - - VectorRef(x, num_cols).setZero(); cholmod_sparse lhs = ss_.CreateSparseMatrixTransposeView(A); event_logger.AddEvent("Setup"); @@ -231,33 +211,23 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( event_logger.AddEvent("Analysis"); if (factor_ == NULL) { - if (per_solve_options.D != NULL) { - A->DeleteRows(num_cols); - } summary.termination_type = LINEAR_SOLVER_FATAL_ERROR; return summary; } summary.termination_type = ss_.Cholesky(&lhs, factor_, &summary.message); if (summary.termination_type != LINEAR_SOLVER_SUCCESS) { - if (per_solve_options.D != NULL) { - A->DeleteRows(num_cols); - } return summary; } - cholmod_dense* rhs = ss_.CreateDenseVector(Atb.data(), num_cols, num_cols); - cholmod_dense* sol = ss_.Solve(factor_, rhs, &summary.message); + cholmod_dense* rhs = ss_.CreateDenseVector(rhs_and_solution, num_cols, num_cols); + cholmod_dense* solution = ss_.Solve(factor_, rhs, &summary.message); event_logger.AddEvent("Solve"); ss_.Free(rhs); - if (per_solve_options.D != NULL) { - A->DeleteRows(num_cols); - } - - if (sol != NULL) { - memcpy(x, sol->x, num_cols * sizeof(*x)); - ss_.Free(sol); + if (solution != NULL) { + memcpy(rhs_and_solution, solution->x, num_cols * sizeof(*rhs_and_solution)); + ss_.Free(solution); } else { summary.termination_type = LINEAR_SOLVER_FAILURE; } @@ -268,9 +238,8 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( #else LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( CompressedRowSparseMatrix* A, - const double* b, const LinearSolver::PerSolveOptions& per_solve_options, - double * x) { + double * rhs_and_solution) { LOG(FATAL) << "No SuiteSparse support in Ceres."; // Unreachable but MSVC does not know this. diff --git a/extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.h b/extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.h index 61111b41b49..adca08a507f 100644 --- a/extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.h +++ b/extern/libmv/third_party/ceres/internal/ceres/sparse_normal_cholesky_solver.h @@ -62,16 +62,14 @@ class SparseNormalCholeskySolver : public CompressedRowSparseMatrixSolver { LinearSolver::Summary SolveImplUsingSuiteSparse( CompressedRowSparseMatrix* A, - const double* b, const LinearSolver::PerSolveOptions& options, - double* x); + double* rhs_and_solution); // Crashes if CSparse is not installed. LinearSolver::Summary SolveImplUsingCXSparse( CompressedRowSparseMatrix* A, - const double* b, const LinearSolver::PerSolveOptions& options, - double* x); + double* rhs_and_solution); SuiteSparse ss_; // Cached factorization @@ -80,7 +78,8 @@ class SparseNormalCholeskySolver : public CompressedRowSparseMatrixSolver { CXSparse cxsparse_; // Cached factorization cs_dis* cxsparse_factor_; - + scoped_ptr<CompressedRowSparseMatrix> outer_product_; + vector<int> pattern_; const LinearSolver::Options options_; CERES_DISALLOW_COPY_AND_ASSIGN(SparseNormalCholeskySolver); }; |