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#ifndef CERES_INTERNAL_COMPRESSED_ROW_SPARSE_MATRIX_H_
32#define CERES_INTERNAL_COMPRESSED_ROW_SPARSE_MATRIX_H_
33
34#include <vector>
35#include "ceres/internal/macros.h"
36#include "ceres/internal/port.h"
37#include "ceres/sparse_matrix.h"
38#include "ceres/types.h"
39#include "glog/logging.h"
40
41namespace ceres {
42
43struct CRSMatrix;
44
45namespace internal {
46
47class TripletSparseMatrix;
48
49class CompressedRowSparseMatrix : public SparseMatrix {
50 public:
51  // Build a matrix with the same content as the TripletSparseMatrix
52  // m. TripletSparseMatrix objects are easier to construct
53  // incrementally, so we use them to initialize SparseMatrix
54  // objects.
55  //
56  // We assume that m does not have any repeated entries.
57  explicit CompressedRowSparseMatrix(const TripletSparseMatrix& m);
58
59  // Use this constructor only if you know what you are doing. This
60  // creates a "blank" matrix with the appropriate amount of memory
61  // allocated. However, the object itself is in an inconsistent state
62  // as the rows and cols matrices do not match the values of
63  // num_rows, num_cols and max_num_nonzeros.
64  //
65  // The use case for this constructor is that when the user knows the
66  // size of the matrix to begin with and wants to update the layout
67  // manually, instead of going via the indirect route of first
68  // constructing a TripletSparseMatrix, which leads to more than
69  // double the peak memory usage.
70  CompressedRowSparseMatrix(int num_rows,
71                            int num_cols,
72                            int max_num_nonzeros);
73
74  // Build a square sparse diagonal matrix with num_rows rows and
75  // columns. The diagonal m(i,i) = diagonal(i);
76  CompressedRowSparseMatrix(const double* diagonal, int num_rows);
77
78  virtual ~CompressedRowSparseMatrix();
79
80  // SparseMatrix interface.
81  virtual void SetZero();
82  virtual void RightMultiply(const double* x, double* y) const;
83  virtual void LeftMultiply(const double* x, double* y) const;
84  virtual void SquaredColumnNorm(double* x) const;
85  virtual void ScaleColumns(const double* scale);
86
87  virtual void ToDenseMatrix(Matrix* dense_matrix) const;
88  virtual void ToTextFile(FILE* file) const;
89  virtual int num_rows() const { return num_rows_; }
90  virtual int num_cols() const { return num_cols_; }
91  virtual int num_nonzeros() const { return rows_[num_rows_]; }
92  virtual const double* values() const { return &values_[0]; }
93  virtual double* mutable_values() { return &values_[0]; }
94
95  // Delete the bottom delta_rows.
96  // num_rows -= delta_rows
97  void DeleteRows(int delta_rows);
98
99  // Append the contents of m to the bottom of this matrix. m must
100  // have the same number of columns as this matrix.
101  void AppendRows(const CompressedRowSparseMatrix& m);
102
103  void ToCRSMatrix(CRSMatrix* matrix) const;
104
105  // Low level access methods that expose the structure of the matrix.
106  const int* cols() const { return &cols_[0]; }
107  int* mutable_cols() { return &cols_[0]; }
108
109  const int* rows() const { return &rows_[0]; }
110  int* mutable_rows() { return &rows_[0]; }
111
112  const vector<int>& row_blocks() const { return row_blocks_; }
113  vector<int>* mutable_row_blocks() { return &row_blocks_; }
114
115  const vector<int>& col_blocks() const { return col_blocks_; }
116  vector<int>* mutable_col_blocks() { return &col_blocks_; }
117
118  // Destructive array resizing method.
119  void SetMaxNumNonZeros(int num_nonzeros);
120
121  // Non-destructive array resizing method.
122  void set_num_rows(const int num_rows) { num_rows_ = num_rows; }
123  void set_num_cols(const int num_cols) { num_cols_ = num_cols; }
124
125  void SolveLowerTriangularInPlace(double* solution) const;
126  void SolveLowerTriangularTransposeInPlace(double* solution) const;
127
128  CompressedRowSparseMatrix* Transpose() const;
129
130  static CompressedRowSparseMatrix* CreateBlockDiagonalMatrix(
131      const double* diagonal,
132      const vector<int>& blocks);
133
134  // Compute the sparsity structure of the product m.transpose() * m
135  // and create a CompressedRowSparseMatrix corresponding to it.
136  //
137  // Also compute a "program" vector, which for every term in the
138  // outer product points to the entry in the values array of the
139  // result matrix where it should be accumulated.
140  //
141  // This program is used by the ComputeOuterProduct function below to
142  // compute the outer product.
143  //
144  // Since the entries of the program are the same for rows with the
145  // same sparsity structure, the program only stores the result for
146  // one row per row block. The ComputeOuterProduct function reuses
147  // this information for each row in the row block.
148  static CompressedRowSparseMatrix* CreateOuterProductMatrixAndProgram(
149      const CompressedRowSparseMatrix& m,
150      vector<int>* program);
151
152  // Compute the values array for the expression m.transpose() * m,
153  // where the matrix used to store the result and a program have been
154  // created using the CreateOuterProductMatrixAndProgram function
155  // above.
156  static void ComputeOuterProduct(const CompressedRowSparseMatrix& m,
157                                  const vector<int>& program,
158                                  CompressedRowSparseMatrix* result);
159
160 private:
161  int num_rows_;
162  int num_cols_;
163  vector<int> rows_;
164  vector<int> cols_;
165  vector<double> values_;
166
167  // If the matrix has an underlying block structure, then it can also
168  // carry with it row and column block sizes. This is auxilliary and
169  // optional information for use by algorithms operating on the
170  // matrix. The class itself does not make use of this information in
171  // any way.
172  vector<int> row_blocks_;
173  vector<int> col_blocks_;
174
175  CERES_DISALLOW_COPY_AND_ASSIGN(CompressedRowSparseMatrix);
176};
177
178}  // namespace internal
179}  // namespace ceres
180
181#endif  // CERES_INTERNAL_COMPRESSED_ROW_SPARSE_MATRIX_H_
182