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/
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
#ifndef CERES_NO_SUITESPARSE
33
#include "ceres/suitesparse.h"
36
#include "ceres/compressed_row_sparse_matrix.h"
37
#include "ceres/triplet_sparse_matrix.h"
41
cholmod_sparse* SuiteSparse::CreateSparseMatrix(TripletSparseMatrix* A) {
42
cholmod_triplet triplet;
44
triplet.nrow = A->num_rows();
45
triplet.ncol = A->num_cols();
46
triplet.nzmax = A->max_num_nonzeros();
47
triplet.nnz = A->num_nonzeros();
48
triplet.i = reinterpret_cast<void*>(A->mutable_rows());
49
triplet.j = reinterpret_cast<void*>(A->mutable_cols());
50
triplet.x = reinterpret_cast<void*>(A->mutable_values());
51
triplet.stype = 0; // Matrix is not symmetric.
52
triplet.itype = CHOLMOD_INT;
53
triplet.xtype = CHOLMOD_REAL;
54
triplet.dtype = CHOLMOD_DOUBLE;
56
return cholmod_triplet_to_sparse(&triplet, triplet.nnz, &cc_);
60
cholmod_sparse* SuiteSparse::CreateSparseMatrixTranspose(
61
TripletSparseMatrix* A) {
62
cholmod_triplet triplet;
64
triplet.ncol = A->num_rows(); // swap row and columns
65
triplet.nrow = A->num_cols();
66
triplet.nzmax = A->max_num_nonzeros();
67
triplet.nnz = A->num_nonzeros();
69
// swap rows and columns
70
triplet.j = reinterpret_cast<void*>(A->mutable_rows());
71
triplet.i = reinterpret_cast<void*>(A->mutable_cols());
72
triplet.x = reinterpret_cast<void*>(A->mutable_values());
73
triplet.stype = 0; // Matrix is not symmetric.
74
triplet.itype = CHOLMOD_INT;
75
triplet.xtype = CHOLMOD_REAL;
76
triplet.dtype = CHOLMOD_DOUBLE;
78
return cholmod_triplet_to_sparse(&triplet, triplet.nnz, &cc_);
81
cholmod_sparse* SuiteSparse::CreateSparseMatrixTransposeView(
82
CompressedRowSparseMatrix* A) {
83
cholmod_sparse* m = new cholmod_sparse_struct;
84
m->nrow = A->num_cols();
85
m->ncol = A->num_rows();
86
m->nzmax = A->num_nonzeros();
88
m->p = reinterpret_cast<void*>(A->mutable_rows());
89
m->i = reinterpret_cast<void*>(A->mutable_cols());
90
m->x = reinterpret_cast<void*>(A->mutable_values());
92
m->stype = 0; // Matrix is not symmetric.
93
m->itype = CHOLMOD_INT;
94
m->xtype = CHOLMOD_REAL;
95
m->dtype = CHOLMOD_DOUBLE;
102
cholmod_dense* SuiteSparse::CreateDenseVector(const double* x,
105
CHECK_LE(in_size, out_size);
106
cholmod_dense* v = cholmod_zeros(out_size, 1, CHOLMOD_REAL, &cc_);
108
memcpy(v->x, x, in_size*sizeof(*x));
113
cholmod_factor* SuiteSparse::AnalyzeCholesky(cholmod_sparse* A) {
114
cholmod_factor* factor = cholmod_analyze(A, &cc_);
115
CHECK_EQ(cc_.status, CHOLMOD_OK)
116
<< "Cholmod symbolic analysis failed " << cc_.status;
117
CHECK_NOTNULL(factor);
121
bool SuiteSparse::Cholesky(cholmod_sparse* A, cholmod_factor* L) {
125
cc_.quick_return_if_not_posdef = 1;
126
int status = cholmod_factorize(A, L, &cc_);
127
switch (cc_.status) {
128
case CHOLMOD_NOT_INSTALLED:
129
LOG(WARNING) << "Cholmod failure: method not installed.";
131
case CHOLMOD_OUT_OF_MEMORY:
132
LOG(WARNING) << "Cholmod failure: out of memory.";
134
case CHOLMOD_TOO_LARGE:
135
LOG(WARNING) << "Cholmod failure: integer overflow occured.";
137
case CHOLMOD_INVALID:
138
LOG(WARNING) << "Cholmod failure: invalid input.";
140
case CHOLMOD_NOT_POSDEF:
141
// TODO(sameeragarwal): These two warnings require more
142
// sophisticated handling going forward. For now we will be
143
// strict and treat them as failures.
144
LOG(WARNING) << "Cholmod warning: matrix not positive definite.";
147
LOG(WARNING) << "Cholmod warning: D for LDL' or diag(L) or "
148
<< "LL' has tiny absolute value.";
154
LOG(WARNING) << "Cholmod failure: cholmod_factorize returned zero "
155
<< "but cholmod_common::status is CHOLMOD_OK."
156
<< "Please report this to ceres-solver@googlegroups.com.";
159
LOG(WARNING) << "Unknown cholmod return code. "
160
<< "Please report this to ceres-solver@googlegroups.com.";
166
cholmod_dense* SuiteSparse::Solve(cholmod_factor* L,
168
if (cc_.status != CHOLMOD_OK) {
169
LOG(WARNING) << "CHOLMOD status NOT OK";
173
return cholmod_solve(CHOLMOD_A, L, b, &cc_);
176
cholmod_dense* SuiteSparse::SolveCholesky(cholmod_sparse* A,
183
if (Cholesky(A, L)) {
190
} // namespace internal
193
#endif // CERES_NO_SUITESPARSE