1 // Ceres Solver - A fast non-linear least squares minimizer 2 // Copyright 2010, 2011, 2012 Google Inc. All rights reserved. 3 // http://code.google.com/p/ceres-solver/ 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 (at) google.com (Sameer Agarwal) 30 31 // This include must come before any #ifndef check on Ceres compile options. 32 #include "ceres/internal/port.h" 33 34 #ifndef CERES_NO_SUITESPARSE 35 #include "ceres/suitesparse.h" 36 37 #include <vector> 38 #include "cholmod.h" 39 #include "ceres/compressed_col_sparse_matrix_utils.h" 40 #include "ceres/compressed_row_sparse_matrix.h" 41 #include "ceres/linear_solver.h" 42 #include "ceres/triplet_sparse_matrix.h" 43 44 namespace ceres { 45 namespace internal { 46 47 SuiteSparse::SuiteSparse() { 48 cholmod_start(&cc_); 49 } 50 51 SuiteSparse::~SuiteSparse() { 52 cholmod_finish(&cc_); 53 } 54 55 cholmod_sparse* SuiteSparse::CreateSparseMatrix(TripletSparseMatrix* A) { 56 cholmod_triplet triplet; 57 58 triplet.nrow = A->num_rows(); 59 triplet.ncol = A->num_cols(); 60 triplet.nzmax = A->max_num_nonzeros(); 61 triplet.nnz = A->num_nonzeros(); 62 triplet.i = reinterpret_cast<void*>(A->mutable_rows()); 63 triplet.j = reinterpret_cast<void*>(A->mutable_cols()); 64 triplet.x = reinterpret_cast<void*>(A->mutable_values()); 65 triplet.stype = 0; // Matrix is not symmetric. 66 triplet.itype = CHOLMOD_INT; 67 triplet.xtype = CHOLMOD_REAL; 68 triplet.dtype = CHOLMOD_DOUBLE; 69 70 return cholmod_triplet_to_sparse(&triplet, triplet.nnz, &cc_); 71 } 72 73 74 cholmod_sparse* SuiteSparse::CreateSparseMatrixTranspose( 75 TripletSparseMatrix* A) { 76 cholmod_triplet triplet; 77 78 triplet.ncol = A->num_rows(); // swap row and columns 79 triplet.nrow = A->num_cols(); 80 triplet.nzmax = A->max_num_nonzeros(); 81 triplet.nnz = A->num_nonzeros(); 82 83 // swap rows and columns 84 triplet.j = reinterpret_cast<void*>(A->mutable_rows()); 85 triplet.i = reinterpret_cast<void*>(A->mutable_cols()); 86 triplet.x = reinterpret_cast<void*>(A->mutable_values()); 87 triplet.stype = 0; // Matrix is not symmetric. 88 triplet.itype = CHOLMOD_INT; 89 triplet.xtype = CHOLMOD_REAL; 90 triplet.dtype = CHOLMOD_DOUBLE; 91 92 return cholmod_triplet_to_sparse(&triplet, triplet.nnz, &cc_); 93 } 94 95 cholmod_sparse SuiteSparse::CreateSparseMatrixTransposeView( 96 CompressedRowSparseMatrix* A) { 97 cholmod_sparse m; 98 m.nrow = A->num_cols(); 99 m.ncol = A->num_rows(); 100 m.nzmax = A->num_nonzeros(); 101 m.nz = NULL; 102 m.p = reinterpret_cast<void*>(A->mutable_rows()); 103 m.i = reinterpret_cast<void*>(A->mutable_cols()); 104 m.x = reinterpret_cast<void*>(A->mutable_values()); 105 m.z = NULL; 106 m.stype = 0; // Matrix is not symmetric. 107 m.itype = CHOLMOD_INT; 108 m.xtype = CHOLMOD_REAL; 109 m.dtype = CHOLMOD_DOUBLE; 110 m.sorted = 1; 111 m.packed = 1; 112 113 return m; 114 } 115 116 cholmod_dense* SuiteSparse::CreateDenseVector(const double* x, 117 int in_size, 118 int out_size) { 119 CHECK_LE(in_size, out_size); 120 cholmod_dense* v = cholmod_zeros(out_size, 1, CHOLMOD_REAL, &cc_); 121 if (x != NULL) { 122 memcpy(v->x, x, in_size*sizeof(*x)); 123 } 124 return v; 125 } 126 127 cholmod_factor* SuiteSparse::AnalyzeCholesky(cholmod_sparse* A, 128 string* message) { 129 // Cholmod can try multiple re-ordering strategies to find a fill 130 // reducing ordering. Here we just tell it use AMD with automatic 131 // matrix dependence choice of supernodal versus simplicial 132 // factorization. 133 cc_.nmethods = 1; 134 cc_.method[0].ordering = CHOLMOD_AMD; 135 cc_.supernodal = CHOLMOD_AUTO; 136 137 cholmod_factor* factor = cholmod_analyze(A, &cc_); 138 if (VLOG_IS_ON(2)) { 139 cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_); 140 } 141 142 if (cc_.status != CHOLMOD_OK) { 143 *message = StringPrintf("cholmod_analyze failed. error code: %d", 144 cc_.status); 145 return NULL; 146 } 147 148 return CHECK_NOTNULL(factor); 149 } 150 151 cholmod_factor* SuiteSparse::BlockAnalyzeCholesky( 152 cholmod_sparse* A, 153 const vector<int>& row_blocks, 154 const vector<int>& col_blocks, 155 string* message) { 156 vector<int> ordering; 157 if (!BlockAMDOrdering(A, row_blocks, col_blocks, &ordering)) { 158 return NULL; 159 } 160 return AnalyzeCholeskyWithUserOrdering(A, ordering, message); 161 } 162 163 cholmod_factor* SuiteSparse::AnalyzeCholeskyWithUserOrdering( 164 cholmod_sparse* A, 165 const vector<int>& ordering, 166 string* message) { 167 CHECK_EQ(ordering.size(), A->nrow); 168 169 cc_.nmethods = 1; 170 cc_.method[0].ordering = CHOLMOD_GIVEN; 171 172 cholmod_factor* factor = 173 cholmod_analyze_p(A, const_cast<int*>(&ordering[0]), NULL, 0, &cc_); 174 if (VLOG_IS_ON(2)) { 175 cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_); 176 } 177 if (cc_.status != CHOLMOD_OK) { 178 *message = StringPrintf("cholmod_analyze failed. error code: %d", 179 cc_.status); 180 return NULL; 181 } 182 183 return CHECK_NOTNULL(factor); 184 } 185 186 cholmod_factor* SuiteSparse::AnalyzeCholeskyWithNaturalOrdering( 187 cholmod_sparse* A, 188 string* message) { 189 cc_.nmethods = 1; 190 cc_.method[0].ordering = CHOLMOD_NATURAL; 191 cc_.postorder = 0; 192 193 cholmod_factor* factor = cholmod_analyze(A, &cc_); 194 if (VLOG_IS_ON(2)) { 195 cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_); 196 } 197 if (cc_.status != CHOLMOD_OK) { 198 *message = StringPrintf("cholmod_analyze failed. error code: %d", 199 cc_.status); 200 return NULL; 201 } 202 203 return CHECK_NOTNULL(factor); 204 } 205 206 bool SuiteSparse::BlockAMDOrdering(const cholmod_sparse* A, 207 const vector<int>& row_blocks, 208 const vector<int>& col_blocks, 209 vector<int>* ordering) { 210 const int num_row_blocks = row_blocks.size(); 211 const int num_col_blocks = col_blocks.size(); 212 213 // Arrays storing the compressed column structure of the matrix 214 // incoding the block sparsity of A. 215 vector<int> block_cols; 216 vector<int> block_rows; 217 218 CompressedColumnScalarMatrixToBlockMatrix(reinterpret_cast<const int*>(A->i), 219 reinterpret_cast<const int*>(A->p), 220 row_blocks, 221 col_blocks, 222 &block_rows, 223 &block_cols); 224 225 cholmod_sparse_struct block_matrix; 226 block_matrix.nrow = num_row_blocks; 227 block_matrix.ncol = num_col_blocks; 228 block_matrix.nzmax = block_rows.size(); 229 block_matrix.p = reinterpret_cast<void*>(&block_cols[0]); 230 block_matrix.i = reinterpret_cast<void*>(&block_rows[0]); 231 block_matrix.x = NULL; 232 block_matrix.stype = A->stype; 233 block_matrix.itype = CHOLMOD_INT; 234 block_matrix.xtype = CHOLMOD_PATTERN; 235 block_matrix.dtype = CHOLMOD_DOUBLE; 236 block_matrix.sorted = 1; 237 block_matrix.packed = 1; 238 239 vector<int> block_ordering(num_row_blocks); 240 if (!cholmod_amd(&block_matrix, NULL, 0, &block_ordering[0], &cc_)) { 241 return false; 242 } 243 244 BlockOrderingToScalarOrdering(row_blocks, block_ordering, ordering); 245 return true; 246 } 247 248 LinearSolverTerminationType SuiteSparse::Cholesky(cholmod_sparse* A, 249 cholmod_factor* L, 250 string* message) { 251 CHECK_NOTNULL(A); 252 CHECK_NOTNULL(L); 253 254 // Save the current print level and silence CHOLMOD, otherwise 255 // CHOLMOD is prone to dumping stuff to stderr, which can be 256 // distracting when the error (matrix is indefinite) is not a fatal 257 // failure. 258 const int old_print_level = cc_.print; 259 cc_.print = 0; 260 261 cc_.quick_return_if_not_posdef = 1; 262 int cholmod_status = cholmod_factorize(A, L, &cc_); 263 cc_.print = old_print_level; 264 265 // TODO(sameeragarwal): This switch statement is not consistent. It 266 // treats all kinds of CHOLMOD failures as warnings. Some of these 267 // like out of memory are definitely not warnings. The problem is 268 // that the return value Cholesky is two valued, but the state of 269 // the linear solver is really three valued. SUCCESS, 270 // NON_FATAL_FAILURE (e.g., indefinite matrix) and FATAL_FAILURE 271 // (e.g. out of memory). 272 switch (cc_.status) { 273 case CHOLMOD_NOT_INSTALLED: 274 *message = "CHOLMOD failure: Method not installed."; 275 return LINEAR_SOLVER_FATAL_ERROR; 276 case CHOLMOD_OUT_OF_MEMORY: 277 *message = "CHOLMOD failure: Out of memory."; 278 return LINEAR_SOLVER_FATAL_ERROR; 279 case CHOLMOD_TOO_LARGE: 280 *message = "CHOLMOD failure: Integer overflow occured."; 281 return LINEAR_SOLVER_FATAL_ERROR; 282 case CHOLMOD_INVALID: 283 *message = "CHOLMOD failure: Invalid input."; 284 return LINEAR_SOLVER_FATAL_ERROR; 285 case CHOLMOD_NOT_POSDEF: 286 *message = "CHOLMOD warning: Matrix not positive definite."; 287 return LINEAR_SOLVER_FAILURE; 288 case CHOLMOD_DSMALL: 289 *message = "CHOLMOD warning: D for LDL' or diag(L) or " 290 "LL' has tiny absolute value."; 291 return LINEAR_SOLVER_FAILURE; 292 case CHOLMOD_OK: 293 if (cholmod_status != 0) { 294 return LINEAR_SOLVER_SUCCESS; 295 } 296 297 *message = "CHOLMOD failure: cholmod_factorize returned false " 298 "but cholmod_common::status is CHOLMOD_OK." 299 "Please report this to ceres-solver (at) googlegroups.com."; 300 return LINEAR_SOLVER_FATAL_ERROR; 301 default: 302 *message = 303 StringPrintf("Unknown cholmod return code: %d. " 304 "Please report this to ceres-solver (at) googlegroups.com.", 305 cc_.status); 306 return LINEAR_SOLVER_FATAL_ERROR; 307 } 308 309 return LINEAR_SOLVER_FATAL_ERROR; 310 } 311 312 cholmod_dense* SuiteSparse::Solve(cholmod_factor* L, 313 cholmod_dense* b, 314 string* message) { 315 if (cc_.status != CHOLMOD_OK) { 316 *message = "cholmod_solve failed. CHOLMOD status is not CHOLMOD_OK"; 317 return NULL; 318 } 319 320 return cholmod_solve(CHOLMOD_A, L, b, &cc_); 321 } 322 323 bool SuiteSparse::ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix, 324 int* ordering) { 325 return cholmod_amd(matrix, NULL, 0, ordering, &cc_); 326 } 327 328 bool SuiteSparse::ConstrainedApproximateMinimumDegreeOrdering( 329 cholmod_sparse* matrix, 330 int* constraints, 331 int* ordering) { 332 #ifndef CERES_NO_CAMD 333 return cholmod_camd(matrix, NULL, 0, constraints, ordering, &cc_); 334 #else 335 LOG(FATAL) << "Congratulations you have found a bug in Ceres." 336 << "Ceres Solver was compiled with SuiteSparse " 337 << "version 4.1.0 or less. Calling this function " 338 << "in that case is a bug. Please contact the" 339 << "the Ceres Solver developers."; 340 return false; 341 #endif 342 } 343 344 } // namespace internal 345 } // namespace ceres 346 347 #endif // CERES_NO_SUITESPARSE 348
