diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc index 439270ed3..34f31ab6c 100644 --- a/internal/ceres/compressed_row_sparse_matrix.cc +++ b/internal/ceres/compressed_row_sparse_matrix.cc @@ -125,7 +125,7 @@ CompressedRowSparseMatrix::CompressedRowSparseMatrix( // Find the cumulative sum of the row counts. for (int i = 1; i < num_rows_ + 1; ++i) { - rows_[i] += rows_[i-1]; + rows_[i] += rows_[i - 1]; } CHECK_EQ(num_nonzeros(), m.num_nonzeros()); @@ -217,8 +217,8 @@ void CompressedRowSparseMatrix::DeleteRows(int delta_rows) { num_rows_ -= delta_rows; rows_.resize(num_rows_ + 1); - // Walk the list of row blocks untill we reach the new number of - // rows and then drop the rest of the row blocks. + // Walk the list of row blocks until we reach the new number of rows + // and the drop the rest of the row blocks. int num_row_blocks = 0; int num_rows = 0; while (num_row_blocks < row_blocks_.size() && num_rows < num_rows_) { @@ -380,6 +380,155 @@ CompressedRowSparseMatrix* CompressedRowSparseMatrix::Transpose() const { return transpose; } +namespace { +// A ProductTerm is a term in the outer product of a matrix with +// itself. +struct ProductTerm { + ProductTerm(const int row, const int col, const int index) + : row(row), col(col), index(index) { + } + + bool operator<(const ProductTerm& right) const { + if (row == right.row) { + if (col == right.col) { + return index < right.index; + } + return col < right.col; + } + return row < right.row; + } + + int row; + int col; + int index; +}; + +CompressedRowSparseMatrix* +CompressAndFillProgram(const int num_rows, + const int num_cols, + const vector& product, + vector* program) { + CHECK_GT(product.size(), 0); + + // Count the number of unique product term, which in turn is the + // number of non-zeros in the outer product. + int num_nonzeros = 1; + for (int i = 1; i < product.size(); ++i) { + if (product[i].row != product[i - 1].row || + product[i].col != product[i - 1].col) { + ++num_nonzeros; + } + } + + CompressedRowSparseMatrix* matrix = + new CompressedRowSparseMatrix(num_rows, num_cols, num_nonzeros); + + int* crsm_rows = matrix->mutable_rows(); + std::fill(crsm_rows, crsm_rows + num_rows + 1, 0); + int* crsm_cols = matrix->mutable_cols(); + std::fill(crsm_cols, crsm_cols + num_nonzeros, 0); + + CHECK_NOTNULL(program)->clear(); + program->resize(product.size()); + + // Iterate over the sorted product terms. This means each row is + // filled one at a time, and we are able to assign a position in the + // values array to each term. + // + // If terms repeat, i.e., they contribute to the same entry in the + // result matrix), then they do not affect the sparsity structure of + // the result matrix. + int nnz = 0; + crsm_cols[0] = product[0].col; + crsm_rows[product[0].row + 1]++; + (*program)[product[0].index] = nnz; + for (int i = 1; i < product.size(); ++i) { + const ProductTerm& previous = product[i - 1]; + const ProductTerm& current = product[i]; + + // Sparsity structure is updated only if the term is not a repeat. + if (previous.row != current.row || previous.col != current.col) { + crsm_cols[++nnz] = current.col; + crsm_rows[current.row + 1]++; + } + + // All terms get assigned the position in the values array where + // their value is accumulated. + (*program)[current.index] = nnz; + } + + for (int i = 1; i < num_rows + 1; ++i) { + crsm_rows[i] += crsm_rows[i - 1]; + } + + return matrix; +} + +} // namespace + +CompressedRowSparseMatrix* +CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram( + const CompressedRowSparseMatrix& m, + vector* program) { + CHECK_NOTNULL(program)->clear(); + CHECK_GT(m.num_nonzeros(), 0) << "Congratulations, " + << "you found a bug in Ceres. Please report it."; + + vector product; + const vector& row_blocks = m.row_blocks(); + int row_block_begin = 0; + // Iterate over row blocks + for (int row_block = 0; row_block < row_blocks.size(); ++row_block) { + const int row_block_end = row_block_begin + row_blocks[row_block]; + // Compute the outer product terms for just one row per row block. + const int r = row_block_begin; + // Compute the lower triangular part of the product. + for (int idx1 = m.rows()[r]; idx1 < m.rows()[r + 1]; ++idx1) { + for (int idx2 = m.rows()[r]; idx2 <= idx1; ++idx2) { + product.push_back(ProductTerm(m.cols()[idx1], m.cols()[idx2], product.size())); + } + } + row_block_begin = row_block_end; + } + CHECK_EQ(row_block_begin, m.num_rows()); + sort(product.begin(), product.end()); + return CompressAndFillProgram(m.num_cols(), m.num_cols(), product, program); +} + +void CompressedRowSparseMatrix::ComputeOuterProduct( + const CompressedRowSparseMatrix& m, + const vector& program, + CompressedRowSparseMatrix* result) { + result->SetZero(); + double* values = result->mutable_values(); + const vector& row_blocks = m.row_blocks(); + + int cursor = 0; + int row_block_begin = 0; + const double* m_values = m.values(); + const int* m_rows = m.rows(); + // Iterate over row blocks. + for (int row_block = 0; row_block < row_blocks.size(); ++row_block) { + const int row_block_end = row_block_begin + row_blocks[row_block]; + const int saved_cursor = cursor; + for (int r = row_block_begin; r < row_block_end; ++r) { + // Reuse the program segment for each row in this row block. + cursor = saved_cursor; + const int row_begin = m_rows[r]; + const int row_end = m_rows[r + 1]; + for (int idx1 = row_begin; idx1 < row_end; ++idx1) { + const double v1 = m_values[idx1]; + for (int idx2 = row_begin; idx2 <= idx1; ++idx2, ++cursor) { + values[program[cursor]] += v1 * m_values[idx2]; + } + } + } + row_block_begin = row_block_end; + } + + CHECK_EQ(row_block_begin, m.num_rows()); + CHECK_EQ(cursor, program.size()); +} } // namespace internal } // namespace ceres diff --git a/internal/ceres/compressed_row_sparse_matrix.h b/internal/ceres/compressed_row_sparse_matrix.h index c5721eb88..06b86896f 100644 --- a/internal/ceres/compressed_row_sparse_matrix.h +++ b/internal/ceres/compressed_row_sparse_matrix.h @@ -128,6 +128,32 @@ class CompressedRowSparseMatrix : public SparseMatrix { const double* diagonal, const vector& blocks); + // Compute the sparsity structure of the product m.transpose() * m + // and create a CompressedRowSparseMatrix corresponding to it. + // + // Also compute a "program" vector, which for every term in the + // outer product points to the entry in the values array of the + // result matrix where it should be accumulated. + // + // This program is used by the ComputeOuterProduct function below to + // compute the outer product. + // + // Since the entries of the program are the same for rows with the + // same sparsity structure, the program only stores the result for + // one row per row block. The ComputeOuterProduct function reuses + // this information for each row in the row block. + static CompressedRowSparseMatrix* CreateOuterProductMatrixAndProgram( + const CompressedRowSparseMatrix& m, + vector* program); + + // Compute the values array for the expression m.transpose() * m, + // where the matrix used to store the result and a program have been + // created using the CreateOuterProductMatrixAndProgram function + // above. + static void ComputeOuterProduct(const CompressedRowSparseMatrix& m, + const vector& program, + CompressedRowSparseMatrix* result); + private: int num_rows_; int num_cols_; diff --git a/internal/ceres/compressed_row_sparse_matrix_test.cc b/internal/ceres/compressed_row_sparse_matrix_test.cc index 4d20083a3..ea737b8d0 100644 --- a/internal/ceres/compressed_row_sparse_matrix_test.cc +++ b/internal/ceres/compressed_row_sparse_matrix_test.cc @@ -30,11 +30,14 @@ #include "ceres/compressed_row_sparse_matrix.h" +#include #include "ceres/casts.h" #include "ceres/crs_matrix.h" +#include "ceres/cxsparse.h" #include "ceres/internal/eigen.h" #include "ceres/internal/scoped_ptr.h" #include "ceres/linear_least_squares_problems.h" +#include "ceres/random.h" #include "ceres/triplet_sparse_matrix.h" #include "glog/logging.h" #include "gtest/gtest.h" @@ -381,5 +384,170 @@ TEST(CompressedRowSparseMatrix, Transpose) { EXPECT_NEAR((dense_matrix - dense_transpose.transpose()).norm(), 0.0, 1e-14); } +#ifndef CERES_NO_CXSPARSE + +struct RandomMatrixOptions { + int num_row_blocks; + int min_row_block_size; + int max_row_block_size; + int num_col_blocks; + int min_col_block_size; + int max_col_block_size; + double block_density; +}; + +CompressedRowSparseMatrix* CreateRandomCompressedRowSparseMatrix( + const RandomMatrixOptions& options) { + vector row_blocks; + for (int i = 0; i < options.num_row_blocks; ++i) { + const int delta_block_size = + Uniform(options.max_row_block_size - options.min_row_block_size); + row_blocks.push_back(options.min_row_block_size + delta_block_size); + } + + vector col_blocks; + for (int i = 0; i < options.num_col_blocks; ++i) { + const int delta_block_size = + Uniform(options.max_col_block_size - options.min_col_block_size); + col_blocks.push_back(options.min_col_block_size + delta_block_size); + } + + vector rows; + vector cols; + vector values; + + while (values.size() == 0) { + int row_block_begin = 0; + for (int r = 0; r < options.num_row_blocks; ++r) { + int col_block_begin = 0; + for (int c = 0; c < options.num_col_blocks; ++c) { + if (RandDouble() <= options.block_density) { + for (int i = 0; i < row_blocks[r]; ++i) { + for (int j = 0; j < col_blocks[c]; ++j) { + rows.push_back(row_block_begin + i); + cols.push_back(col_block_begin + j); + values.push_back(RandNormal()); + } + } + } + col_block_begin += col_blocks[c]; + } + row_block_begin += row_blocks[r]; + } + } + + const int num_rows = std::accumulate(row_blocks.begin(), row_blocks.end(), 0); + const int num_cols = std::accumulate(col_blocks.begin(), col_blocks.end(), 0); + const int num_nonzeros = values.size(); + + TripletSparseMatrix tsm(num_rows, num_cols, num_nonzeros); + std::copy(rows.begin(), rows.end(), tsm.mutable_rows()); + std::copy(cols.begin(), cols.end(), tsm.mutable_cols()); + std::copy(values.begin(), values.end(), tsm.mutable_values()); + tsm.set_num_nonzeros(num_nonzeros); + CompressedRowSparseMatrix* matrix = new CompressedRowSparseMatrix(tsm); + (*matrix->mutable_row_blocks()) = row_blocks; + (*matrix->mutable_col_blocks()) = col_blocks; + return matrix; +} + +void ToDenseMatrix(const cs_di* matrix, Matrix* dense_matrix) { + dense_matrix->resize(matrix->m, matrix->n); + dense_matrix->setZero(); + + for (int c = 0; c < matrix->n; ++c) { + for (int idx = matrix->p[c]; idx < matrix->p[c + 1]; ++idx) { + const int r = matrix->i[idx]; + (*dense_matrix)(r, c) = matrix->x[idx]; + } + } +} + +TEST(CompressedRowSparseMatrix, ComputeOuterProduct) { + // "Randomly generated seed." + SetRandomState(29823); + int kMaxNumRowBlocks = 10; + int kMaxNumColBlocks = 10; + int kNumTrials = 10; + + CXSparse cxsparse; + const double kTolerance = 1e-18; + + // Create a random matrix, compute its outer product using CXSParse + // and ComputeOuterProduct. Convert both matrices to dense matrices + // and compare their upper triangular parts. They should be within + // kTolerance of each other. + for (int num_row_blocks = 1; + num_row_blocks < kMaxNumRowBlocks; + ++num_row_blocks) { + for (int num_col_blocks = 1; + num_col_blocks < kMaxNumColBlocks; + ++num_col_blocks) { + for (int trial = 0; trial < kNumTrials; ++trial) { + + + RandomMatrixOptions options; + options.num_row_blocks = num_row_blocks; + options.num_col_blocks = num_col_blocks; + options.min_row_block_size = 1; + options.max_row_block_size = 5; + options.min_col_block_size = 1; + options.max_col_block_size = 10; + options.block_density = std::max(0.1, RandDouble()); + + VLOG(2) << "num row blocks: " << options.num_row_blocks; + VLOG(2) << "num col blocks: " << options.num_col_blocks; + VLOG(2) << "min row block size: " << options.min_row_block_size; + VLOG(2) << "max row block size: " << options.max_row_block_size; + VLOG(2) << "min col block size: " << options.min_col_block_size; + VLOG(2) << "max col block size: " << options.max_col_block_size; + VLOG(2) << "block density: " << options.block_density; + + scoped_ptr matrix( + CreateRandomCompressedRowSparseMatrix(options)); + + cs_di cs_matrix_transpose = cxsparse.CreateSparseMatrixTransposeView(matrix.get()); + cs_di* cs_matrix = cxsparse.TransposeMatrix(&cs_matrix_transpose); + cs_di* expected_outer_product = + cxsparse.MatrixMatrixMultiply(&cs_matrix_transpose, cs_matrix); + + vector program; + scoped_ptr outer_product( + CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram( + *matrix, &program)); + CompressedRowSparseMatrix::ComputeOuterProduct(*matrix, + program, + outer_product.get()); + + cs_di actual_outer_product = + cxsparse.CreateSparseMatrixTransposeView(outer_product.get()); + + ASSERT_EQ(actual_outer_product.m, actual_outer_product.n); + ASSERT_EQ(expected_outer_product->m, expected_outer_product->n); + ASSERT_EQ(actual_outer_product.m, expected_outer_product->m); + + Matrix actual_matrix; + Matrix expected_matrix; + + ToDenseMatrix(expected_outer_product, &expected_matrix); + expected_matrix.triangularView().setZero(); + + ToDenseMatrix(&actual_outer_product, &actual_matrix); + const double diff_norm = (actual_matrix - expected_matrix).norm() / expected_matrix.norm(); + ASSERT_NEAR(diff_norm, 0.0, kTolerance) + << "expected: \n" + << expected_matrix + << "\nactual: \n" + << actual_matrix; + + cxsparse.Free(cs_matrix); + cxsparse.Free(expected_outer_product); + } + } + } +} + +#endif // CERES_NO_CXSPARSE + } // namespace internal } // namespace ceres diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc index b2dba440e..880adf7eb 100644 --- a/internal/ceres/sparse_normal_cholesky_solver.cc +++ b/internal/ceres/sparse_normal_cholesky_solver.cc @@ -77,37 +77,10 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImpl( const double* b, const LinearSolver::PerSolveOptions& per_solve_options, double * x) { - switch (options_.sparse_linear_algebra_library_type) { - case SUITE_SPARSE: - return SolveImplUsingSuiteSparse(A, b, per_solve_options, x); - case CX_SPARSE: - return SolveImplUsingCXSparse(A, b, per_solve_options, x); - default: - LOG(FATAL) << "Unknown sparse linear algebra library : " - << options_.sparse_linear_algebra_library_type; - } - - LOG(FATAL) << "Unknown sparse linear algebra library : " - << options_.sparse_linear_algebra_library_type; - return LinearSolver::Summary(); -} - -#ifndef CERES_NO_CXSPARSE -LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( - CompressedRowSparseMatrix* A, - const double* b, - const LinearSolver::PerSolveOptions& per_solve_options, - double * x) { - EventLogger event_logger("SparseNormalCholeskySolver::CXSparse::Solve"); - - LinearSolver::Summary summary; - summary.num_iterations = 1; - summary.termination_type = LINEAR_SOLVER_SUCCESS; - summary.message = "Success."; const int num_cols = A->num_cols(); - Vector Atb = Vector::Zero(num_cols); - A->LeftMultiply(b, Atb.data()); + VectorRef(x, num_cols).setZero(); + A->LeftMultiply(b, x); if (per_solve_options.D != NULL) { // Temporarily append a diagonal block to the A matrix, but undo @@ -123,10 +96,37 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( A->AppendRows(*regularizer); } - VectorRef(x, num_cols).setZero(); + LinearSolver::Summary summary; + switch (options_.sparse_linear_algebra_library_type) { + case SUITE_SPARSE: + summary = SolveImplUsingSuiteSparse(A, per_solve_options, x); + break; + case CX_SPARSE: + summary = SolveImplUsingCXSparse(A, per_solve_options, x); + break; + default: + LOG(FATAL) << "Unknown sparse linear algebra library : " + << options_.sparse_linear_algebra_library_type; + } - // Wrap the augmented Jacobian in a compressed sparse column matrix. - cs_di At = cxsparse_.CreateSparseMatrixTransposeView(A); + if (per_solve_options.D != NULL) { + A->DeleteRows(num_cols); + } + + return summary; +} + +#ifndef CERES_NO_CXSPARSE +LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( + CompressedRowSparseMatrix* A, + const LinearSolver::PerSolveOptions& per_solve_options, + double * rhs_and_solution) { + EventLogger event_logger("SparseNormalCholeskySolver::CXSparse::Solve"); + + LinearSolver::Summary summary; + summary.num_iterations = 1; + summary.termination_type = LINEAR_SOLVER_SUCCESS; + summary.message = "Success."; // Compute the normal equations. J'J delta = J'f and solve them // using a sparse Cholesky factorization. Notice that when compared @@ -135,13 +135,18 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( // factorized. CHOLMOD/SuiteSparse on the other hand can just work // off of Jt to compute the Cholesky factorization of the normal // equations. - cs_di* A2 = cxsparse_.TransposeMatrix(&At); - cs_di* AtA = cxsparse_.MatrixMatrixMultiply(&At, A2); - - cxsparse_.Free(A2); - if (per_solve_options.D != NULL) { - A->DeleteRows(num_cols); + if (outer_product_.get() == NULL) { + outer_product_.reset( + CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram( + *A, &pattern_)); } + + CompressedRowSparseMatrix::ComputeOuterProduct( + *A, pattern_, outer_product_.get()); + cs_di AtA_view = + cxsparse_.CreateSparseMatrixTransposeView(outer_product_.get()); + cs_di* AtA = &AtA_view; + event_logger.AddEvent("Setup"); // Compute symbolic factorization if not available. @@ -160,23 +165,17 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( summary.termination_type = LINEAR_SOLVER_FATAL_ERROR; summary.message = "CXSparse failure. Unable to find symbolic factorization."; - } else if (cxsparse_.SolveCholesky(AtA, cxsparse_factor_, Atb.data())) { - VectorRef(x, Atb.rows()) = Atb; - } else { + } else if (!cxsparse_.SolveCholesky(AtA, cxsparse_factor_, rhs_and_solution)) { summary.termination_type = LINEAR_SOLVER_FAILURE; } event_logger.AddEvent("Solve"); - - cxsparse_.Free(AtA); - event_logger.AddEvent("Teardown"); return summary; } #else LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( CompressedRowSparseMatrix* A, - const double* b, const LinearSolver::PerSolveOptions& per_solve_options, - double * x) { + double * rhs_and_solution) { LOG(FATAL) << "No CXSparse support in Ceres."; // Unreachable but MSVC does not know this. @@ -187,9 +186,8 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingCXSparse( #ifndef CERES_NO_SUITESPARSE LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( CompressedRowSparseMatrix* A, - const double* b, const LinearSolver::PerSolveOptions& per_solve_options, - double * x) { + double * rhs_and_solution) { EventLogger event_logger("SparseNormalCholeskySolver::SuiteSparse::Solve"); LinearSolver::Summary summary; summary.termination_type = LINEAR_SOLVER_SUCCESS; @@ -197,24 +195,6 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( summary.message = "Success."; const int num_cols = A->num_cols(); - Vector Atb = Vector::Zero(num_cols); - A->LeftMultiply(b, Atb.data()); - - if (per_solve_options.D != NULL) { - // Temporarily append a diagonal block to the A matrix, but undo - // it before returning the matrix to the user. - scoped_ptr regularizer; - if (A->col_blocks().size() > 0) { - regularizer.reset(CompressedRowSparseMatrix::CreateBlockDiagonalMatrix( - per_solve_options.D, A->col_blocks())); - } else { - regularizer.reset(new CompressedRowSparseMatrix( - per_solve_options.D, num_cols)); - } - A->AppendRows(*regularizer); - } - - VectorRef(x, num_cols).setZero(); cholmod_sparse lhs = ss_.CreateSparseMatrixTransposeView(A); event_logger.AddEvent("Setup"); @@ -231,33 +211,23 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( event_logger.AddEvent("Analysis"); if (factor_ == NULL) { - if (per_solve_options.D != NULL) { - A->DeleteRows(num_cols); - } summary.termination_type = LINEAR_SOLVER_FATAL_ERROR; return summary; } summary.termination_type = ss_.Cholesky(&lhs, factor_, &summary.message); if (summary.termination_type != LINEAR_SOLVER_SUCCESS) { - if (per_solve_options.D != NULL) { - A->DeleteRows(num_cols); - } return summary; } - cholmod_dense* rhs = ss_.CreateDenseVector(Atb.data(), num_cols, num_cols); - cholmod_dense* sol = ss_.Solve(factor_, rhs, &summary.message); + cholmod_dense* rhs = ss_.CreateDenseVector(rhs_and_solution, num_cols, num_cols); + cholmod_dense* solution = ss_.Solve(factor_, rhs, &summary.message); event_logger.AddEvent("Solve"); ss_.Free(rhs); - if (per_solve_options.D != NULL) { - A->DeleteRows(num_cols); - } - - if (sol != NULL) { - memcpy(x, sol->x, num_cols * sizeof(*x)); - ss_.Free(sol); + if (solution != NULL) { + memcpy(rhs_and_solution, solution->x, num_cols * sizeof(*rhs_and_solution)); + ss_.Free(solution); } else { summary.termination_type = LINEAR_SOLVER_FAILURE; } @@ -268,9 +238,8 @@ LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( #else LinearSolver::Summary SparseNormalCholeskySolver::SolveImplUsingSuiteSparse( CompressedRowSparseMatrix* A, - const double* b, const LinearSolver::PerSolveOptions& per_solve_options, - double * x) { + double * rhs_and_solution) { LOG(FATAL) << "No SuiteSparse support in Ceres."; // Unreachable but MSVC does not know this. diff --git a/internal/ceres/sparse_normal_cholesky_solver.h b/internal/ceres/sparse_normal_cholesky_solver.h index 61111b41b..adca08a50 100644 --- a/internal/ceres/sparse_normal_cholesky_solver.h +++ b/internal/ceres/sparse_normal_cholesky_solver.h @@ -62,16 +62,14 @@ class SparseNormalCholeskySolver : public CompressedRowSparseMatrixSolver { LinearSolver::Summary SolveImplUsingSuiteSparse( CompressedRowSparseMatrix* A, - const double* b, const LinearSolver::PerSolveOptions& options, - double* x); + double* rhs_and_solution); // Crashes if CSparse is not installed. LinearSolver::Summary SolveImplUsingCXSparse( CompressedRowSparseMatrix* A, - const double* b, const LinearSolver::PerSolveOptions& options, - double* x); + double* rhs_and_solution); SuiteSparse ss_; // Cached factorization @@ -80,7 +78,8 @@ class SparseNormalCholeskySolver : public CompressedRowSparseMatrixSolver { CXSparse cxsparse_; // Cached factorization cs_dis* cxsparse_factor_; - + scoped_ptr outer_product_; + vector pattern_; const LinearSolver::Options options_; CERES_DISALLOW_COPY_AND_ASSIGN(SparseNormalCholeskySolver); }; diff --git a/internal/ceres/unsymmetric_linear_solver_test.cc b/internal/ceres/unsymmetric_linear_solver_test.cc index 949fdd1f5..624db68bb 100644 --- a/internal/ceres/unsymmetric_linear_solver_test.cc +++ b/internal/ceres/unsymmetric_linear_solver_test.cc @@ -57,7 +57,7 @@ class UnsymmetricLinearSolverTest : public ::testing::Test { } void TestSolver(const LinearSolver::Options& options) { - scoped_ptr solver(LinearSolver::Create(options)); + LinearSolver::PerSolveOptions per_solve_options; LinearSolver::Summary unregularized_solve_summary; @@ -84,13 +84,17 @@ class UnsymmetricLinearSolverTest : public ::testing::Test { } else { LOG(FATAL) << "Unknown linear solver : " << options.type; } + // Unregularized + scoped_ptr solver(LinearSolver::Create(options)); unregularized_solve_summary = solver->Solve(transformed_A.get(), b_.get(), per_solve_options, x_unregularized.data()); + // Sparsity structure is changing, reset the solver. + solver.reset(LinearSolver::Create(options)); // Regularized solution per_solve_options.D = D_.get(); regularized_solve_summary =