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@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
44namespace ceres {
45namespace internal {
46
47SuiteSparse::SuiteSparse() {
48  cholmod_start(&cc_);
49}
50
51SuiteSparse::~SuiteSparse() {
52  cholmod_finish(&cc_);
53}
54
55cholmod_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
74cholmod_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
95cholmod_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
116cholmod_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
127cholmod_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
151cholmod_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
163cholmod_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
186cholmod_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
206bool 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
248LinearSolverTerminationType 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@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@googlegroups.com.",
305                       cc_.status);
306      return LINEAR_SOLVER_FATAL_ERROR;
307  }
308
309  return LINEAR_SOLVER_FATAL_ERROR;
310}
311
312cholmod_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
323bool SuiteSparse::ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix,
324                                                   int* ordering) {
325  return cholmod_amd(matrix, NULL, 0, ordering, &cc_);
326}
327
328bool 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