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 #define NEW_SPARSE_MATRIX 1 33 #define ZERO_TESTING_JACOBI 1 34 35 36 #include "Array.h" 37 38 template <class T> 39 struct MatrixEntry 40 { MatrixEntryMatrixEntry41 MatrixEntry( void ) { N =-1; Value = 0; } MatrixEntryMatrixEntry42 MatrixEntry( int i ) { N = i; Value = 0; } MatrixEntryMatrixEntry43 MatrixEntry( int i , T v ) { N = i; Value = v; } 44 int N; 45 T Value; 46 }; 47 48 template<class T> class SparseMatrix 49 { 50 private: 51 bool _contiguous; 52 int _maxEntriesPerRow; 53 void _init( void ); 54 public: 55 int rows; 56 Pointer( int ) rowSizes; 57 Pointer( Pointer( MatrixEntry< T > ) ) m_ppElements; Pointer(MatrixEntry<T>)58 Pointer( MatrixEntry< T > ) operator[] ( int idx ) { return m_ppElements[idx]; } ConstPointer(MatrixEntry<T>)59 ConstPointer( MatrixEntry< T > ) operator[] ( int idx ) const { return m_ppElements[idx]; } 60 61 SparseMatrix( void ); 62 SparseMatrix( int rows ); 63 SparseMatrix( int rows , int maxEntriesPerRow ); 64 void Resize( int rows ); 65 void Resize( int rows , int maxEntriesPerRow ); 66 void SetRowSize( int row , int count ); 67 int Entries( void ) const; 68 69 SparseMatrix( const SparseMatrix& M ); 70 ~SparseMatrix(); 71 72 void SetZero(); 73 74 SparseMatrix<T>& operator = (const SparseMatrix<T>& M); 75 76 SparseMatrix<T> operator * (const T& V) const; 77 SparseMatrix<T>& operator *= (const T& V); 78 79 template< class T2 > void Multiply( ConstPointer( T2 ) in , Pointer( T2 ) out , int threads=1 ) const; 80 template< class T2 > void MultiplyAndAddAverage( ConstPointer( T2 ) in , Pointer( T2 ) out , int threads=1 ) const; 81 82 bool write( FILE* fp ) const; 83 bool write( const char* fileName ) const; 84 bool read( FILE* fp ); 85 bool read( const char* fileName ); 86 87 template< class T2 > void getDiagonal( Pointer( T2 ) diagonal , int threads=1 ) const; 88 template< class T2 > static int SolveJacobi( const SparseMatrix<T>& M , ConstPointer( T2 ) b , Pointer( T2 ) x , Pointer( T2 ) Mx , T2 sor , int threads=1 ); 89 template< class T2 > static int SolveJacobi( const SparseMatrix<T>& M , ConstPointer( T2 ) diagonal , ConstPointer( T2 ) b , Pointer( T2 ) x , Pointer( T2 ) Mx , T2 sor , int threads=1 ); 90 template< class T2 > static int SolveGS( const SparseMatrix<T>& M , ConstPointer( T2 ) b , Pointer( T2 ) x , bool forward ); 91 template< class T2 > static int SolveGS( const SparseMatrix<T>& M , ConstPointer( T2 ) diagonal , ConstPointer( T2 ) b , Pointer( T2 ) x , bool forward ); 92 template< class T2 > static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseMatrix<T>& M , ConstPointer( T2 ) diagonal , ConstPointer( T2 ) b , Pointer( T2 ) x , bool forward , int threads=1 ); 93 template< class T2 > static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseMatrix<T>& M , ConstPointer( T2 ) b , Pointer( T2 ) x , bool forward , int threads=1 ); 94 template< class T2 > static int SolveCG( const SparseMatrix<T>& M , ConstPointer( T2 ) b , int iters , Pointer( T2 ) x , T2 eps=1e-8 , int reset=1 , bool addDCTerm=false , bool solveNormal=false , int threads=1 ); 95 }; 96 97 98 #if !NEW_SPARSE_MATRIX 99 template< class T2 > 100 struct MapReduceVector 101 { 102 private: 103 int _dim; 104 public: 105 std::vector< T2* > out; MapReduceVectorMapReduceVector106 MapReduceVector( void ) { _dim = 0; } ~MapReduceVectorMapReduceVector107 ~MapReduceVector( void ) 108 { 109 if( _dim ) for( int t=0 ; t<int(out.size()) ; t++ ) delete[] out[t]; 110 out.resize( 0 ); 111 } 112 T2* operator[]( int t ) { return out[t]; } 113 const T2* operator[]( int t ) const { return out[t]; } threadsMapReduceVector114 int threads( void ) const { return int( out.size() ); } resizeMapReduceVector115 void resize( int threads , int dim ) 116 { 117 if( threads!=out.size() || _dim<dim ) 118 { 119 for( int t=0 ; t<int(out.size()) ; t++ ) delete[] out[t]; 120 out.resize( threads ); 121 for( int t=0 ; t<int(out.size()) ; t++ ) out[t] = new T2[dim]; 122 _dim = dim; 123 } 124 } 125 126 }; 127 128 template< class T > 129 class SparseSymmetricMatrix : public SparseMatrix< T > 130 { 131 public: 132 133 template< class T2 > 134 Vector< T2 > operator * ( const Vector<T2>& V ) const; 135 136 template< class T2 > 137 Vector< T2 > Multiply( const Vector<T2>& V ) const; 138 139 template< class T2 > 140 void Multiply( const Vector<T2>& In, Vector<T2>& Out , bool addDCTerm=false ) const; 141 142 template< class T2 > 143 void Multiply( const Vector<T2>& In, Vector<T2>& Out , MapReduceVector< T2 >& OutScratch , bool addDCTerm=false ) const; 144 145 template< class T2 > 146 void Multiply( const Vector<T2>& In, Vector<T2>& Out , std::vector< T2* >& OutScratch , const std::vector< int >& bounds ) const; 147 148 template< class T2 > 149 static int SolveCG( const SparseSymmetricMatrix<T>& M , const Vector<T2>& b , int iters , Vector<T2>& x , T2 eps=1e-8 , int reset=1 , int threads=0 , bool addDCTerm=false , bool solveNormal=false ); 150 151 template< class T2 > 152 static int SolveCG( const SparseSymmetricMatrix<T>& M , const Vector<T2>& b , int iters , Vector<T2>& x , MapReduceVector<T2>& scratch , T2 eps=1e-8 , int reset=1 , bool addDCTerm=false , bool solveNormal=false ); 153 #ifdef WIN32 154 template< class T2 > 155 static int SolveCGAtomic( const SparseSymmetricMatrix<T>& M , const Vector<T2>& b , int iters , Vector<T2>& x , T2 eps=1e-8 , int reset=1 , int threads=0 , bool solveNormal=false ); 156 #endif // WIN32 157 template< class T2 > 158 static int SolveJacobi( const SparseSymmetricMatrix<T>& M , const Vector<T2>& diagonal , const Vector<T2>& b , Vector<T2>& x , MapReduceVector<T2>& scratch , Vector<T2>& Mx , T2 sor , int reset ); 159 template< class T2 > 160 static int SolveJacobi( const SparseSymmetricMatrix<T>& M , const Vector<T2>& b , int iters , Vector<T2>& x , MapReduceVector<T2>& scratch , T2 sor=T2(1.) , int reset=1 ); 161 template< class T2 > 162 static int SolveJacobi( const SparseSymmetricMatrix<T>& M , const Vector<T2>& diagonal , const Vector<T2>& b , Vector<T2>& x , Vector<T2>& Mx , T2 sor , int reset ); 163 template< class T2 > 164 static int SolveJacobi( const SparseSymmetricMatrix<T>& M , const Vector<T2>& b , int iters , Vector<T2>& x , T2 sor=T2(1.) , int reset=1 ); 165 166 enum 167 { 168 ORDERING_UPPER_TRIANGULAR , 169 ORDERING_LOWER_TRIANGULAR , 170 ORDERING_NONE 171 }; 172 template< class T2 > 173 static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseSymmetricMatrix<T>& M , const Vector<T2>& diagonal , const Vector<T2>& b , Vector<T2>& x , MapReduceVector<T2>& scratch , Vector<T2>& Mx , Vector<T2>& dx , bool forward , int reset ); 174 template< class T2 > 175 static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseSymmetricMatrix<T>& M , const Vector<T2>& b , int iters , Vector<T2>& x , MapReduceVector<T2>& scratch , bool forward , int reset=1 ); 176 177 template< class T2 > 178 static int SolveGS( const SparseSymmetricMatrix<T>& M , const Vector<T2>& diagonal , const Vector<T2>& b , Vector<T2>& x , MapReduceVector<T2>& scratch , Vector<T2>& Mx , Vector<T2>& dx , bool forward , int reset , int ordering ); 179 template< class T2 > 180 static int SolveGS( const SparseSymmetricMatrix<T>& M , const Vector<T2>& b , int iters , Vector<T2>& x , MapReduceVector<T2>& scratch , bool forward , int reset=1 , int ordering=ORDERING_NONE ); 181 template< class T2 > 182 static int SolveGS( const SparseSymmetricMatrix<T>& M , const Vector<T2>& diagonal , const Vector<T2>& b , Vector<T2>& x , Vector<T2>& Mx , Vector<T2>& dx , bool forward , int reset , int ordering ); 183 template< class T2 > 184 static int SolveGS( const SparseSymmetricMatrix<T>& M , const Vector<T2>& b , int iters , Vector<T2>& x , bool forward , int reset=1 , int ordering=ORDERING_NONE ); 185 186 template< class T2 > 187 void getDiagonal( Vector< T2 >& diagonal , int threads=1 ) const; 188 }; 189 #endif // !NEW_SPARSE_MATRIX 190 191 #include "SparseMatrix.inl" 192 193 #endif 194 195