1 /* 2 HMat-OSS (HMatrix library, open source software) 3 4 Copyright (C) 2014-2015 Airbus Group SAS 5 6 This program is free software; you can redistribute it and/or 7 modify it under the terms of the GNU General Public License 8 as published by the Free Software Foundation; either version 2 9 of the License, or (at your option) any later version. 10 11 This program is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 GNU General Public License for more details. 15 16 You should have received a copy of the GNU General Public License 17 along with this program; if not, write to the Free Software 18 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 19 20 http://github.com/jeromerobert/hmat-oss 21 */ 22 23 #ifndef _RK_MATRIX_HPP 24 #define _RK_MATRIX_HPP 25 /* Implementation of Rk-matrices */ 26 #include <vector> 27 #include <algorithm> 28 #include <utility> 29 30 #include "full_matrix.hpp" 31 #include "compression.hpp" 32 #include "common/my_assert.h" 33 34 namespace hmat { 35 36 template<typename T> class HMatrix; 37 template<typename T, template <typename> class F> class AssemblyFunction; 38 class ClusterData; 39 class IndexSet; 40 41 /** Control the approximation of Rk-matrices. 42 43 In the case where k != 0, we do an approximation with a fixed rank k, 44 otherwise the approximation is adaptive, it stops when the 45 singular value is less than an error relative to the 46 sum of the singular values in the SVD. 47 */ 48 class RkApproximationControl { 49 public: 50 double coarseningEpsilon; /// Tolerance for the coarsening 51 int compressionMinLeafSize; 52 53 /** Initialization with impossible values by default 54 */ RkApproximationControl()55 RkApproximationControl() : coarseningEpsilon(-1.), 56 compressionMinLeafSize(100) {} 57 }; 58 59 60 template<typename T> class RkMatrix { 61 // Type double precision associated to T 62 typedef typename Types<T>::dp dp_t; 63 /** Swaps members of two RkMatrix instances. 64 Since rows and cols are constant, they cannot be swaped and 65 the other instance must have the same members. 66 67 \param other other RkMatrix instance. 68 */ 69 void swap(RkMatrix<T>& other); 70 /** Recompress an RkMatrix in place with a modified Gram-Schmidt algorithm. 71 72 @warning The previous rk->a and rk->b are no longer valid after this function. 73 \param epsilon is the accuracy of the recompression 74 \param initialPivotA/B is the number of orthogonal columns in panels a and b 75 */ 76 void mGSTruncate(double epsilon, int initialPivotA=0, int initialPivotB=0); 77 public: 78 /** @brief A hook which can be called at the begining of formatedAddParts */ 79 static bool (*formatedAddPartsHook)(RkMatrix<T> * me, double epsilon, const T* alpha, const RkMatrix<T>* const * parts, const int n); 80 const IndexSet *rows; 81 const IndexSet *cols; 82 // A B^t 83 ScalarArray<T>* a; 84 ScalarArray<T>* b; 85 86 /// Control of the approximation. See \a RkApproximationControl for more 87 /// details. 88 static RkApproximationControl approx; 89 90 /** Construction of a RkMatrix . 91 92 A Rk-matrix is a compressed representation of a matrix of size 93 n*m by a matrix R = AB^t with A a n*k matrix and B is a m*k matrix. 94 The represented matrix R corresponds to a set of indices _rows and _cols. 95 96 \param _a matrix A 97 \param _rows indices of the rows (of size n) 98 \param _b matrix B (not B^t) 99 \param _cols indices of the columns (of size k) 100 */ 101 RkMatrix(ScalarArray<T>* _a, const IndexSet* _rows, 102 ScalarArray<T>* _b, const IndexSet* _cols); 103 ~RkMatrix(); 104 rank() const105 int rank() const { 106 return a ? a->cols : 0; 107 } 108 109 /** Returns a pointer to a new RkMatrix representing a subset of indices. 110 The pointer is supposed to be read-only (for efficiency reasons). 111 112 \param subRows subset of rows 113 \param subCols subset of columns 114 \return pointer to a new matrix with subRows and subCols. 115 */ 116 const RkMatrix* subset(const IndexSet* subRows, const IndexSet* subCols) const; 117 RkMatrix* truncatedSubset(const IndexSet* subRows, const IndexSet* subCols, double epsilon) const; 118 /** Returns the compression ratio (stored_elements, total_elements). 119 */ 120 size_t compressedSize(); 121 122 size_t uncompressedSize(); 123 124 /** Returns a pointer to a new FullMatrix M = AB^t (uncompressed) 125 */ 126 FullMatrix<T>* eval() const; 127 /** Returns a pointer to a new ScalarArray M = AB^t (uncompressed) or fill an existing one 128 */ 129 ScalarArray<T>* evalArray(ScalarArray<T> *result=NULL) const ; 130 /** Recompress an RkMatrix in place. 131 132 @warning The previous rk->a and rk->b are no longer valid after this function. 133 \param initialPivotA/B is the number of orthogonal columns in panels a and b 134 */ 135 void truncate(double epsilon, int initialPivotA=0, int initialPivotB=0); 136 /** Add randomness to the RkMatrix */ 137 void addRand(double epsilon); 138 /*! \brief Return square of the Frobenius norm of the matrix. 139 140 \return the matrix norm. 141 */ 142 double normSqr() const; 143 /** this <- this + alpha * mat 144 145 \param alpha 146 \param mat 147 */ 148 void axpy(double epsilon, T alpha, const FullMatrix<T>* mat); 149 /** this <- this + alpha * mat 150 151 \param alpha 152 \param mat 153 */ 154 void axpy(double epsilon, T alpha, const RkMatrix<T>* mat); 155 /** Adds a list of RkMatrix to a RkMatrix. 156 157 In this function, RkMatrix may include some 158 only indices, that is to say be subsets of 159 evidence of this. 160 161 \param units The list RkMatrix adding. 162 \param n Number of matrices to add 163 \param epsilon truncate epsilon (negative value disable truncate) 164 \param hook true to call formattedAddPartsHook, false to do not call 165 \return truncate(*this + parts[0] + parts[1] + ... + parts[n-1]) 166 */ 167 void formattedAddParts(double epsilon, const T* alpha, const RkMatrix<T>* const * parts, const int n, 168 bool hook = true); 169 /** Adds a list of MatrixXd (solid matrices) to RkMatrix. 170 171 In this function, MatrixXd may cover a portion of 172 index only, that is to say be subsets of indices 173 this. 174 The MatrixXd are converted RkMatrix before the addition, which is 175 with RkMatrix :: formattedAddParts. 176 177 \param units The list MatrixXd adding. 178 \param n Number of matrices to add 179 \return truncate (*this + parts[0] + parts[1] + ... + parts[n-1]) 180 */ 181 void formattedAddParts(double epsilon, const T* alpha, const FullMatrix<T>* const * parts, int n); 182 183 /*! \brief Add a product of HMatrix to an RkMatrix 184 */ 185 void gemmRk(double epsilon, char transA, char transB, T alpha, const HMatrix<T>* a, const HMatrix<T>* b); 186 187 /** Multiplication by a scalar. 188 189 \param alpha the scalar 190 */ 191 void scale(T alpha); 192 /** \brief Transpose in place. 193 */ 194 void transpose(); 195 void clear(); 196 197 /** Copy RkMatrix into this. 198 */ 199 void copy(const RkMatrix<T>* o); 200 201 /** Return a copy of this. 202 */ 203 RkMatrix<T>* copy() const; 204 205 /** Compute y <- alpha * op(A) * y + beta * y if side is Side::LEFT 206 or y <- alpha * y * op(A) + beta * y if side is Side::RIGHT 207 with x and y ScalarArray<T>*. 208 209 The arguments are similar to BLAS GEMV. 210 */ 211 void gemv(char trans, T alpha, const ScalarArray<T>* x, T beta, ScalarArray<T>* y, Side side = Side::LEFT) const; 212 213 /** Right multiplication of RkMatrix by a matrix. 214 215 \param transR 'N' or 'T' depending on whether R is transposed or not 216 \param Beam 'N' or 'T' depending on whether M is transposed or not 217 \return R * M 218 */ 219 static RkMatrix<T>* multiplyRkFull(char transR, char transM, 220 const RkMatrix<T>* rk, const FullMatrix<T>* m); 221 222 /** Left multiplication of RkMatrix by a matrix. 223 224 \param transR 'N' or 'T' depending on whether R is transposed or not 225 \param Beam 'N' or 'T' depending on whether M is transposed or not 226 \return R * M 227 228 */ 229 static RkMatrix<T>* multiplyFullRk(char transM, char transR, 230 const FullMatrix<T>* m, 231 const RkMatrix<T>* rk); 232 /* These functions are added to manage the particular case of the product by 233 H-matrix, which is treated by decomposing the product into the succession of 234 products by a vector, the result being a RkMatrix. 235 */ 236 /** Right multiplying of RkMatrix by HMatrix. 237 238 The product is made by multiplying the vector by vector. 239 240 \param R rk 241 \param h H 242 \param transRk 'N' or 'T' depending on whether or not R is transposed 243 \param transh 'N' or 'T' depending on whether or not H is transposed 244 \return R * H 245 */ 246 static RkMatrix<T>* multiplyRkH(char transRk, char transH, const RkMatrix<T>* rk, const HMatrix<T>* h); 247 /** Left multiplying a RkMatrix by HMatrix. 248 249 The product is made by multiplying vector by vector. 250 251 \param h H 252 \param R rk 253 \param transR 'N' or 'T' depending on whether or not R is transposed 254 \param transh 'N' or 'T' depending on whether or not H is transposed 255 \return H * R 256 */ 257 static RkMatrix<T>* multiplyHRk(char transH, char transR, const HMatrix<T>* h, const RkMatrix* rk); 258 /** Multiplying a RkMatrix by a RkMatrix 259 260 The product is made by multiplying vector by vector. 261 262 \param a 263 \param b 264 \return A * B 265 */ 266 static RkMatrix<T>* multiplyRkRk(char transA, char transB, const RkMatrix<T>* a, const RkMatrix<T>* b, double epsilon); 267 /*! \brief in situ multiplication of the matrix by the diagonal of the matrix given as argument 268 269 \param d D matrix which we just considered the diagonal 270 \param inverse D or D^{-1} 271 \param left multiplication (side = Side::LEFT) or right (side = Side::RIGHT) of this by D 272 */ 273 void multiplyWithDiagOrDiagInv(const HMatrix<T> * d, bool inverse, Side side = Side::RIGHT); 274 /*! \brief Triggers an assertion if there are NaNs in the RkMatrix. 275 */ 276 void checkNan() const; 277 278 /*! Conjugate the content of the complex matrix */ 279 void conjugate(); 280 281 /** Return the memory size of the a*b product */ 282 static size_t computeRkRkMemorySize(char transA, char transB, const RkMatrix<T>* a, const RkMatrix<T>* b); 283 284 /** Simpler accessor for the data. 285 286 There is only 1 type to because we cant modify the data. 287 */ 288 T get(int i, int j) const ; 289 290 /*! \brief Write the RkMatrix data 'a' and 'b' in a stream (FILE*, unix fd, ...) 291 */ 292 void writeArray(hmat_iostream writeFunc, void * userData) const; 293 }; 294 295 } // end namespace hmat 296 297 #endif 298