1 // Ceres Solver - A fast non-linear least squares minimizer
2 // Copyright 2017 Google Inc. All rights reserved.
3 // http://ceres-solver.org/
5 // Redistribution and use in source and binary forms, with or without
6 // modification, are permitted provided that the following conditions are met:
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.
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.
29 // Author: sameeragarwal@google.com (Sameer Agarwal)
31 // A simple C++ interface to the SuiteSparse and CHOLMOD libraries.
33 #ifndef CERES_INTERNAL_SUITESPARSE_H_
34 #define CERES_INTERNAL_SUITESPARSE_H_
36 // This include must come before any #ifndef check on Ceres compile options.
37 #include "ceres/internal/port.h"
39 #ifndef CERES_NO_SUITESPARSE
44 #include "SuiteSparseQR.hpp"
45 #include "ceres/linear_solver.h"
46 #include "ceres/sparse_cholesky.h"
48 #include "glog/logging.h"
50 // Before SuiteSparse version 4.2.0, cholmod_camd was only enabled
51 // if SuiteSparse was compiled with Metis support. This makes
52 // calling and linking into cholmod_camd problematic even though it
53 // has nothing to do with Metis. This has been fixed reliably in
56 // The fix was actually committed in 4.1.0, but there is
57 // some confusion about a silent update to the tar ball, so we are
58 // being conservative and choosing the next minor version where
60 #if (SUITESPARSE_VERSION < 4002)
64 // UF_long is deprecated but SuiteSparse_long is only available in
65 // newer versions of SuiteSparse. So for older versions of
66 // SuiteSparse, we define SuiteSparse_long to be the same as UF_long,
67 // which is what recent versions of SuiteSparse do anyways.
68 #ifndef SuiteSparse_long
69 #define SuiteSparse_long UF_long
75 class CompressedRowSparseMatrix;
76 class TripletSparseMatrix;
78 // The raw CHOLMOD and SuiteSparseQR libraries have a slightly
79 // cumbersome c like calling format. This object abstracts it away and
80 // provides the user with a simpler interface. The methods here cannot
81 // be static as a cholmod_common object serves as a global variable
82 // for all cholmod function calls.
88 // Functions for building cholmod_sparse objects from sparse
89 // matrices stored in triplet form. The matrix A is not
90 // modifed. Called owns the result.
91 cholmod_sparse* CreateSparseMatrix(TripletSparseMatrix* A);
93 // This function works like CreateSparseMatrix, except that the
94 // return value corresponds to A' rather than A.
95 cholmod_sparse* CreateSparseMatrixTranspose(TripletSparseMatrix* A);
97 // Create a cholmod_sparse wrapper around the contents of A. This is
98 // a shallow object, which refers to the contents of A and does not
99 // use the SuiteSparse machinery to allocate memory.
100 cholmod_sparse CreateSparseMatrixTransposeView(CompressedRowSparseMatrix* A);
102 // Given a vector x, build a cholmod_dense vector of size out_size
103 // with the first in_size entries copied from x. If x is NULL, then
104 // an all zeros vector is returned. Caller owns the result.
105 cholmod_dense* CreateDenseVector(const double* x, int in_size, int out_size);
107 // The matrix A is scaled using the matrix whose diagonal is the
108 // vector scale. mode describes how scaling is applied. Possible
109 // values are CHOLMOD_ROW for row scaling - diag(scale) * A,
110 // CHOLMOD_COL for column scaling - A * diag(scale) and CHOLMOD_SYM
111 // for symmetric scaling which scales both the rows and the columns
112 // - diag(scale) * A * diag(scale).
113 void Scale(cholmod_dense* scale, int mode, cholmod_sparse* A) {
114 cholmod_scale(scale, mode, A, &cc_);
117 // Create and return a matrix m = A * A'. Caller owns the
118 // result. The matrix A is not modified.
119 cholmod_sparse* AATranspose(cholmod_sparse* A) {
120 cholmod_sparse*m = cholmod_aat(A, NULL, A->nrow, 1, &cc_);
121 m->stype = 1; // Pay attention to the upper triangular part.
125 // y = alpha * A * x + beta * y. Only y is modified.
126 void SparseDenseMultiply(cholmod_sparse* A, double alpha, double beta,
127 cholmod_dense* x, cholmod_dense* y) {
128 double alpha_[2] = {alpha, 0};
129 double beta_[2] = {beta, 0};
130 cholmod_sdmult(A, 0, alpha_, beta_, x, y, &cc_);
133 // Find an ordering of A or AA' (if A is unsymmetric) that minimizes
134 // the fill-in in the Cholesky factorization of the corresponding
135 // matrix. This is done by using the AMD algorithm.
137 // Using this ordering, the symbolic Cholesky factorization of A (or
138 // AA') is computed and returned.
140 // A is not modified, only the pattern of non-zeros of A is used,
141 // the actual numerical values in A are of no consequence.
143 // message contains an explanation of the failures if any.
145 // Caller owns the result.
146 cholmod_factor* AnalyzeCholesky(cholmod_sparse* A, std::string* message);
148 cholmod_factor* BlockAnalyzeCholesky(cholmod_sparse* A,
149 const std::vector<int>& row_blocks,
150 const std::vector<int>& col_blocks,
151 std::string* message);
153 // If A is symmetric, then compute the symbolic Cholesky
154 // factorization of A(ordering, ordering). If A is unsymmetric, then
155 // compute the symbolic factorization of
156 // A(ordering,:) A(ordering,:)'.
158 // A is not modified, only the pattern of non-zeros of A is used,
159 // the actual numerical values in A are of no consequence.
161 // message contains an explanation of the failures if any.
163 // Caller owns the result.
164 cholmod_factor* AnalyzeCholeskyWithUserOrdering(
166 const std::vector<int>& ordering,
167 std::string* message);
169 // Perform a symbolic factorization of A without re-ordering A. No
170 // postordering of the elimination tree is performed. This ensures
171 // that the symbolic factor does not introduce an extra permutation
172 // on the matrix. See the documentation for CHOLMOD for more details.
174 // message contains an explanation of the failures if any.
175 cholmod_factor* AnalyzeCholeskyWithNaturalOrdering(cholmod_sparse* A,
176 std::string* message);
178 // Use the symbolic factorization in L, to find the numerical
179 // factorization for the matrix A or AA^T. Return true if
180 // successful, false otherwise. L contains the numeric factorization
183 // message contains an explanation of the failures if any.
184 LinearSolverTerminationType Cholesky(cholmod_sparse* A,
186 std::string* message);
188 // Given a Cholesky factorization of a matrix A = LL^T, solve the
189 // linear system Ax = b, and return the result. If the Solve fails
190 // NULL is returned. Caller owns the result.
192 // message contains an explanation of the failures if any.
193 cholmod_dense* Solve(cholmod_factor* L, cholmod_dense* b, std::string* message);
195 // By virtue of the modeling layer in Ceres being block oriented,
196 // all the matrices used by Ceres are also block oriented. When
197 // doing sparse direct factorization of these matrices the
198 // fill-reducing ordering algorithms (in particular AMD) can either
199 // be run on the block or the scalar form of these matrices. The two
200 // SuiteSparse::AnalyzeCholesky methods allows the the client to
201 // compute the symbolic factorization of a matrix by either using
202 // AMD on the matrix or a user provided ordering of the rows.
204 // But since the underlying matrices are block oriented, it is worth
205 // running AMD on just the block structre of these matrices and then
206 // lifting these block orderings to a full scalar ordering. This
207 // preserves the block structure of the permuted matrix, and exposes
208 // more of the super-nodal structure of the matrix to the numerical
209 // factorization routines.
211 // Find the block oriented AMD ordering of a matrix A, whose row and
212 // column blocks are given by row_blocks, and col_blocks
213 // respectively. The matrix may or may not be symmetric. The entries
214 // of col_blocks do not need to sum to the number of columns in
215 // A. If this is the case, only the first sum(col_blocks) are used
216 // to compute the ordering.
217 bool BlockAMDOrdering(const cholmod_sparse* A,
218 const std::vector<int>& row_blocks,
219 const std::vector<int>& col_blocks,
220 std::vector<int>* ordering);
222 // Find a fill reducing approximate minimum degree
223 // ordering. ordering is expected to be large enough to hold the
225 bool ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix, int* ordering);
228 // Before SuiteSparse version 4.2.0, cholmod_camd was only enabled
229 // if SuiteSparse was compiled with Metis support. This makes
230 // calling and linking into cholmod_camd problematic even though it
231 // has nothing to do with Metis. This has been fixed reliably in
234 // The fix was actually committed in 4.1.0, but there is
235 // some confusion about a silent update to the tar ball, so we are
236 // being conservative and choosing the next minor version where
237 // things are stable.
238 static bool IsConstrainedApproximateMinimumDegreeOrderingAvailable() {
239 return (SUITESPARSE_VERSION > 4001);
242 // Find a fill reducing approximate minimum degree
243 // ordering. constraints is an array which associates with each
244 // column of the matrix an elimination group. i.e., all columns in
245 // group 0 are eliminated first, all columns in group 1 are
246 // eliminated next etc. This function finds a fill reducing ordering
247 // that obeys these constraints.
249 // Calling ApproximateMinimumDegreeOrdering is equivalent to calling
250 // ConstrainedApproximateMinimumDegreeOrdering with a constraint
251 // array that puts all columns in the same elimination group.
253 // If CERES_NO_CAMD is defined then calling this function will
254 // result in a crash.
255 bool ConstrainedApproximateMinimumDegreeOrdering(cholmod_sparse* matrix,
259 void Free(cholmod_sparse* m) { cholmod_free_sparse(&m, &cc_); }
260 void Free(cholmod_dense* m) { cholmod_free_dense(&m, &cc_); }
261 void Free(cholmod_factor* m) { cholmod_free_factor(&m, &cc_); }
263 void Print(cholmod_sparse* m, const std::string& name) {
264 cholmod_print_sparse(m, const_cast<char*>(name.c_str()), &cc_);
267 void Print(cholmod_dense* m, const std::string& name) {
268 cholmod_print_dense(m, const_cast<char*>(name.c_str()), &cc_);
271 void Print(cholmod_triplet* m, const std::string& name) {
272 cholmod_print_triplet(m, const_cast<char*>(name.c_str()), &cc_);
275 cholmod_common* mutable_cc() { return &cc_; }
281 class SuiteSparseCholesky : public SparseCholesky {
283 static SuiteSparseCholesky* Create(const OrderingType ordering_type);
285 // SparseCholesky interface.
286 virtual ~SuiteSparseCholesky();
287 virtual CompressedRowSparseMatrix::StorageType StorageType() const;
288 virtual LinearSolverTerminationType Factorize(
289 CompressedRowSparseMatrix* lhs, std::string* message);
290 virtual LinearSolverTerminationType Solve(const double* rhs,
292 std::string* message);
294 SuiteSparseCholesky(const OrderingType ordering_type);
296 const OrderingType ordering_type_;
298 cholmod_factor* factor_;
301 } // namespace internal
304 #else // CERES_NO_SUITESPARSE
306 typedef void cholmod_factor;
313 // Defining this static function even when SuiteSparse is not
314 // available, allows client code to check for the presence of CAMD
315 // without checking for the absence of the CERES_NO_CAMD symbol.
317 // This is safer because the symbol maybe missing due to a user
318 // accidently not including suitesparse.h in their code when
319 // checking for the symbol.
320 static bool IsConstrainedApproximateMinimumDegreeOrderingAvailable() {
324 void Free(void* arg) {}
327 } // namespace internal
330 #endif // CERES_NO_SUITESPARSE
332 #endif // CERES_INTERNAL_SUITESPARSE_H_