1 #ifndef SPARSEMATRIX_H 2 #define SPARSEMATRIX_H 3 #include <exception> 4 #include "MatrixLibrary.h" 5 #include <algorithm> 6 7 namespace ATC_matrix { 8 9 /** 10 * @struct TRI_COORD 11 * @brief Triplet SparseMatrix entry 12 */ 13 template <typename T> 14 struct TRI_COORD 15 { 16 TRI_COORD<T>(INDEX row=0, INDEX col=0); 17 TRI_COORD<T>(INDEX row, INDEX col, T val, bool add_to=0); 18 INDEX i, j; 19 T v; 20 bool add; 21 }; 22 23 template<typename T> 24 void ParMultAB(MPI_Comm comm, const SparseMatrix<T>& A, const Matrix<T>& B, DenseMatrix<T>& C); 25 26 /** 27 * @class SparseMatrix 28 * @brief Stores data in triplet format or CRS format 29 */ 30 template<typename T> 31 class SparseMatrix : public Matrix<T> 32 { 33 //* SparseMatrix-Vector multiplication (S * v) 34 friend DenseVector<T> operator*<T>(const SparseMatrix<T> &A, const Vector<T>& x); 35 //* SparseMatrix-DenseMatrix multiplication (S * F) 36 friend DenseMatrix<T> operator*<T>(const SparseMatrix<T> &A, const Matrix<T>& D); 37 //* SparseMatrix-DiagonalMatrix multiplication (S * D) 38 friend SparseMatrix<T> operator*<T>(const SparseMatrix<T> &A, const DiagonalMatrix<T>& D); 39 //* SparseMatrix-SparseMatrix multiplication (S * S) 40 friend SparseMatrix<T> operator*<T>(const SparseMatrix<T> &A, const SparseMatrix<T> &B); 41 //* computes the product of a SparseMatrix tranpose with a SparseVector (M'*v). 42 friend SparseVector<T> operator*<T>(const SparseMatrix<T> &M, const SparseVector<T> &v); 43 //* computes the product of a SparseMatrix tranpose with a SparseVector (M'*v). 44 friend SparseVector<T> operator*<T>(const SparseVector<T> &v, const SparseMatrix<T> &M); 45 46 template<typename U> 47 friend void ParMultAB(MPI_Comm comm, const SparseMatrix<U>& A, const Matrix<U>& B, DenseMatrix<U>& C); 48 49 public: 50 SparseMatrix(INDEX rows=0, INDEX cols=0); 51 SparseMatrix(const SparseMatrix<T>& c); 52 SparseMatrix(const DenseMatrix<T>& c); 53 SparseMatrix(INDEX* rows, INDEX* cols, T* vals, INDEX size, 54 INDEX nRows, INDEX nCols, INDEX nRowsCRS); ~SparseMatrix()55 virtual ~SparseMatrix() { _delete(); } 56 57 //* General index by value (requires a binary search on the row) 58 T operator()(INDEX i, INDEX j) const; 59 //* General index by reference (requires a binary search on the row) 60 T& operator()(INDEX i, INDEX j); 61 //* General flat index by value operator (by nth nonzero) 62 T operator[](INDEX i) const; 63 //* General flat index by reference operator (by nth nonzero) 64 T& operator[](INDEX i); 65 //* sets a value to index i,j 66 void set(INDEX i, INDEX j, T v); 67 //* adds a value to index i,j 68 void add(INDEX i, INDEX j, T v); 69 //* return a triplet value of the ith nonzero 70 TRIPLET<T> triplet(INDEX i) const; 71 72 //* full reset - completely wipes out all SparseMatrix data 73 void reset(INDEX rows=0, INDEX cols=0, bool zero=true); 74 //* only changes the bounds of the matrix, no deletion 75 void resize(INDEX rows=0, INDEX cols=0, bool zero=true); 76 //* reset - from DenseMatrix - this will be SLOW 77 void reset(const DenseMatrix<T>& D, double TOL=-1.0); 78 //* copy data 79 void copy(const T * ptr, INDEX rows=0, INDEX cols=0); 80 81 void dense_copy(DenseMatrix<T>& D) const; 82 DenseMatrix<T> dense_copy(void) const; 83 84 //* returns true if the matrix has no nonzero elements 85 bool empty() const; 86 //* returns the user-specified number of rows 87 INDEX nRows() const; 88 INDEX nRowsCRS() const; 89 //* returns the user-specified number of cols 90 INDEX nCols() const; 91 //* returns the number of non-zero elements 92 INDEX size() const; 93 //* returns the number of non-zeros in a row 94 INDEX RowSize(INDEX r) const; 95 //* returns a pointer to the CRS list of rows 96 inline INDEX* rows() const; 97 //* returns a pointer to the CRS list of cols 98 inline INDEX* cols() const; 99 //* returns a pointer to the nonzero data 100 inline T* ptr() const; 101 //* checks if the index i,j falls in the user-specified range 102 bool in_range(INDEX i, INDEX j) const; 103 //* check if the total matrix has a value set for an index pair 104 bool has_entry(INDEX i, INDEX j) const; 105 //* check if the uncompressed part of the matrix has a value set for an index pair 106 bool has_entry_uncompressed(INDEX i, INDEX j) const; 107 //* check if the compressed part matrix has a value set for an index pair 108 bool has_entry_compressed(INDEX i, INDEX j) const; 109 //* check if the matrix has been compressed at least once 110 bool has_template(void) const; 111 112 /* 113 * \section assignment operators 114 */ 115 //* copies SparseMatrix R to this 116 SparseMatrix<T>& operator=(const SparseMatrix &R); 117 //* sets all nonzero values to a constant 118 SparseMatrix<T>& operator=(const T v); 119 //* scales all nonzero values by a constant 120 SparseMatrix<T>& operator*=(const T &a); 121 //* calls operator*= of base class 122 SparseMatrix<T>& operator*=(const SparseMatrix<T> &a); 123 // Adds two matrices together. 124 SparseMatrix<T>& operator+=(const SparseMatrix & R); 125 126 /* 127 * \section Multiplication operations 128 */ 129 130 //----------------------------------------------------------------------------- 131 // multiply sparse matrix by a vector 132 //----------------------------------------------------------------------------- MultMv(const Vector<T> & v,DenseVector<T> & c)133 virtual void MultMv(const Vector<T>& v, DenseVector<T>& c) const 134 { 135 compress(*this); 136 GCK(*this, v, this->nCols() != v.size(), "SparseMatrix * Vector") 137 138 // resize c if necessary 139 if (c.size() != this->nRows()) { 140 c.resize(this->nRows()); 141 c.zero(); 142 } 143 144 INDEX i, j; 145 for (i = 0; i < this->_nRowsCRS; i++) 146 for (j = this->_ia[i]; j < this->_ia[i + 1]; j++) 147 c(i) += this->_val[j] * v(this->_ja[j]); 148 } 149 150 //----------------------------------------------------------------------------- 151 // multiply sparse matrix by dense matrix 152 //----------------------------------------------------------------------------- MultAB(const Matrix<T> & B,DenseMatrix<T> & C)153 virtual void MultAB(const Matrix<T>& B, DenseMatrix<T>& C) const 154 { 155 GCK(*this, B, this->nCols() != B.nRows(), "SparseMatrix * DenseMatrix") 156 157 const INDEX J = B.nCols(); 158 159 INDEX i, ik, j; 160 for (i = 0; i < this->_nRowsCRS; i++) 161 for (ik = this->_ia[i]; ik < this->_ia[i + 1]; ik++) 162 for (j = 0; j < J; j++) 163 C(i, j) += this->_val[ik] * B(this->_ja[ik], j); // C(i,j) = S(i,k) * B(k, j) 164 } 165 166 //----------------------------------------------------------------------------- 167 // Multiplies this SparseMatrix transposed times a vector 168 //----------------------------------------------------------------------------- transMat(const Vector<T> & x)169 virtual DenseVector<T> transMat(const Vector<T> &x) const 170 { 171 compress(*this); 172 DenseVector<T> y(nCols(), true); 173 GCK(*this, x, nRows()!=x.size(),"operator *: Sparse matrix incompatible with Vector.") 174 175 INDEX i, ij; 176 for(i=0; i<_nRowsCRS; i++) 177 for(ij=_ia[i]; ij<_ia[i+1]; ij++) 178 y(_ja[ij]) += _val[ij]*x(i); 179 return y; 180 } 181 182 //----------------------------------------------------------------------------- 183 // Matrix Transpose/DenseMatrix multiply 184 //----------------------------------------------------------------------------- transMat(const DenseMatrix<T> & D)185 virtual DenseMatrix<T> transMat(const DenseMatrix<T> &D) const 186 { 187 compress(*this); 188 GCK(*this, D, nRows()!=D.nRows(),"transMat: Sparse matrix incompatible with DenseMatrix.") 189 DenseMatrix<T> C(nCols(), D.nCols(), true); // initialized to zero 190 INDEX j, k, ki; 191 192 for (k=0; k<_nRowsCRS; k++) 193 for (ki=_ia[k]; ki<_ia[k+1]; ki++) 194 for (j=0; j<D.nCols(); j++) 195 C(_ja[ki], j) += _val[ki]*D(k,j); // C(i,j) = S(k,i) * D(k, j) 196 197 return C; 198 } 199 200 //----------------------------------------------------------------------------- 201 // Matrix Transpose/SparseMatrix multiply - IS THIS REALLY NEEDED?? 202 //----------------------------------------------------------------------------- transMat(const SparseMatrix<T> & D)203 virtual DenseMatrix<T> transMat(const SparseMatrix<T> &D) const 204 { 205 compress(*this); 206 GCK(*this, D, nRows()!=D.nRows(),"transMat: Sparse matrix incompatible with DenseMatrix.") 207 DenseMatrix<T> C(nCols(), D.nCols(), true); // initialized to zero 208 209 INDEX k, ki, kj; 210 for (k=0; k<_nRowsCRS; k++) 211 for (kj=D._ia[k]; kj<D._ia[k+1]; kj++) 212 for (ki=_ia[k]; ki<_ia[k+1]; ki++) 213 C(_ja[ki], D._ja[kj]) += _val[ki]*D._val[kj]; // C(i,j) = S(k,i)*D(k,j) 214 215 return C; 216 } 217 218 SparseMatrix<T> transpose() const; 219 SparseMatrix<T>& row_scale(const Vector<T> &v); 220 SparseMatrix<T>& col_scale(const Vector<T> &v); 221 DenseVector<T> col_sum() const; 222 DenseVector<INDEX> column_count() const; 223 DiagonalMatrix<T> diag() const; 224 DiagonalMatrix<T> row_sum_lump() const; 225 void row(INDEX i, DenseVector<T>& row, DenseVector<INDEX>& indx) const; 226 void weighted_least_squares(const SparseMatrix<T> &N, const DiagonalMatrix<T> &D); 227 void set_all_elements_to(const T &v); 228 229 T row_max(INDEX row) const; 230 T row_min(INDEX row) const; 231 232 /* 233 * \section I/O functions 234 */ 235 //* outputs this SparseMatrix to a formatted string 236 std::string to_string() const; 237 using Matrix<T>::matlab; 238 //* writes a command to recreate this matrix in matlab to a stream 239 void matlab(std::ostream &o, const std::string &name="S") const; 240 //* prints a row histogram for each row 241 void print_row_histogram(const std::string &name, INDEX nbins = 10) const; 242 //* prints a histogram of the values in a row 243 void print_row_histogram(INDEX row, INDEX nbins) const; 244 //* prints the current triplets 245 void print_triplets() const; 246 //! Writes the matrix to a binary file (after a compress). 247 void binary_write(std::fstream& f) const; 248 //! Reads a SparseMatrix from a binary file. (wipes out any original data) 249 void binary_read(std::fstream& f); 250 //* Dump templated type to disk; operation not safe for all types 251 void write_restart(FILE *f) const; 252 253 /* 254 * \section Utility functions 255 */ 256 //* converts all triplets and merges with CRS 257 void compress(); 258 //* converts T to CRS 259 static void compress(const SparseMatrix<T> &C); 260 //* sorts and returns the # of unique triplets 261 INDEX CountUniqueTriplets(); 262 263 private: 264 //* creates a CRS structure 265 void _create(INDEX size, INDEX nrows); 266 //* clears all memory and nulls references 267 void _delete(); 268 //* copies all data from another SparseMatrix 269 void _copy(const SparseMatrix<T> &C); 270 //* general sparse matrix assignment 271 void _set_equal(const Matrix<T> &r); 272 //* returns the first row with a nonzero in it (from the CRS structure only) 273 int _first_nonzero_row_crs() const; 274 275 /* 276 * \section CRS storage variables 277 */ 278 protected: 279 T * _val; // matrix non-zeros 280 INDEX *_ia, *_ja; // ptrs to rows, column indexes 281 INDEX _size, _nRowsCRS; // # of non-zeros, rows 282 bool hasTemplate_; 283 284 void copy(const SparseMatrix<T> &C); 285 286 //* new (unsorted triplet values - won't intersect CRS values) 287 mutable std::vector<TRI_COORD<T> > _tri; 288 /* 289 * \section User specified variables 290 */ 291 INDEX _nRows, _nCols; 292 static T _zero; 293 }; 294 295 } // end namespace 296 297 #include "SparseMatrix-inl.h" 298 #endif 299 300