10ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Ceres Solver - A fast non-linear least squares minimizer 20ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Copyright 2010, 2011, 2012 Google Inc. All rights reserved. 30ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// http://code.google.com/p/ceres-solver/ 40ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// 50ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Redistribution and use in source and binary forms, with or without 60ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// modification, are permitted provided that the following conditions are met: 70ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// 80ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// * Redistributions of source code must retain the above copyright notice, 90ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// this list of conditions and the following disclaimer. 100ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// * Redistributions in binary form must reproduce the above copyright notice, 110ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// this list of conditions and the following disclaimer in the documentation 120ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// and/or other materials provided with the distribution. 130ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// * Neither the name of Google Inc. nor the names of its contributors may be 140ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// used to endorse or promote products derived from this software without 150ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// specific prior written permission. 160ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// 170ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 180ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 190ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 200ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 210ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 220ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 230ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 240ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 250ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 260ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 270ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// POSSIBILITY OF SUCH DAMAGE. 280ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// 290ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Author: sameeragarwal@google.com (Sameer Agarwal) 300ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// 310ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// A simple C++ interface to the SuiteSparse and CHOLMOD libraries. 320ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 330ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#ifndef CERES_INTERNAL_SUITESPARSE_H_ 340ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#define CERES_INTERNAL_SUITESPARSE_H_ 350ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 3679397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez// This include must come before any #ifndef check on Ceres compile options. 3779397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez#include "ceres/internal/port.h" 381d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 390ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#ifndef CERES_NO_SUITESPARSE 400ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 410ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include <cstring> 420ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include <string> 430ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include <vector> 440ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 450ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include "ceres/internal/port.h" 4679397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez#include "ceres/linear_solver.h" 471d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include "cholmod.h" 481d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include "glog/logging.h" 49399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger#include "SuiteSparseQR.hpp" 501d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 511d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// Before SuiteSparse version 4.2.0, cholmod_camd was only enabled 521d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// if SuiteSparse was compiled with Metis support. This makes 531d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// calling and linking into cholmod_camd problematic even though it 541d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// has nothing to do with Metis. This has been fixed reliably in 551d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 4.2.0. 561d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// 571d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// The fix was actually committed in 4.1.0, but there is 581d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// some confusion about a silent update to the tar ball, so we are 591d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// being conservative and choosing the next minor version where 601d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling// things are stable. 611d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#if (SUITESPARSE_VERSION < 4002) 621d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#define CERES_NO_CAMD 631d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#endif 640ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 65399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger// UF_long is deprecated but SuiteSparse_long is only available in 66399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger// newer versions of SuiteSparse. So for older versions of 67399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger// SuiteSparse, we define SuiteSparse_long to be the same as UF_long, 68399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger// which is what recent versions of SuiteSparse do anyways. 69399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger#ifndef SuiteSparse_long 70399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger#define SuiteSparse_long UF_long 71399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger#endif 72399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger 730ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongnamespace ceres { 740ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongnamespace internal { 750ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 760ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongclass CompressedRowSparseMatrix; 770ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongclass TripletSparseMatrix; 780ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 790ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// The raw CHOLMOD and SuiteSparseQR libraries have a slightly 800ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// cumbersome c like calling format. This object abstracts it away and 810ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// provides the user with a simpler interface. The methods here cannot 820ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// be static as a cholmod_common object serves as a global variable 830ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// for all cholmod function calls. 840ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongclass SuiteSparse { 850ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong public: 861d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling SuiteSparse(); 871d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling ~SuiteSparse(); 880ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 890ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Functions for building cholmod_sparse objects from sparse 900ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // matrices stored in triplet form. The matrix A is not 910ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // modifed. Called owns the result. 920ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong cholmod_sparse* CreateSparseMatrix(TripletSparseMatrix* A); 930ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 940ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // This function works like CreateSparseMatrix, except that the 950ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // return value corresponds to A' rather than A. 960ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong cholmod_sparse* CreateSparseMatrixTranspose(TripletSparseMatrix* A); 970ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 980ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Create a cholmod_sparse wrapper around the contents of A. This is 990ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // a shallow object, which refers to the contents of A and does not 1001d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // use the SuiteSparse machinery to allocate memory. 1011d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling cholmod_sparse CreateSparseMatrixTransposeView(CompressedRowSparseMatrix* A); 1020ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1030ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Given a vector x, build a cholmod_dense vector of size out_size 1040ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // with the first in_size entries copied from x. If x is NULL, then 1050ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // an all zeros vector is returned. Caller owns the result. 1060ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong cholmod_dense* CreateDenseVector(const double* x, int in_size, int out_size); 1070ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1080ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // The matrix A is scaled using the matrix whose diagonal is the 1090ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // vector scale. mode describes how scaling is applied. Possible 1100ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // values are CHOLMOD_ROW for row scaling - diag(scale) * A, 1110ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // CHOLMOD_COL for column scaling - A * diag(scale) and CHOLMOD_SYM 1120ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // for symmetric scaling which scales both the rows and the columns 1130ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // - diag(scale) * A * diag(scale). 1140ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong void Scale(cholmod_dense* scale, int mode, cholmod_sparse* A) { 1150ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong cholmod_scale(scale, mode, A, &cc_); 1160ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong } 1170ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1180ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Create and return a matrix m = A * A'. Caller owns the 1190ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // result. The matrix A is not modified. 1200ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong cholmod_sparse* AATranspose(cholmod_sparse* A) { 1210ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong cholmod_sparse*m = cholmod_aat(A, NULL, A->nrow, 1, &cc_); 1220ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong m->stype = 1; // Pay attention to the upper triangular part. 1230ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong return m; 1240ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong } 1250ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1260ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // y = alpha * A * x + beta * y. Only y is modified. 1270ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong void SparseDenseMultiply(cholmod_sparse* A, double alpha, double beta, 1280ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong cholmod_dense* x, cholmod_dense* y) { 1290ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong double alpha_[2] = {alpha, 0}; 1300ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong double beta_[2] = {beta, 0}; 1310ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong cholmod_sdmult(A, 0, alpha_, beta_, x, y, &cc_); 1320ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong } 1330ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1340ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Find an ordering of A or AA' (if A is unsymmetric) that minimizes 1350ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // the fill-in in the Cholesky factorization of the corresponding 1360ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // matrix. This is done by using the AMD algorithm. 1370ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 1380ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Using this ordering, the symbolic Cholesky factorization of A (or 1390ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // AA') is computed and returned. 1400ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 1410ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // A is not modified, only the pattern of non-zeros of A is used, 1420ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // the actual numerical values in A are of no consequence. 1430ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 14479397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // message contains an explanation of the failures if any. 14579397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // 1460ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Caller owns the result. 14779397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez cholmod_factor* AnalyzeCholesky(cholmod_sparse* A, string* message); 1480ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1490ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong cholmod_factor* BlockAnalyzeCholesky(cholmod_sparse* A, 1500ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong const vector<int>& row_blocks, 15179397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez const vector<int>& col_blocks, 15279397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez string* message); 1530ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1540ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // If A is symmetric, then compute the symbolic Cholesky 1550ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // factorization of A(ordering, ordering). If A is unsymmetric, then 1560ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // compute the symbolic factorization of 1570ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // A(ordering,:) A(ordering,:)'. 1580ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 1590ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // A is not modified, only the pattern of non-zeros of A is used, 1600ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // the actual numerical values in A are of no consequence. 1610ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 16279397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // message contains an explanation of the failures if any. 16379397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // 1640ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Caller owns the result. 1650ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong cholmod_factor* AnalyzeCholeskyWithUserOrdering(cholmod_sparse* A, 16679397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez const vector<int>& ordering, 16779397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez string* message); 1680ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1691d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Perform a symbolic factorization of A without re-ordering A. No 1701d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // postordering of the elimination tree is performed. This ensures 1711d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // that the symbolic factor does not introduce an extra permutation 1721d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // on the matrix. See the documentation for CHOLMOD for more details. 17379397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // 17479397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // message contains an explanation of the failures if any. 17579397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez cholmod_factor* AnalyzeCholeskyWithNaturalOrdering(cholmod_sparse* A, 17679397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez string* message); 1771d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 1780ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Use the symbolic factorization in L, to find the numerical 1790ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // factorization for the matrix A or AA^T. Return true if 1800ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // successful, false otherwise. L contains the numeric factorization 1810ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // on return. 18279397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // 18379397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // message contains an explanation of the failures if any. 18479397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez LinearSolverTerminationType Cholesky(cholmod_sparse* A, 18579397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez cholmod_factor* L, 18679397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez string* message); 1870ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1880ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Given a Cholesky factorization of a matrix A = LL^T, solve the 1890ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // linear system Ax = b, and return the result. If the Solve fails 1900ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // NULL is returned. Caller owns the result. 19179397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // 19279397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // message contains an explanation of the failures if any. 19379397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez cholmod_dense* Solve(cholmod_factor* L, cholmod_dense* b, string* message); 1940ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1950ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // By virtue of the modeling layer in Ceres being block oriented, 1960ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // all the matrices used by Ceres are also block oriented. When 1970ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // doing sparse direct factorization of these matrices the 1980ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // fill-reducing ordering algorithms (in particular AMD) can either 1990ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // be run on the block or the scalar form of these matrices. The two 2000ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // SuiteSparse::AnalyzeCholesky methods allows the the client to 2010ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // compute the symbolic factorization of a matrix by either using 2020ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // AMD on the matrix or a user provided ordering of the rows. 2030ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2040ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // But since the underlying matrices are block oriented, it is worth 2050ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // running AMD on just the block structre of these matrices and then 2060ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // lifting these block orderings to a full scalar ordering. This 2070ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // preserves the block structure of the permuted matrix, and exposes 2080ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // more of the super-nodal structure of the matrix to the numerical 2090ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // factorization routines. 2100ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2110ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Find the block oriented AMD ordering of a matrix A, whose row and 2120ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // column blocks are given by row_blocks, and col_blocks 2130ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // respectively. The matrix may or may not be symmetric. The entries 2140ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // of col_blocks do not need to sum to the number of columns in 2150ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // A. If this is the case, only the first sum(col_blocks) are used 2160ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // to compute the ordering. 2170ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong bool BlockAMDOrdering(const cholmod_sparse* A, 2180ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong const vector<int>& row_blocks, 2190ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong const vector<int>& col_blocks, 2200ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong vector<int>* ordering); 2210ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2221d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Find a fill reducing approximate minimum degree 2231d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // ordering. ordering is expected to be large enough to hold the 2241d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // ordering. 22579397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez bool ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix, int* ordering); 2261d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 2271d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 2281d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Before SuiteSparse version 4.2.0, cholmod_camd was only enabled 2291d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // if SuiteSparse was compiled with Metis support. This makes 2301d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // calling and linking into cholmod_camd problematic even though it 2311d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // has nothing to do with Metis. This has been fixed reliably in 2321d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 4.2.0. 2331d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2341d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // The fix was actually committed in 4.1.0, but there is 2351d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // some confusion about a silent update to the tar ball, so we are 2361d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // being conservative and choosing the next minor version where 2371d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // things are stable. 2381d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling static bool IsConstrainedApproximateMinimumDegreeOrderingAvailable() { 2391d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling return (SUITESPARSE_VERSION>4001); 2401d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling } 2411d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 2421d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Find a fill reducing approximate minimum degree 2431d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // ordering. constraints is an array which associates with each 2441d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // column of the matrix an elimination group. i.e., all columns in 2451d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // group 0 are eliminated first, all columns in group 1 are 2461d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // eliminated next etc. This function finds a fill reducing ordering 2471d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // that obeys these constraints. 2481d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // 2491d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // Calling ApproximateMinimumDegreeOrdering is equivalent to calling 2501d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // ConstrainedApproximateMinimumDegreeOrdering with a constraint 2511d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // array that puts all columns in the same elimination group. 2520ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2531d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // If CERES_NO_CAMD is defined then calling this function will 2541d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // result in a crash. 25579397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez bool ConstrainedApproximateMinimumDegreeOrdering(cholmod_sparse* matrix, 2561d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling int* constraints, 2571d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling int* ordering); 2580ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2590ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong void Free(cholmod_sparse* m) { cholmod_free_sparse(&m, &cc_); } 2600ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong void Free(cholmod_dense* m) { cholmod_free_dense(&m, &cc_); } 2610ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong void Free(cholmod_factor* m) { cholmod_free_factor(&m, &cc_); } 2620ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2630ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong void Print(cholmod_sparse* m, const string& name) { 2640ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong cholmod_print_sparse(m, const_cast<char*>(name.c_str()), &cc_); 2650ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong } 2660ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2670ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong void Print(cholmod_dense* m, const string& name) { 2680ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong cholmod_print_dense(m, const_cast<char*>(name.c_str()), &cc_); 2690ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong } 2700ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2710ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong void Print(cholmod_triplet* m, const string& name) { 2720ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong cholmod_print_triplet(m, const_cast<char*>(name.c_str()), &cc_); 2730ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong } 2740ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2750ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong cholmod_common* mutable_cc() { return &cc_; } 2760ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2770ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong private: 2780ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong cholmod_common cc_; 2790ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong}; 2800ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2810ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong} // namespace internal 2820ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong} // namespace ceres 2830ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 284399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger#else // CERES_NO_SUITESPARSE 285399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger 286399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettingertypedef void cholmod_factor; 287399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger 28879397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandezclass SuiteSparse { 28979397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez public: 29079397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // Defining this static function even when SuiteSparse is not 29179397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // available, allows client code to check for the presence of CAMD 29279397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // without checking for the absence of the CERES_NO_CAMD symbol. 29379397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // 29479397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // This is safer because the symbol maybe missing due to a user 29579397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // accidently not including suitesparse.h in their code when 29679397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // checking for the symbol. 29779397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez static bool IsConstrainedApproximateMinimumDegreeOrderingAvailable() { 29879397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez return false; 29979397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez } 30079397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez 30179397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez void Free(void*) {}; 30279397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez}; 30379397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez 3040ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#endif // CERES_NO_SUITESPARSE 3050ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 3060ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#endif // CERES_INTERNAL_SUITESPARSE_H_ 307