[fee079] | 1 | /*
|
---|
| 2 | * MatrixContent.hpp
|
---|
| 3 | *
|
---|
| 4 | * Created on: Jul 2, 2010
|
---|
| 5 | * Author: crueger
|
---|
| 6 | */
|
---|
| 7 |
|
---|
| 8 | #ifndef MATRIXCONTENT_HPP_
|
---|
| 9 | #define MATRIXCONTENT_HPP_
|
---|
| 10 |
|
---|
[56f73b] | 11 | // include config.h
|
---|
| 12 | #ifdef HAVE_CONFIG_H
|
---|
| 13 | #include <config.h>
|
---|
| 14 | #endif
|
---|
| 15 |
|
---|
| 16 |
|
---|
[3bc926] | 17 | /** MatrixContent is a wrapper for gsl_matrix.
|
---|
| 18 | *
|
---|
| 19 | * The GNU Scientific Library contaisn some very well written routines for
|
---|
| 20 | * linear algebra problems. However, it's syntax is C and hence it does not
|
---|
| 21 | * lend itself to readable written code, i.e. C = A * B, where A, B, and C
|
---|
| 22 | * are all matrices. Writing code this way is very convenient, concise and
|
---|
| 23 | * also in the same manner as one would type in MatLab.
|
---|
| 24 | * However, with C++ and its feature to overload functions and its operator
|
---|
| 25 | * functions such syntax is easy to create.
|
---|
| 26 | *
|
---|
| 27 | * Hence, this class is a C++ wrapper to gsl_matrix. There already some other
|
---|
| 28 | * libraries, however they fall short for the following reasons:
|
---|
| 29 | * -# gslwrap: not very well written and uses floats instead of doubles
|
---|
| 30 | * -# gslcpp: last change is from 2007 and only very few commits
|
---|
| 31 | * -# o2scl: basically, the same functionality as gsl only written in C++,
|
---|
| 32 | * however it uses GPLv3 license which is inconvenient for MoleCuilder.
|
---|
| 33 | *
|
---|
| 34 | * <h1>Howto use MatrixContent</h1>
|
---|
| 35 | *
|
---|
| 36 | * One should not use MatrixContent directly but either have it as a member
|
---|
| 37 | * variable in a specialized class or inherit from it.
|
---|
| 38 | *
|
---|
[fee079] | 39 | */
|
---|
| 40 |
|
---|
| 41 | #include <gsl/gsl_matrix.h>
|
---|
| 42 |
|
---|
[b4cf2b] | 43 | #include <iosfwd>
|
---|
| 44 |
|
---|
[3da2fb] | 45 | #include "LinearAlgebra/MatrixVector_ops.hpp"
|
---|
| 46 |
|
---|
[60dada] | 47 | class VectorContent;
|
---|
[3bc926] | 48 |
|
---|
[8e9ce1] | 49 | /** Dummy structure to create a unique signature.
|
---|
| 50 | *
|
---|
| 51 | */
|
---|
| 52 | struct MatrixBaseCase{};
|
---|
| 53 |
|
---|
[3bc926] | 54 | class MatrixContent
|
---|
| 55 | {
|
---|
[b4cf2b] | 56 | friend std::ostream & operator<<(std::ostream &ost, const MatrixContent &mat);
|
---|
[cca9ef] | 57 | friend class RealSpaceMatrix;
|
---|
[0d4424] | 58 | friend class LinearSystemOfEquations;
|
---|
| 59 |
|
---|
[3da2fb] | 60 | // matrix vector products
|
---|
| 61 | friend VectorContent const operator*(const VectorContent& vec, const MatrixContent& mat);
|
---|
| 62 | friend VectorContent const operator*(const MatrixContent& mat, const VectorContent& vec);
|
---|
| 63 |
|
---|
| 64 | // unit tests
|
---|
[0d4424] | 65 | friend class MatrixContentTest;
|
---|
| 66 | friend class MatrixContentSymmetricTest;
|
---|
[3bc926] | 67 | public:
|
---|
| 68 | MatrixContent(size_t rows, size_t columns);
|
---|
[8e9ce1] | 69 | MatrixContent(size_t _rows, size_t _columns, MatrixBaseCase);
|
---|
[3bc926] | 70 | MatrixContent(size_t rows, size_t columns, const double *src);
|
---|
| 71 | MatrixContent(gsl_matrix *&src);
|
---|
| 72 | MatrixContent(const MatrixContent &src);
|
---|
| 73 | MatrixContent(const MatrixContent *src);
|
---|
[8e9ce1] | 74 | virtual ~MatrixContent();
|
---|
[3bc926] | 75 |
|
---|
| 76 | // Accessing
|
---|
| 77 | double &at(size_t i, size_t j);
|
---|
| 78 | const double at(size_t i, size_t j) const;
|
---|
| 79 | void set(size_t i, size_t j, const double value);
|
---|
[bbf1bd] | 80 |
|
---|
| 81 | // Initializing
|
---|
[0d4424] | 82 | void setFromDoubleArray(double * x);
|
---|
[bbf1bd] | 83 | void setIdentity();
|
---|
| 84 | void setValue(double _value);
|
---|
| 85 | void setZero();
|
---|
[0d4424] | 86 |
|
---|
| 87 | // Exchanging elements
|
---|
| 88 | bool SwapRows(size_t i, size_t j);
|
---|
| 89 | bool SwapColumns(size_t i, size_t j);
|
---|
| 90 | bool SwapRowColumn(size_t i, size_t j);
|
---|
[3bc926] | 91 |
|
---|
| 92 | // Transformations
|
---|
| 93 | MatrixContent transpose() const;
|
---|
| 94 | MatrixContent &transpose();
|
---|
| 95 | gsl_vector* transformToEigenbasis();
|
---|
[e828c0] | 96 | void sortEigenbasis(gsl_vector *eigenvalues);
|
---|
[3bc926] | 97 |
|
---|
[0d4424] | 98 | // Properties
|
---|
| 99 | bool IsNull() const;
|
---|
| 100 | bool IsPositive() const;
|
---|
| 101 | bool IsNegative() const;
|
---|
| 102 | bool IsNonNegative() const;
|
---|
| 103 | bool IsPositiveDefinite() const;
|
---|
| 104 | double Determinant() const;
|
---|
| 105 |
|
---|
[3bc926] | 106 | // Operators
|
---|
| 107 | MatrixContent &operator=(const MatrixContent &src);
|
---|
| 108 | const MatrixContent &operator+=(const MatrixContent &rhs);
|
---|
| 109 | const MatrixContent &operator-=(const MatrixContent &rhs);
|
---|
| 110 | const MatrixContent &operator*=(const MatrixContent &rhs);
|
---|
| 111 | const MatrixContent operator*(const MatrixContent &rhs) const;
|
---|
[d85c28] | 112 | const MatrixContent &operator&=(const MatrixContent &rhs);
|
---|
| 113 | const MatrixContent operator&(const MatrixContent &rhs) const;
|
---|
[3bc926] | 114 | const MatrixContent &operator*=(const double factor);
|
---|
| 115 | bool operator==(const MatrixContent &rhs) const;
|
---|
| 116 |
|
---|
[35fbef] | 117 | const size_t getRows() const;
|
---|
| 118 | const size_t getColumns() const;
|
---|
[60dada] | 119 |
|
---|
| 120 | VectorContent *getColumnVector(size_t column) const;
|
---|
| 121 | VectorContent *getRowVector(size_t row) const;
|
---|
| 122 | VectorContent *getDiagonalVector() const;
|
---|
| 123 |
|
---|
[0d4424] | 124 | protected:
|
---|
| 125 | double *Pointer(size_t m, size_t n);
|
---|
| 126 | const double *const_Pointer(size_t m, size_t n) const;
|
---|
[3bc926] | 127 |
|
---|
[fee079] | 128 | gsl_matrix * content;
|
---|
[3bc926] | 129 | const size_t rows; // store for internal purposes
|
---|
| 130 | const size_t columns; // store for internal purposes
|
---|
[8e9ce1] | 131 |
|
---|
| 132 | private:
|
---|
[fee079] | 133 | };
|
---|
| 134 |
|
---|
[3bc926] | 135 | const MatrixContent operator*(const double factor,const MatrixContent& mat);
|
---|
| 136 | const MatrixContent operator*(const MatrixContent &mat,const double factor);
|
---|
[b4cf2b] | 137 | std::ostream & operator<<(std::ostream &ost, const MatrixContent &mat);
|
---|
[3bc926] | 138 |
|
---|
[8e9ce1] | 139 | /** Matrix view.
|
---|
| 140 | * Extension of MatrixContent to contain not a gsl_matrix but only a view on a
|
---|
| 141 | * gsl_matrix
|
---|
| 142 | *
|
---|
| 143 | * We need the above MatrixBaseCase here:
|
---|
| 144 | * content, i.e. the gsl_matrix, must not be allocated, as it is just a view
|
---|
| 145 | * with an internal gsl_matrix_view. Hence, MatrixBaseCase is just dummy class
|
---|
| 146 | * to give the constructor a unique signature.
|
---|
| 147 | */
|
---|
| 148 | struct MatrixViewContent : public MatrixContent
|
---|
| 149 | {
|
---|
| 150 | MatrixViewContent(gsl_matrix_view _view) :
|
---|
| 151 | MatrixContent(_view.matrix.size1, _view.matrix.size2, MatrixBaseCase()),
|
---|
| 152 | view(_view)
|
---|
| 153 | {
|
---|
| 154 | content = &view.matrix;
|
---|
| 155 | }
|
---|
| 156 | ~MatrixViewContent(){
|
---|
| 157 | content=0;
|
---|
| 158 | }
|
---|
| 159 | gsl_matrix_view view;
|
---|
| 160 | };
|
---|
| 161 |
|
---|
[fee079] | 162 | #endif /* MATRIXCONTENT_HPP_ */
|
---|