suitesparse.h revision 1d2624a10e2c559f8ba9ef89eaa30832c0a83a96
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// A simple C++ interface to the SuiteSparse and CHOLMOD libraries. 32 33#ifndef CERES_INTERNAL_SUITESPARSE_H_ 34#define CERES_INTERNAL_SUITESPARSE_H_ 35 36 37#ifndef CERES_NO_SUITESPARSE 38 39#include <cstring> 40#include <string> 41#include <vector> 42 43#include "ceres/internal/port.h" 44#include "cholmod.h" 45#include "glog/logging.h" 46 47// Before SuiteSparse version 4.2.0, cholmod_camd was only enabled 48// if SuiteSparse was compiled with Metis support. This makes 49// calling and linking into cholmod_camd problematic even though it 50// has nothing to do with Metis. This has been fixed reliably in 51// 4.2.0. 52// 53// The fix was actually committed in 4.1.0, but there is 54// some confusion about a silent update to the tar ball, so we are 55// being conservative and choosing the next minor version where 56// things are stable. 57#if (SUITESPARSE_VERSION < 4002) 58#define CERES_NO_CAMD 59#endif 60 61namespace ceres { 62namespace internal { 63 64class CompressedRowSparseMatrix; 65class TripletSparseMatrix; 66 67// The raw CHOLMOD and SuiteSparseQR libraries have a slightly 68// cumbersome c like calling format. This object abstracts it away and 69// provides the user with a simpler interface. The methods here cannot 70// be static as a cholmod_common object serves as a global variable 71// for all cholmod function calls. 72class SuiteSparse { 73 public: 74 SuiteSparse(); 75 ~SuiteSparse(); 76 77 // Functions for building cholmod_sparse objects from sparse 78 // matrices stored in triplet form. The matrix A is not 79 // modifed. Called owns the result. 80 cholmod_sparse* CreateSparseMatrix(TripletSparseMatrix* A); 81 82 // This function works like CreateSparseMatrix, except that the 83 // return value corresponds to A' rather than A. 84 cholmod_sparse* CreateSparseMatrixTranspose(TripletSparseMatrix* A); 85 86 // Create a cholmod_sparse wrapper around the contents of A. This is 87 // a shallow object, which refers to the contents of A and does not 88 // use the SuiteSparse machinery to allocate memory. 89 cholmod_sparse CreateSparseMatrixTransposeView(CompressedRowSparseMatrix* A); 90 91 // Given a vector x, build a cholmod_dense vector of size out_size 92 // with the first in_size entries copied from x. If x is NULL, then 93 // an all zeros vector is returned. Caller owns the result. 94 cholmod_dense* CreateDenseVector(const double* x, int in_size, int out_size); 95 96 // The matrix A is scaled using the matrix whose diagonal is the 97 // vector scale. mode describes how scaling is applied. Possible 98 // values are CHOLMOD_ROW for row scaling - diag(scale) * A, 99 // CHOLMOD_COL for column scaling - A * diag(scale) and CHOLMOD_SYM 100 // for symmetric scaling which scales both the rows and the columns 101 // - diag(scale) * A * diag(scale). 102 void Scale(cholmod_dense* scale, int mode, cholmod_sparse* A) { 103 cholmod_scale(scale, mode, A, &cc_); 104 } 105 106 // Create and return a matrix m = A * A'. Caller owns the 107 // result. The matrix A is not modified. 108 cholmod_sparse* AATranspose(cholmod_sparse* A) { 109 cholmod_sparse*m = cholmod_aat(A, NULL, A->nrow, 1, &cc_); 110 m->stype = 1; // Pay attention to the upper triangular part. 111 return m; 112 } 113 114 // y = alpha * A * x + beta * y. Only y is modified. 115 void SparseDenseMultiply(cholmod_sparse* A, double alpha, double beta, 116 cholmod_dense* x, cholmod_dense* y) { 117 double alpha_[2] = {alpha, 0}; 118 double beta_[2] = {beta, 0}; 119 cholmod_sdmult(A, 0, alpha_, beta_, x, y, &cc_); 120 } 121 122 // Find an ordering of A or AA' (if A is unsymmetric) that minimizes 123 // the fill-in in the Cholesky factorization of the corresponding 124 // matrix. This is done by using the AMD algorithm. 125 // 126 // Using this ordering, the symbolic Cholesky factorization of A (or 127 // AA') is computed and returned. 128 // 129 // A is not modified, only the pattern of non-zeros of A is used, 130 // the actual numerical values in A are of no consequence. 131 // 132 // Caller owns the result. 133 cholmod_factor* AnalyzeCholesky(cholmod_sparse* A); 134 135 cholmod_factor* BlockAnalyzeCholesky(cholmod_sparse* A, 136 const vector<int>& row_blocks, 137 const vector<int>& col_blocks); 138 139 // If A is symmetric, then compute the symbolic Cholesky 140 // factorization of A(ordering, ordering). If A is unsymmetric, then 141 // compute the symbolic factorization of 142 // A(ordering,:) A(ordering,:)'. 143 // 144 // A is not modified, only the pattern of non-zeros of A is used, 145 // the actual numerical values in A are of no consequence. 146 // 147 // Caller owns the result. 148 cholmod_factor* AnalyzeCholeskyWithUserOrdering(cholmod_sparse* A, 149 const vector<int>& ordering); 150 151 // Perform a symbolic factorization of A without re-ordering A. No 152 // postordering of the elimination tree is performed. This ensures 153 // that the symbolic factor does not introduce an extra permutation 154 // on the matrix. See the documentation for CHOLMOD for more details. 155 cholmod_factor* AnalyzeCholeskyWithNaturalOrdering(cholmod_sparse* A); 156 157 // Use the symbolic factorization in L, to find the numerical 158 // factorization for the matrix A or AA^T. Return true if 159 // successful, false otherwise. L contains the numeric factorization 160 // on return. 161 bool Cholesky(cholmod_sparse* A, cholmod_factor* L); 162 163 // Given a Cholesky factorization of a matrix A = LL^T, solve the 164 // linear system Ax = b, and return the result. If the Solve fails 165 // NULL is returned. Caller owns the result. 166 cholmod_dense* Solve(cholmod_factor* L, cholmod_dense* b); 167 168 // Combine the calls to Cholesky and Solve into a single call. If 169 // the cholesky factorization or the solve fails, return 170 // NULL. Caller owns the result. 171 cholmod_dense* SolveCholesky(cholmod_sparse* A, 172 cholmod_factor* L, 173 cholmod_dense* b); 174 175 // By virtue of the modeling layer in Ceres being block oriented, 176 // all the matrices used by Ceres are also block oriented. When 177 // doing sparse direct factorization of these matrices the 178 // fill-reducing ordering algorithms (in particular AMD) can either 179 // be run on the block or the scalar form of these matrices. The two 180 // SuiteSparse::AnalyzeCholesky methods allows the the client to 181 // compute the symbolic factorization of a matrix by either using 182 // AMD on the matrix or a user provided ordering of the rows. 183 // 184 // But since the underlying matrices are block oriented, it is worth 185 // running AMD on just the block structre of these matrices and then 186 // lifting these block orderings to a full scalar ordering. This 187 // preserves the block structure of the permuted matrix, and exposes 188 // more of the super-nodal structure of the matrix to the numerical 189 // factorization routines. 190 // 191 // Find the block oriented AMD ordering of a matrix A, whose row and 192 // column blocks are given by row_blocks, and col_blocks 193 // respectively. The matrix may or may not be symmetric. The entries 194 // of col_blocks do not need to sum to the number of columns in 195 // A. If this is the case, only the first sum(col_blocks) are used 196 // to compute the ordering. 197 bool BlockAMDOrdering(const cholmod_sparse* A, 198 const vector<int>& row_blocks, 199 const vector<int>& col_blocks, 200 vector<int>* ordering); 201 202 // Find a fill reducing approximate minimum degree 203 // ordering. ordering is expected to be large enough to hold the 204 // ordering. 205 void ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix, int* ordering); 206 207 208 // Before SuiteSparse version 4.2.0, cholmod_camd was only enabled 209 // if SuiteSparse was compiled with Metis support. This makes 210 // calling and linking into cholmod_camd problematic even though it 211 // has nothing to do with Metis. This has been fixed reliably in 212 // 4.2.0. 213 // 214 // The fix was actually committed in 4.1.0, but there is 215 // some confusion about a silent update to the tar ball, so we are 216 // being conservative and choosing the next minor version where 217 // things are stable. 218 static bool IsConstrainedApproximateMinimumDegreeOrderingAvailable() { 219 return (SUITESPARSE_VERSION>4001); 220 } 221 222 // Find a fill reducing approximate minimum degree 223 // ordering. constraints is an array which associates with each 224 // column of the matrix an elimination group. i.e., all columns in 225 // group 0 are eliminated first, all columns in group 1 are 226 // eliminated next etc. This function finds a fill reducing ordering 227 // that obeys these constraints. 228 // 229 // Calling ApproximateMinimumDegreeOrdering is equivalent to calling 230 // ConstrainedApproximateMinimumDegreeOrdering with a constraint 231 // array that puts all columns in the same elimination group. 232 // 233 // If CERES_NO_CAMD is defined then calling this function will 234 // result in a crash. 235 void ConstrainedApproximateMinimumDegreeOrdering(cholmod_sparse* matrix, 236 int* constraints, 237 int* ordering); 238 239 void Free(cholmod_sparse* m) { cholmod_free_sparse(&m, &cc_); } 240 void Free(cholmod_dense* m) { cholmod_free_dense(&m, &cc_); } 241 void Free(cholmod_factor* m) { cholmod_free_factor(&m, &cc_); } 242 243 void Print(cholmod_sparse* m, const string& name) { 244 cholmod_print_sparse(m, const_cast<char*>(name.c_str()), &cc_); 245 } 246 247 void Print(cholmod_dense* m, const string& name) { 248 cholmod_print_dense(m, const_cast<char*>(name.c_str()), &cc_); 249 } 250 251 void Print(cholmod_triplet* m, const string& name) { 252 cholmod_print_triplet(m, const_cast<char*>(name.c_str()), &cc_); 253 } 254 255 cholmod_common* mutable_cc() { return &cc_; } 256 257 private: 258 cholmod_common cc_; 259}; 260 261} // namespace internal 262} // namespace ceres 263 264#endif // CERES_NO_SUITESPARSE 265 266#endif // CERES_INTERNAL_SUITESPARSE_H_ 267