Mercurial > hg > kwantix
annotate include/linalg.hpp @ 12:6e06eb6ec448
Trivial change in include guards.
| author | Jordi Gutiérrez Hermoso <jordigh@gmail.com> |
|---|---|
| date | Mon, 07 Jul 2008 13:07:11 -0500 |
| parents | d0076d9b2ef1 |
| children | 0531194a431d |
| rev | line source |
|---|---|
| 0 | 1 /*!\file linalg.hpp |
| 2 * \brief Wrapper linear algebra classes for the GSL. | |
| 3 * | |
| 4 * This header file puts some C++ wrappers around the GSL vector and | |
| 5 * matrix structures plus providing a C++ interface to some pertinent | |
| 6 * BLAS and general linear algebra routines. | |
| 7 */ | |
| 8 | |
|
12
6e06eb6ec448
Trivial change in include guards.
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
11
diff
changeset
|
9 #ifndef __LINALG_HPP__ |
|
6e06eb6ec448
Trivial change in include guards.
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
11
diff
changeset
|
10 #define __LINALG_HPP__ |
| 0 | 11 |
| 12 #include <gsl/gsl_matrix.h> | |
| 13 #include <gsl/gsl_vector.h> | |
| 14 #include <gsl/gsl_permutation.h> | |
| 15 #include "error.hpp" | |
| 16 | |
| 17 /// Linear algebra namespace | |
| 18 namespace linalg{ | |
| 19 using namespace error_handling; | |
| 20 | |
| 21 class slice; | |
| 22 class vector; | |
| 23 class vector_view; | |
| 24 | |
| 25 /// A wrapper class for GNU Scientific Library matrices | |
| 26 class matrix{ | |
| 27 private: | |
| 28 class LUmatrix; | |
| 29 | |
| 30 public: | |
| 31 //Constructors | |
| 32 | |
| 33 /// Allocate 1x1 matrix zero matrix. | |
| 34 matrix(); | |
| 35 /*! \brief Allocate m by n matrix, filled with some value. | |
| 36 \param m - Number of rows. | |
| 37 \param n - Number of columns. | |
| 38 \param fillvalue - Value to fill matrix with, default is 0. | |
| 39 */ | |
| 40 matrix(const size_t m, const size_t n, const double fillvalue = 0) ; | |
| 41 /// Assign a matrix from a GSL matrix pointer. | |
| 42 matrix(gsl_matrix *M); | |
| 43 | |
| 44 /// Copy constructor | |
| 45 matrix(const matrix& M); | |
| 46 /// Assignment. | |
| 47 matrix operator=(const matrix& M); | |
| 48 | |
| 49 //Destructors | |
| 50 | |
| 51 /// Clears the memory allocated by the GSL | |
| 52 ~matrix(); | |
| 53 | |
| 54 /// Number of decimal digits to output. | |
| 55 size_t precision() const; | |
| 56 /// Set the number of decimal digits to output. | |
| 57 static void set_precision(size_t p); | |
| 58 | |
| 59 /// Number of rows. | |
| 60 size_t rows() const; | |
| 61 /// Number of columns. | |
| 62 size_t cols() const; | |
| 63 | |
| 64 /*! \brief Fortran-style parenthetical indexing (hence Octave-style too) | |
| 65 * Indices start from 1. | |
| 66 * @param i - Row number. | |
| 67 * @param j - Column number. | |
| 68 * @return A reference to the matrix element. | |
| 69 */ | |
| 70 double& operator()(const size_t i, const size_t j); | |
| 71 | |
| 72 /// Constant version of previous function. | |
| 73 const double& operator()(const size_t i, const size_t j) const; | |
| 74 | |
| 75 /*! \brief Indexing for vectorisation. | |
| 76 * The GSL provides limited vectorisation routines which can be | |
| 77 * useful for operating on blocks of matrices all at once. | |
| 78 * @param i - Row to be sliced. | |
| 79 * @param b - Slice range. | |
| 80 * @return A vector_view of the sliced row. | |
| 81 */ | |
| 82 vector_view operator()(const size_t i, const slice &b); | |
| 83 /// Constant version of previous function. | |
| 84 const vector_view operator()(const size_t i, const slice &b) const; | |
| 85 | |
| 86 /*! \brief Indexing for vectorisation. | |
| 87 * The GSL provides limited vectorisation routines which can be | |
| 88 * useful for operating on blocks of matrices all at once. | |
| 89 * @param a - Slice range | |
| 90 * @param j - column to be sliced. | |
| 91 * @return A vector_view of the sliced column. | |
| 92 */ | |
| 93 vector_view operator()(const slice &a, const size_t j); | |
| 94 /// Constant version of previous function. | |
| 95 const vector_view operator()(const slice &a, const size_t j) const; | |
| 96 | |
| 97 //Arithmetic operations | |
| 98 ///Scale the matrix. | |
| 99 matrix operator*(const double a) const; | |
| 100 /// Matrix addition | |
| 101 matrix operator+(const matrix& N) const; | |
| 102 /// Matrix multiplication | |
| 103 matrix operator*(const matrix& N) const; | |
| 104 /// Matrix subtraction | |
| 105 matrix operator-(const matrix& N) const; | |
| 106 /// Matrix vector product (gaxpy operation). | |
| 107 vector operator*(const vector& v) const; //Mv, where v is treated | |
| 108 //as a column vector. | |
| 109 | |
| 110 //More complex operations. | |
| 111 | |
| 112 /*! \brief Inverse. | |
| 113 * | |
| 114 * This function returns the inverse matrix, computed via an LU | |
| 115 * factorisation. If the matrix is ill-conditioned or the inverse | |
| 116 * does not exist, it will happily return a matrix full of NaN. | |
| 117 * \returns The inverse matrix. | |
| 118 */ | |
| 119 matrix inv() const; | |
| 120 /*! \brief Tranpose. | |
| 121 * Returns a copy of the original matrix, transposed. The original | |
| 122 * matrix is not modified. | |
| 123 * \returns Transposed copy of matrix. | |
| 124 */ | |
| 125 matrix T() const; | |
| 126 /*! \brief Trace. | |
| 127 * The trace of a square matrix. | |
| 128 * \returns The trace. | |
| 129 */ | |
| 130 double tr() const; | |
| 131 /*! \brief Determinant. | |
| 132 * The determinant computed via an LU factorisation | |
| 133 * \returns The determinant. | |
| 134 */ | |
| 135 double det() const; | |
| 136 | |
| 137 /*! \brief Solves Mv = w for v with LU factorisation. | |
| 138 * Given another vector, will return the solution vector v to the | |
| 139 * equation Mv = w. | |
| 140 * \param w - Another vector, must have the same size as number of | |
| 141 * rows in the matrix. | |
| 142 * \returns The solution v to the equation Mv = w. | |
| 143 */ | |
| 144 vector inv(const vector& w) const; | |
| 145 | |
| 146 /*! \brief L2 condition number, using SVD. | |
| 147 * | |
| 148 * This function computes the L2 condition number of a matrix, the | |
| 149 * largest singular value divided by the smallest. It directly | |
| 150 * uses the SVD decomposition (as provided by the GSL), and it may | |
| 151 * take a long time for large matrices. | |
| 152 * \returns The L2 condition number. | |
| 153 */ | |
| 154 double cond() const; | |
| 155 | |
| 156 friend class vector_view; | |
| 157 private: | |
| 158 /*! \brief LU decomposition, pivots in U. | |
| 159 * | |
| 160 * This function will compute the LU factorisation of a matrix, | |
| 161 * where the pivots are in the upper triangular matrix U and L has | |
| 162 * ones along the diagonal (not stored). For efficiency, L and U | |
| 163 * are both internally stored within a single matrix. | |
| 164 * \returns A pointer to the private class LUmatrix. | |
| 165 */ | |
| 166 LUmatrix* LU() const; | |
| 167 | |
| 168 gsl_matrix * A; | |
| 169 | |
| 170 ///Precision | |
| 171 static size_t precsn; | |
| 172 /// Machine epsilon | |
| 173 static const double eps; | |
| 174 | |
| 175 // Internal members, all can change, none of them affect the | |
| 176 // actual matrix per se. | |
| 177 mutable LUmatrix* LUptr; | |
| 178 mutable bool LUfactored; | |
| 179 | |
| 180 mutable double condition_number; //L2 condition number, obtained by svd. | |
| 181 mutable gsl_vector* SVD; //Matrix's singular values. | |
| 182 mutable bool SVDfactored; | |
| 183 | |
| 184 void SVDfactor() const; | |
| 185 | |
| 186 /// A private matrix member class for matrices factored in LU form. | |
| 187 class LUmatrix{ | |
| 188 public: | |
| 189 gsl_matrix* matrix_ptr(); | |
| 190 gsl_permutation* perm_ptr(); | |
| 191 int sgn(); | |
| 192 int* sgn_ptr(); | |
| 193 | |
| 194 LUmatrix(gsl_matrix* M); | |
| 195 LUmatrix(const LUmatrix& LU); | |
| 196 LUmatrix operator=(const LUmatrix& LU); | |
| 197 ~LUmatrix(); | |
| 198 private: | |
| 199 // Private interface to the GSL. | |
| 200 gsl_permutation* p; | |
| 201 gsl_matrix* A; | |
| 202 int signum; | |
| 203 LUmatrix(); | |
| 204 };//End LUmatrix | |
| 205 }; //End matrix | |
| 206 | |
| 207 /// A wrapper class for GSL vectors | |
| 208 class vector{ | |
| 209 public: | |
| 210 //Constructors | |
| 211 | |
| 212 ///Allocate zero vector of size one. | |
| 213 vector(); | |
| 214 /*! \brief Allocate GSL vector of size n, filled with fillvalue | |
| 215 * \param n - Size of new vector | |
| 216 * \param fillvalue - Value for filling vector, default is 0. | |
| 217 */ | |
| 218 vector(const size_t n, const double fillvalue = 0); | |
| 219 /// Associate a vector to a GSL vector pointer. | |
| 220 vector(gsl_vector *y); | |
| 221 /// Associate a vector to a constant GSL vector pointer. | |
| 222 vector(const gsl_vector *y); | |
| 223 | |
| 224 /// Copy contstructor. | |
| 225 vector(const vector& y); | |
| 226 /// Assignment. | |
| 227 vector& operator=(const vector &y); | |
| 228 | |
| 229 //Destructor | |
| 230 /// Clear memory assigned to vector by GSL. | |
| 231 ~vector(); | |
| 232 | |
| 233 //Number of decimal digits to output. | |
| 234 /// Number of decimal digits to output | |
| 235 size_t precision() const; | |
| 236 /// Set the number of decimal digits to output | |
| 237 static void set_precision(size_t p); | |
| 238 | |
| 239 ///Number of elements | |
| 240 size_t size() const; | |
| 241 | |
| 242 /*! \brief Fortran-style parenthetical indexing (hence | |
| 243 * Octave-style too). Indices start at 1. | |
| 244 * \param i - Index number. | |
| 245 * \return A reference to the vector element. | |
| 246 */ | |
| 247 double& operator()(const size_t i); | |
| 248 | |
| 249 /// Constant version of previous function. | |
| 250 const double& operator()(const size_t i) const; | |
| 251 | |
| 252 /*! \brief Indexing for vectorisation. | |
| 253 * The GSL provides limited vectorisation routines which can be | |
| 254 * useful for operating on blocks of vectors all at once. | |
| 255 * @param v - Slice range. | |
| 256 * @return A vector_view of the slice. | |
| 257 */ | |
| 258 vector_view operator()(const slice &v); | |
| 259 ///Constant version of previous function. | |
| 260 const vector_view operator()(const slice &v) const; | |
| 261 | |
| 262 //Arithmetic operations: | |
| 263 ///Scale the vector. | |
| 264 vector operator*(const double a) const; | |
| 265 ///Vector addition. | |
| 266 vector operator+(const vector& w) const; | |
| 267 ///Vector subtration. | |
| 268 vector operator-(const vector& w) const; | |
|
11
d0076d9b2ef1
Conclude precomputation revision
Jordi Guitérrez Hermoso <jordigh@gmail.com>
parents:
0
diff
changeset
|
269 ///Dot product. |
| 0 | 270 double operator*(const vector& w) const; |
| 271 ///Computes vM where v is treated as a row vector. | |
| 272 vector operator*(const matrix& M) const; | |
| 273 | |
| 274 //Comparison | |
| 275 | |
| 276 /// Compares vectors elementwise, using machine epsilon. | |
| 277 bool operator==(const vector& w) const; | |
| 278 /*! \brief Lexicographical order, used for putting into STL sets or | |
| 279 * maps. | |
| 280 * | |
| 281 * Lexicographical ordering of vectors. Vectors of smaller | |
| 282 * dimension are smaller. | |
| 283 */ | |
| 284 bool operator<(const vector& w) const; | |
| 285 | |
| 286 //More complex operations | |
| 287 ///Euclidean norm. | |
| 288 double norm() const; | |
| 289 | |
| 290 friend class vector_view; | |
| 291 friend class matrix; | |
| 292 private: | |
| 293 /// Pointer to associated GSL vector. | |
| 294 gsl_vector * x; | |
| 295 /// Output precision. | |
| 296 static size_t precsn; | |
| 297 /// Machine epsilon | |
| 298 static const double eps; | |
| 299 }; | |
| 300 | |
| 301 | |
| 302 /*! \brief A vector that doesn't own its data; rather, points to data owned | |
| 303 * by another vector. | |
| 304 */ | |
| 305 class vector_view : public vector{ | |
| 306 friend class vector; | |
| 307 friend class matrix; | |
| 308 public: | |
| 309 ~vector_view(); | |
| 310 /// Assignment to a vector; will share same data in memory as that | |
| 311 /// other vector. | |
| 312 vector_view& operator=(const vector& w); | |
| 313 /// Assignment, useful for modifying chunks of vectors all at once. | |
| 314 vector_view& operator=(const vector_view& w); | |
| 315 private: | |
| 316 /// vector_views are invisible to the user and can only | |
| 317 /// be created using slices and matrix or vector indexing. | |
| 318 vector_view(); | |
| 319 /// See above. | |
| 320 vector_view(gsl_vector* y, const slice& s); | |
| 321 /// See above. | |
| 322 vector_view(gsl_matrix* A, const slice& a, const size_t j); | |
| 323 /// See above. | |
| 324 vector_view(gsl_matrix* A, const size_t i, const slice& b); | |
| 325 }; | |
| 326 | |
| 327 | |
| 328 /*! \brief Vector slices corresponding to GNU Octave ranges. | |
| 329 * | |
| 330 * Vector slices are useful for certain vectorisation routines and | |
| 331 * for referring to rows, columns, or equally spaced chunks of a | |
| 332 * vector or a matrix. Set a begininning and an end for the | |
| 333 * stride, and a step size, and this will be the elements that can | |
| 334 * be indexed by this slice in a vector or matrix. E.g. if the | |
| 335 * beginning is 2, the end 8, and the stride is 2, then this | |
| 336 * refers to the slice [2, 4, 6, 8]. The default stride is 1. | |
| 337 * | |
| 338 * Unlike GNU Octave slices, there is no shorthand for referring | |
| 339 * to the entirety of a column or row of a matrix. Must instead | |
| 340 * create a slice for that purpose as s(1,n) where n is the number | |
| 341 * of rows or columns respectively. | |
| 342 */ | |
| 343 class slice{ | |
| 344 public: | |
| 345 /*! \brief Create a vector slice. | |
| 346 * | |
| 347 * \param a - Beginning of the slice. | |
| 348 * \param b - End of the slice. | |
| 349 * \param k - Stride (default is 1). | |
| 350 */ | |
| 351 slice(size_t a, size_t b, size_t k=1); | |
| 352 | |
| 353 /*! For setting the slice parameters anew. | |
| 354 * | |
| 355 * Indices start from 1. | |
| 356 * \param a - Beginning of the slice. | |
| 357 * \param b - End of the slice. | |
| 358 * \param k - Stride (default is 1). | |
| 359 */ | |
| 360 slice operator()(size_t a, size_t b, size_t k=1); | |
| 361 slice set(size_t a, size_t b, size_t k=1) ; | |
| 362 | |
| 363 size_t begin()const {return beg;}; | |
| 364 size_t end() const{return fin;}; | |
| 365 size_t stride() const{return str;}; | |
| 366 | |
| 367 private: | |
| 368 /// Empty and private default constructor. | |
| 369 slice(){}; | |
| 370 size_t beg,fin; | |
| 371 size_t str; | |
| 372 }; | |
| 373 | |
| 374 /// Useful alias, vectors are also points in space. | |
| 375 typedef vector point; //Useful alias. | |
| 376 } | |
| 377 | |
| 378 namespace linalg{ //Non-member functions. | |
| 379 | |
| 380 // I/O | |
| 381 | |
| 382 ///Stream insertion operator | |
| 383 std::ostream& operator<<(std::ostream& os, const vector &v); | |
| 384 ///Stream extraction operator | |
| 385 vector operator>>(std::istream& is, vector& v); | |
| 386 ///Stream insertion operator | |
| 387 std::ostream& operator<<(std::ostream& os, const matrix &M); | |
| 388 ///Stream extraction operator | |
| 389 matrix operator>>(std::istream& is, matrix& v); | |
| 390 | |
| 391 //Some arithmetic functions for comfortable syntax. | |
| 392 | |
| 393 /// Scale a vector. | |
| 394 vector operator*(double a, const vector& v); | |
| 395 /// Euclidean norm of a vector. | |
| 396 double norm(const vector& v); | |
| 397 /// Scale a matrix. | |
| 398 matrix operator*(double a, const matrix& M); | |
| 399 /// Matrix inverse, computed with LU factorisation. | |
| 400 matrix inv(const matrix& A); | |
| 401 /// Return copy of transposed matrix. | |
| 402 matrix T(const matrix& A); | |
| 403 /// Trace. | |
| 404 double tr(const matrix& A); | |
| 405 /// Determinant. | |
| 406 double det(matrix& A); | |
| 407 /// L2 condition number, computed with SVD. | |
| 408 double cond(matrix& A); | |
| 409 | |
| 410 } | |
| 411 | |
| 412 //Inlined functions | |
| 413 namespace linalg{ | |
| 414 inline double& matrix::operator()(const size_t i, const size_t j){ | |
| 415 try{ | |
| 416 if(LUfactored) | |
| 417 delete LUptr; | |
| 418 if(SVDfactored) | |
| 419 gsl_vector_free(SVD); | |
| 420 | |
| 421 LUfactored = false; | |
| 422 SVDfactored = false; | |
| 423 return(*gsl_matrix_ptr(A,i-1,j-1)); | |
| 424 } | |
| 425 catch(indexOutOfRange& exc){ | |
| 426 exc.i = i; | |
| 427 exc.j = j; | |
| 428 exc.m = A -> size1; | |
| 429 exc.n = A -> size2; | |
| 430 throw exc; | |
| 431 } | |
| 432 } | |
| 433 inline const double& matrix::operator()(const size_t i, | |
| 434 const size_t j) const{ | |
| 435 try{ | |
| 436 return(*gsl_matrix_const_ptr(A,i-1,j-1)); | |
| 437 } | |
| 438 catch(indexOutOfRange& exc){ | |
| 439 exc.i = i; | |
| 440 exc.j = j; | |
| 441 exc.m = A -> size1; | |
| 442 exc.n = A -> size2; | |
| 443 throw exc; | |
| 444 } | |
| 445 } | |
| 446 | |
| 447 inline double& vector::operator()(const size_t i) { | |
| 448 try{ | |
| 449 return *gsl_vector_ptr(x,i-1); | |
| 450 } | |
| 451 catch(indexOutOfRange& exc){ | |
| 452 exc.i = i; | |
| 453 exc.n = x -> size; | |
| 454 throw exc; | |
| 455 } | |
| 456 } | |
| 457 | |
| 458 inline const double& vector::operator()(const size_t i) const { | |
| 459 try{ | |
| 460 return *gsl_vector_const_ptr(x,i-1); | |
| 461 } | |
| 462 catch(indexOutOfRange& exc){ | |
| 463 exc.i = i; | |
| 464 exc.n = x -> size; | |
| 465 throw exc; | |
| 466 } | |
| 467 } | |
| 468 } | |
|
12
6e06eb6ec448
Trivial change in include guards.
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
11
diff
changeset
|
469 #endif //__LINALG_HPP__ |