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