[fc3b67] | 1 | /*
|
---|
| 2 | * gslmatrix.cpp
|
---|
| 3 | *
|
---|
| 4 | * Created on: Jan 8, 2010
|
---|
| 5 | * Author: heber
|
---|
| 6 | */
|
---|
| 7 |
|
---|
[112b09] | 8 | #include "Helpers/MemDebug.hpp"
|
---|
| 9 |
|
---|
[fc3b67] | 10 | using namespace std;
|
---|
| 11 |
|
---|
| 12 | #include "gslmatrix.hpp"
|
---|
| 13 | #include "helpers.hpp"
|
---|
[0a4f7f] | 14 | #include "Helpers/fast_functions.hpp"
|
---|
[fc3b67] | 15 |
|
---|
| 16 | #include <cassert>
|
---|
| 17 | #include <gsl/gsl_linalg.h>
|
---|
| 18 |
|
---|
| 19 | /** Constructor of class GSLMatrix.
|
---|
| 20 | * Allocates GSL structures
|
---|
| 21 | * \param m dimension of matrix
|
---|
| 22 | */
|
---|
| 23 | GSLMatrix::GSLMatrix(size_t m, size_t n) : rows(m), columns(n)
|
---|
| 24 | {
|
---|
| 25 | matrix = gsl_matrix_calloc(rows, columns);
|
---|
| 26 | };
|
---|
| 27 |
|
---|
| 28 | /** Copy constructor of class GSLMatrix.
|
---|
| 29 | * Allocates GSL structures and copies components from \a *src.
|
---|
| 30 | * \param *src source matrix
|
---|
| 31 | */
|
---|
| 32 | GSLMatrix::GSLMatrix(const GSLMatrix * const src) : rows(src->rows), columns(src->columns)
|
---|
| 33 | {
|
---|
| 34 | matrix = gsl_matrix_alloc(rows, columns);
|
---|
| 35 | gsl_matrix_memcpy (matrix, src->matrix);
|
---|
| 36 | };
|
---|
| 37 |
|
---|
| 38 | /** Destructor of class GSLMatrix.
|
---|
| 39 | * Frees GSL structures
|
---|
| 40 | */
|
---|
| 41 | GSLMatrix::~GSLMatrix()
|
---|
| 42 | {
|
---|
| 43 | gsl_matrix_free(matrix);
|
---|
| 44 | rows = 0;
|
---|
| 45 | columns = 0;
|
---|
| 46 | };
|
---|
| 47 |
|
---|
| 48 | /** Assignment operator.
|
---|
| 49 | * \param &rhs right hand side
|
---|
| 50 | * \return object itself
|
---|
| 51 | */
|
---|
| 52 | GSLMatrix& GSLMatrix::operator=(const GSLMatrix& rhs)
|
---|
| 53 | {
|
---|
| 54 | if (this == &rhs) // not necessary here, but identity assignment check is generally a good idea
|
---|
| 55 | return *this;
|
---|
| 56 | assert(rows == rhs.rows && columns == rhs.columns && "Number of rows and columns do not match!");
|
---|
| 57 |
|
---|
| 58 | gsl_matrix_memcpy (matrix, rhs.matrix);
|
---|
| 59 |
|
---|
| 60 | return *this;
|
---|
| 61 | };
|
---|
| 62 |
|
---|
| 63 | /* ============================ Accessing =============================== */
|
---|
| 64 | /** This function sets the matrix from a double array.
|
---|
| 65 | * Creates a matrix view of the array and performs a memcopy.
|
---|
| 66 | * \param *x array of values (no dimension check is performed)
|
---|
| 67 | */
|
---|
| 68 | void GSLMatrix::SetFromDoubleArray(double * x)
|
---|
| 69 | {
|
---|
| 70 | gsl_matrix_view m = gsl_matrix_view_array (x, rows, columns);
|
---|
| 71 | gsl_matrix_memcpy (matrix, &m.matrix);
|
---|
| 72 | };
|
---|
| 73 |
|
---|
| 74 | /** This function returns the i-th element of a matrix.
|
---|
| 75 | * If \a m or \a n lies outside the allowed range of 0 to GSLMatrix::dimension-1 then the error handler is invoked and 0 is returned.
|
---|
| 76 | * \param m row index
|
---|
| 77 | * \param n colum index
|
---|
| 78 | * \return (m,n)-th element of matrix
|
---|
| 79 | */
|
---|
| 80 | double GSLMatrix::Get(size_t m, size_t n)
|
---|
| 81 | {
|
---|
| 82 | return gsl_matrix_get (matrix, m, n);
|
---|
| 83 | };
|
---|
| 84 |
|
---|
| 85 | /** This function sets the value of the \a m -th element of a matrix to \a x.
|
---|
| 86 | * If \a m or \a n lies outside the allowed range of 0 to GSLMatrix::dimension-1 then the error handler is invoked.
|
---|
| 87 | * \param m row index
|
---|
| 88 | * \param m column index
|
---|
| 89 | * \param x value to set element (m,n) to
|
---|
| 90 | */
|
---|
| 91 | void GSLMatrix::Set(size_t m, size_t n, double x)
|
---|
| 92 | {
|
---|
| 93 | gsl_matrix_set (matrix, m, n, x);
|
---|
| 94 | };
|
---|
| 95 |
|
---|
| 96 | /** These functions return a pointer to the \a m-th element of a matrix.
|
---|
| 97 | * If \a m or \a n lies outside the allowed range of 0 to GSLMatrix::dimension-1 then the error handler is invoked and a null pointer is returned.
|
---|
| 98 | * \param m index
|
---|
| 99 | * \return pointer to \a m-th element
|
---|
| 100 | */
|
---|
| 101 | double *GSLMatrix::Pointer(size_t m, size_t n)
|
---|
| 102 | {
|
---|
| 103 | return gsl_matrix_ptr (matrix, m, n);
|
---|
| 104 | };
|
---|
| 105 |
|
---|
| 106 | /** These functions return a constant pointer to the \a m-th element of a matrix.
|
---|
| 107 | * If \a m or \a n lies outside the allowed range of 0 to GSLMatrix::dimension-1 then the error handler is invoked and a null pointer is returned.
|
---|
| 108 | * \param m index
|
---|
| 109 | * \return const pointer to \a m-th element
|
---|
| 110 | */
|
---|
| 111 | const double *GSLMatrix::const_Pointer(size_t m, size_t n)
|
---|
| 112 | {
|
---|
| 113 | return gsl_matrix_const_ptr (matrix, m, n);
|
---|
| 114 | };
|
---|
| 115 |
|
---|
| 116 | /* ========================== Initializing =============================== */
|
---|
| 117 | /** This function sets all the elements of the matrix to the value \a x.
|
---|
| 118 | * \param *x
|
---|
| 119 | */
|
---|
| 120 | void GSLMatrix::SetAll(double x)
|
---|
| 121 | {
|
---|
| 122 | gsl_matrix_set_all (matrix, x);
|
---|
| 123 | };
|
---|
| 124 |
|
---|
| 125 | /** This function sets all the elements of the matrix to zero.
|
---|
| 126 | */
|
---|
| 127 | void GSLMatrix::SetZero()
|
---|
| 128 | {
|
---|
| 129 | gsl_matrix_set_zero (matrix);
|
---|
| 130 | };
|
---|
| 131 |
|
---|
| 132 | /** This function sets the elements of the matrix to the corresponding elements of the identity matrix, \f$m(i,j) = \delta(i,j)\f$, i.e. a unit diagonal with all off-diagonal elements zero.
|
---|
| 133 | * This applies to both square and rectangular matrices.
|
---|
| 134 | */
|
---|
| 135 | void GSLMatrix::SetIdentity()
|
---|
| 136 | {
|
---|
| 137 | gsl_matrix_set_identity (matrix);
|
---|
| 138 | };
|
---|
| 139 |
|
---|
| 140 | /* ====================== Exchanging elements ============================ */
|
---|
| 141 | /** This function exchanges the \a i-th and \a j-th row of the matrix in-place.
|
---|
| 142 | * \param i i-th row to swap with ...
|
---|
| 143 | * \param j ... j-th row to swap against
|
---|
| 144 | */
|
---|
| 145 | bool GSLMatrix::SwapRows(size_t i, size_t j)
|
---|
| 146 | {
|
---|
| 147 | return (gsl_matrix_swap_rows (matrix, i, j) == GSL_SUCCESS);
|
---|
| 148 | };
|
---|
| 149 |
|
---|
| 150 | /** This function exchanges the \a i-th and \a j-th column of the matrix in-place.
|
---|
| 151 | * \param i i-th column to swap with ...
|
---|
| 152 | * \param j ... j-th column to swap against
|
---|
| 153 | */
|
---|
| 154 | bool GSLMatrix::SwapColumns(size_t i, size_t j)
|
---|
| 155 | {
|
---|
| 156 | return (gsl_matrix_swap_columns (matrix, i, j) == GSL_SUCCESS);
|
---|
| 157 | };
|
---|
| 158 |
|
---|
| 159 | /** This function exchanges the \a i-th row and \a j-th column of the matrix in-place.
|
---|
| 160 | * The matrix must be square for this operation to be possible.
|
---|
| 161 | * \param i i-th row to swap with ...
|
---|
| 162 | * \param j ... j-th column to swap against
|
---|
| 163 | */
|
---|
| 164 | bool GSLMatrix::SwapRowColumn(size_t i, size_t j)
|
---|
| 165 | {
|
---|
| 166 | assert (rows == columns && "The matrix must be square for swapping row against column to be possible.");
|
---|
| 167 | return (gsl_matrix_swap_rowcol (matrix, i, j) == GSL_SUCCESS);
|
---|
| 168 | };
|
---|
| 169 |
|
---|
| 170 | /** This function transposes the matrix.
|
---|
| 171 | * Note that the function is extended to the non-square case, where the matrix is re-allocated and copied.
|
---|
| 172 | */
|
---|
| 173 | bool GSLMatrix::Transpose()
|
---|
| 174 | {
|
---|
| 175 | if (rows == columns)// if square, use GSL
|
---|
| 176 | return (gsl_matrix_transpose (matrix) == GSL_SUCCESS);
|
---|
| 177 | else { // otherwise we have to copy a bit
|
---|
| 178 | gsl_matrix *dest = gsl_matrix_alloc(columns, rows);
|
---|
| 179 | for (size_t i=0;i<rows; i++)
|
---|
| 180 | for (size_t j=0;j<columns;j++) {
|
---|
| 181 | gsl_matrix_set(dest, j,i, gsl_matrix_get(matrix, i,j) );
|
---|
| 182 | }
|
---|
[93c8ed] | 183 | gsl_matrix_free(matrix);
|
---|
[fc3b67] | 184 | matrix = dest;
|
---|
| 185 | flip(rows, columns);
|
---|
| 186 | return true;
|
---|
| 187 | }
|
---|
| 188 | };
|
---|
| 189 |
|
---|
| 190 | /* ============================ Properties ============================== */
|
---|
| 191 | /** Checks whether matrix' elements are strictly null.
|
---|
| 192 | * \return true - is null, false - else
|
---|
| 193 | */
|
---|
| 194 | bool GSLMatrix::IsNull()
|
---|
| 195 | {
|
---|
| 196 | return gsl_matrix_isnull (matrix);
|
---|
| 197 | };
|
---|
| 198 |
|
---|
| 199 | /** Checks whether matrix' elements are strictly positive.
|
---|
| 200 | * \return true - is positive, false - else
|
---|
| 201 | */
|
---|
| 202 | bool GSLMatrix::IsPositive()
|
---|
| 203 | {
|
---|
| 204 | return gsl_matrix_ispos (matrix);
|
---|
| 205 | };
|
---|
| 206 |
|
---|
| 207 | /** Checks whether matrix' elements are strictly negative.
|
---|
| 208 | * \return true - is negative, false - else
|
---|
| 209 | */
|
---|
| 210 | bool GSLMatrix::IsNegative()
|
---|
| 211 | {
|
---|
| 212 | return gsl_matrix_isneg (matrix);
|
---|
| 213 | };
|
---|
| 214 |
|
---|
| 215 | /** Checks whether matrix' elements are strictly non-negative.
|
---|
| 216 | * \return true - is non-negative, false - else
|
---|
| 217 | */
|
---|
| 218 | bool GSLMatrix::IsNonNegative()
|
---|
| 219 | {
|
---|
| 220 | return gsl_matrix_isnonneg (matrix);
|
---|
| 221 | };
|
---|
| 222 |
|
---|
| 223 | /** This function performs a Cholesky decomposition to determine whether matrix is positive definite.
|
---|
| 224 | * We check whether GSL returns GSL_EDOM as error, indicating that decomposition failed due to matrix not being positive-definite.
|
---|
| 225 | * \return true - matrix is positive-definite, false - else
|
---|
| 226 | */
|
---|
| 227 | bool GSLMatrix::IsPositiveDefinite()
|
---|
| 228 | {
|
---|
| 229 | if (rows != columns) // only possible for square matrices.
|
---|
| 230 | return false;
|
---|
| 231 | else
|
---|
| 232 | return (gsl_linalg_cholesky_decomp (matrix) != GSL_EDOM);
|
---|
| 233 | };
|
---|
[865272f] | 234 |
|
---|
| 235 |
|
---|
| 236 | /** Calculates the determinant of the matrix.
|
---|
| 237 | * if matrix is square, uses LU decomposition.
|
---|
| 238 | */
|
---|
| 239 | double GSLMatrix::Determinant()
|
---|
| 240 | {
|
---|
| 241 | int signum = 0;
|
---|
| 242 | assert (rows == columns && "Determinant can only be calculated for square matrices.");
|
---|
| 243 | gsl_permutation *p = gsl_permutation_alloc(rows);
|
---|
| 244 | gsl_linalg_LU_decomp(matrix, p, &signum);
|
---|
| 245 | gsl_permutation_free(p);
|
---|
| 246 | return gsl_linalg_LU_det(matrix, signum);
|
---|
| 247 | };
|
---|
| 248 |
|
---|