Point Cloud Library (PCL)  1.9.1-dev
sparse_matrix.h
1 /*
2 Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
3 All rights reserved.
4 
5 Redistribution and use in source and binary forms, with or without modification,
6 are permitted provided that the following conditions are met:
7 
8 Redistributions of source code must retain the above copyright notice, this list of
9 conditions and the following disclaimer. Redistributions in binary form must reproduce
10 the above copyright notice, this list of conditions and the following disclaimer
11 in the documentation and/or other materials provided with the distribution.
12 
13 Neither the name of the Johns Hopkins University nor the names of its contributors
14 may be used to endorse or promote products derived from this software without specific
15 prior written permission.
16 
17 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
18 EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
19 OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
20 SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
21 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
22 TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
23 BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
25 ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
26 DAMAGE.
27 */
28 
29 #ifndef __SPARSEMATRIX_HPP
30 #define __SPARSEMATRIX_HPP
31 
32 #if defined __GNUC__
33 # pragma GCC system_header
34 #endif
35 
36 #include "vector.h"
37 #include "allocator.h"
38 
39 namespace pcl
40 {
41  namespace poisson
42  {
43 
44  template <class T>
45  struct MatrixEntry
46  {
47  MatrixEntry( void ) { N =-1; Value = 0; }
48  MatrixEntry( int i ) { N = i; Value = 0; }
49  int N;
50  T Value;
51  };
52 
53  template<class T> class SparseMatrix
54  {
55  private:
56  bool _contiguous;
57  int _maxEntriesPerRow;
58  static int UseAlloc;
59  public:
61  static int UseAllocator(void);
62  static void SetAllocator( int blockSize );
63 
64  int rows;
65  int* rowSizes;
67  MatrixEntry< T >* operator[] ( int idx ) { return m_ppElements[idx]; }
68  const MatrixEntry< T >* operator[] ( int idx ) const { return m_ppElements[idx]; }
69 
70  SparseMatrix( void );
71  SparseMatrix( int rows );
72  SparseMatrix( int rows , int maxEntriesPerRow );
73  void Resize( int rows );
74  void Resize( int rows , int maxEntriesPerRow );
75  void SetRowSize( int row , int count );
76  int Entries( void ) const;
77 
78  SparseMatrix( const SparseMatrix& M );
79  ~SparseMatrix();
80 
81  void SetZero();
82  void SetIdentity();
83 
84  SparseMatrix<T>& operator = (const SparseMatrix<T>& M);
85 
86  SparseMatrix<T> operator * (const T& V) const;
87  SparseMatrix<T>& operator *= (const T& V);
88 
89 
90  SparseMatrix<T> operator * (const SparseMatrix<T>& M) const;
91  SparseMatrix<T> Multiply( const SparseMatrix<T>& M ) const;
92  SparseMatrix<T> MultiplyTranspose( const SparseMatrix<T>& Mt ) const;
93 
94  template<class T2>
95  Vector<T2> operator * (const Vector<T2>& V) const;
96  template<class T2>
97  Vector<T2> Multiply( const Vector<T2>& V ) const;
98  template<class T2>
99  void Multiply( const Vector<T2>& In , Vector<T2>& Out , int threads=1 ) const;
100 
101 
102  SparseMatrix<T> Transpose() const;
103 
104  static int Solve (const SparseMatrix<T>& M,const Vector<T>& b, int iters,Vector<T>& solution,const T eps=1e-8);
105 
106  template<class T2>
107  static int SolveSymmetric( const SparseMatrix<T>& M , const Vector<T2>& b , int iters , Vector<T2>& solution , const T2 eps=1e-8 , int reset=1 , int threads=1 );
108 
109  bool write( FILE* fp ) const;
110  bool write( const char* fileName ) const;
111  bool read( FILE* fp );
112  bool read( const char* fileName );
113  };
114 
115  template< class T2 >
117  {
118  private:
119  int _dim;
120  public:
121  std::vector< T2* > out;
122  MapReduceVector( void ) { _dim = 0; }
124  {
125  if( _dim ) for( int t=0 ; t<int(out.size()) ; t++ ) delete[] out[t];
126  out.resize( 0 );
127  }
128  T2* operator[]( int t ) { return out[t]; }
129  const T2* operator[]( int t ) const { return out[t]; }
130  int threads( void ) const { return int( out.size() ); }
131  void resize( int threads , int dim )
132  {
133  if( threads!=out.size() || _dim<dim )
134  {
135  for( int t=0 ; t<int(out.size()) ; t++ ) delete[] out[t];
136  out.resize( threads );
137  for( int t=0 ; t<int(out.size()) ; t++ ) out[t] = new T2[dim];
138  _dim = dim;
139  }
140  }
141 
142  };
143 
144  template< class T >
146  {
147  public:
148 
149  template< class T2 >
150  Vector< T2 > operator * ( const Vector<T2>& V ) const;
151 
152  template< class T2 >
153  Vector< T2 > Multiply( const Vector<T2>& V ) const;
154 
155  template< class T2 >
156  void Multiply( const Vector<T2>& In, Vector<T2>& Out , bool addDCTerm=false ) const;
157 
158  template< class T2 >
159  void Multiply( const Vector<T2>& In, Vector<T2>& Out , MapReduceVector< T2 >& OutScratch , bool addDCTerm=false ) const;
160 
161  template< class T2 >
162  void Multiply( const Vector<T2>& In, Vector<T2>& Out , std::vector< T2* >& OutScratch , const std::vector< int >& bounds ) const;
163 
164  template< class T2 >
165  static int Solve( const SparseSymmetricMatrix<T>& M , const Vector<T2>& b , int iters , Vector<T2>& solution , T2 eps=1e-8 , int reset=1 , int threads=0 , bool addDCTerm=false , bool solveNormal=false );
166 
167  template< class T2 >
168  static int Solve( const SparseSymmetricMatrix<T>& M , const Vector<T2>& b , int iters , Vector<T2>& solution , MapReduceVector<T2>& scratch , T2 eps=1e-8 , int reset=1 , bool addDCTerm=false , bool solveNormal=false );
169 #if defined _WIN32 && !defined __MINGW32__
170  template< class T2 >
171  static int SolveAtomic( const SparseSymmetricMatrix<T>& M , const Vector<T2>& b , int iters , Vector<T2>& solution , T2 eps=1e-8 , int reset=1 , int threads=0 , bool solveNormal=false );
172 #endif // _WIN32 || __MINGW32__
173  template<class T2>
174  static int Solve( const SparseSymmetricMatrix<T>& M , const Vector<T2>& diagonal , const Vector<T2>& b , int iters , Vector<T2>& solution , int reset=1 );
175 
176  template< class T2 >
177  void getDiagonal( Vector< T2 >& diagonal ) const;
178  };
179 
180 
181  }
182 }
183 
184 
185 #include "sparse_matrix.hpp"
186 
187 #endif
188 
MatrixEntry(int i)
Definition: sparse_matrix.h:48
static Allocator< MatrixEntry< T > > internalAllocator
Definition: sparse_matrix.h:60
This file defines compatibility wrappers for low level I/O functions.
Definition: convolution.h:45
Definition: sparse_matrix.h:45
void read(std::istream &stream, Type &value)
Function for reading data from a stream.
Definition: region_xy.h:46
PCL_EXPORTS void Multiply(const double in1[2], const double in2[2], double out[2])
T Value
Definition: sparse_matrix.h:50
int N
Definition: sparse_matrix.h:49
const T2 * operator[](int t) const
PCL_EXPORTS int Solve(const double *eqns, const double *values, double *solutions, int dim)
MatrixEntry< T > ** m_ppElements
Definition: sparse_matrix.h:66
void resize(int threads, int dim)
MatrixEntry(void)
Definition: sparse_matrix.h:47
This templated class assists in memory allocation and is well suited for instances when it is known t...
Definition: allocator.h:50
void write(std::ostream &stream, Type value)
Function for writing data to a stream.
Definition: region_xy.h:63