| [325390] | 1 | /* | 
|---|
| [cca9ef] | 2 | * RealSpaceMatrix.hpp | 
|---|
| [325390] | 3 | * | 
|---|
|  | 4 | *  Created on: Jun 25, 2010 | 
|---|
|  | 5 | *      Author: crueger | 
|---|
|  | 6 | */ | 
|---|
|  | 7 |  | 
|---|
| [cca9ef] | 8 | #ifndef REALSPACEMATRIX_HPP_ | 
|---|
|  | 9 | #define REALSPACEMATRIX_HPP_ | 
|---|
| [325390] | 10 |  | 
|---|
| [986ed3] | 11 | #include <iosfwd> | 
|---|
| [e4fe8d] | 12 | #include "Helpers/defs.hpp" | 
|---|
| [325390] | 13 |  | 
|---|
| [cd82ec] | 14 | /** | 
|---|
|  | 15 | * Simple class to store matrices with a few extra functions | 
|---|
|  | 16 | */ | 
|---|
|  | 17 |  | 
|---|
| [325390] | 18 | class Vector; | 
|---|
| [fee079] | 19 | class MatrixContent; | 
|---|
| [325390] | 20 |  | 
|---|
| [cca9ef] | 21 | /** 3x3 Matrix class. | 
|---|
|  | 22 | * This class has some specific features for 3D space such as rotations or | 
|---|
|  | 23 | * hard-coded determinants. | 
|---|
|  | 24 | */ | 
|---|
|  | 25 | class RealSpaceMatrix | 
|---|
| [325390] | 26 | { | 
|---|
| [cca9ef] | 27 | friend Vector operator*(const RealSpaceMatrix&,const Vector&); | 
|---|
| [325390] | 28 | public: | 
|---|
| [cca9ef] | 29 | RealSpaceMatrix(); | 
|---|
| [cd82ec] | 30 |  | 
|---|
|  | 31 | /** | 
|---|
|  | 32 | * construct a matrix from a 3x3 double array that contains all Elements. | 
|---|
|  | 33 | * | 
|---|
|  | 34 | * The elements are laid out in the following way: | 
|---|
|  | 35 | * array -> matrix | 
|---|
|  | 36 | * 0 -> (0,0) | 
|---|
|  | 37 | * 1 -> (1,0) | 
|---|
|  | 38 | * 2 -> (2,0) | 
|---|
|  | 39 | * 3 -> (0,1) | 
|---|
|  | 40 | * 4 -> (1,1) | 
|---|
|  | 41 | * 5 -> (2,1) | 
|---|
|  | 42 | * 6 -> (0,2) | 
|---|
|  | 43 | * 7 -> (1,2) | 
|---|
|  | 44 | * 8 -> (2,2) | 
|---|
|  | 45 | * | 
|---|
|  | 46 | */ | 
|---|
| [cca9ef] | 47 | RealSpaceMatrix(const double*); | 
|---|
|  | 48 | RealSpaceMatrix(const RealSpaceMatrix&); | 
|---|
|  | 49 | RealSpaceMatrix(const MatrixContent&); | 
|---|
|  | 50 | virtual ~RealSpaceMatrix(); | 
|---|
| [325390] | 51 |  | 
|---|
| [cd82ec] | 52 | /** | 
|---|
|  | 53 | * Set the matrix to a unit matrix. | 
|---|
|  | 54 | */ | 
|---|
| [9eb7580] | 55 | void setIdentity(); | 
|---|
| [1da5f5] | 56 |  | 
|---|
| [a439e5] | 57 | /** | 
|---|
|  | 58 | * Set all matrix entries to zero. | 
|---|
|  | 59 | */ | 
|---|
| [9eb7580] | 60 | void setZero(); | 
|---|
| [a439e5] | 61 |  | 
|---|
| [31fb1d] | 62 | /** | 
|---|
|  | 63 | * Sets all matrix corresponding to a rotation matrix, | 
|---|
|  | 64 | * first around the x-, then the y-, finally the z-axis | 
|---|
|  | 65 | * with given angles. | 
|---|
|  | 66 | */ | 
|---|
| [9eb7580] | 67 | void setRotation(const double x, const double y, const double z); | 
|---|
| [31fb1d] | 68 |  | 
|---|
| [cd82ec] | 69 | /** | 
|---|
|  | 70 | * Access the matrix at index (i,j) | 
|---|
|  | 71 | */ | 
|---|
| [325390] | 72 | double &at(size_t i, size_t j); | 
|---|
| [cd82ec] | 73 | /** | 
|---|
|  | 74 | * Access the matrix at index (i,j) | 
|---|
|  | 75 | */ | 
|---|
| [436f04] | 76 | const double at(size_t i, size_t j) const; | 
|---|
|  | 77 |  | 
|---|
| [cd82ec] | 78 | /** | 
|---|
|  | 79 | * Set the matrix at index (i,j). | 
|---|
|  | 80 | * | 
|---|
|  | 81 | * Slightly faster than at(i,j)=x | 
|---|
|  | 82 | */ | 
|---|
| [436f04] | 83 | void set(size_t i, size_t j, const double value); | 
|---|
| [325390] | 84 |  | 
|---|
| [3dbb9d] | 85 | /** | 
|---|
|  | 86 | * get the ith row of the matrix as a vector | 
|---|
|  | 87 | */ | 
|---|
|  | 88 | Vector &row(size_t); | 
|---|
|  | 89 | const Vector &row(size_t) const; | 
|---|
|  | 90 |  | 
|---|
|  | 91 | /** | 
|---|
|  | 92 | * get the ith column of the matrix as a vector | 
|---|
|  | 93 | */ | 
|---|
|  | 94 | Vector &column(size_t); | 
|---|
|  | 95 | const Vector &column(size_t) const; | 
|---|
|  | 96 |  | 
|---|
|  | 97 | /** | 
|---|
|  | 98 | * get the diagonal of the matrix as a vector | 
|---|
|  | 99 | */ | 
|---|
|  | 100 | Vector &diagonal(); | 
|---|
|  | 101 | const Vector &diagonal() const; | 
|---|
|  | 102 |  | 
|---|
| [cd82ec] | 103 | /** | 
|---|
|  | 104 | * Calculate the determinant of the matrix | 
|---|
|  | 105 | */ | 
|---|
| [cadbc1] | 106 | double determinant() const; | 
|---|
| [325390] | 107 |  | 
|---|
| [cd82ec] | 108 | /** | 
|---|
|  | 109 | * Calculate the inverse of the matrix. | 
|---|
|  | 110 | * | 
|---|
|  | 111 | * Rather costly, so use precomputation as often as possible. | 
|---|
|  | 112 | */ | 
|---|
| [cca9ef] | 113 | RealSpaceMatrix invert() const; | 
|---|
| [325390] | 114 |  | 
|---|
| [a439e5] | 115 | /** | 
|---|
|  | 116 | * Diagonalizes a matrix and sets its rows to the resulting eigenvalues. | 
|---|
|  | 117 | * The eigenvalues are returned as a vector. | 
|---|
|  | 118 | * | 
|---|
|  | 119 | * Rather costly, so use precomputation as often as possible. | 
|---|
|  | 120 | */ | 
|---|
|  | 121 | Vector transformToEigenbasis(); | 
|---|
|  | 122 |  | 
|---|
|  | 123 | /** | 
|---|
|  | 124 | * Calculate the transpose of the matrix. | 
|---|
|  | 125 | */ | 
|---|
| [cca9ef] | 126 | RealSpaceMatrix transpose() const; | 
|---|
|  | 127 | RealSpaceMatrix &transpose(); | 
|---|
| [a439e5] | 128 |  | 
|---|
| [325390] | 129 | // operators | 
|---|
| [cca9ef] | 130 | RealSpaceMatrix &operator=(const RealSpaceMatrix&); | 
|---|
| [325390] | 131 |  | 
|---|
| [cca9ef] | 132 | const RealSpaceMatrix &operator+=(const RealSpaceMatrix&); | 
|---|
|  | 133 | const RealSpaceMatrix &operator-=(const RealSpaceMatrix&); | 
|---|
|  | 134 | const RealSpaceMatrix &operator*=(const RealSpaceMatrix&); | 
|---|
| [325390] | 135 |  | 
|---|
| [cca9ef] | 136 | const RealSpaceMatrix &operator*=(const double); | 
|---|
| [325390] | 137 |  | 
|---|
| [cca9ef] | 138 | const RealSpaceMatrix operator+(const RealSpaceMatrix&) const; | 
|---|
|  | 139 | const RealSpaceMatrix operator-(const RealSpaceMatrix&) const; | 
|---|
|  | 140 | const RealSpaceMatrix operator*(const RealSpaceMatrix&) const; | 
|---|
| [325390] | 141 |  | 
|---|
| [cca9ef] | 142 | bool operator==(const RealSpaceMatrix&) const; | 
|---|
| [0eb2dc] | 143 |  | 
|---|
| [325390] | 144 | private: | 
|---|
| [cca9ef] | 145 | RealSpaceMatrix(MatrixContent*); | 
|---|
| [3dbb9d] | 146 | void createViews(); | 
|---|
| [fee079] | 147 | MatrixContent *content; | 
|---|
| [8e17d6] | 148 | // we keep around some Vector views of the matrix, to return references | 
|---|
| [3dbb9d] | 149 | Vector* rows_ptr[NDIM]; | 
|---|
|  | 150 | Vector* columns_ptr[NDIM]; | 
|---|
|  | 151 | Vector* diagonal_ptr; | 
|---|
| [325390] | 152 | }; | 
|---|
|  | 153 |  | 
|---|
| [cca9ef] | 154 | const RealSpaceMatrix operator*(const double,const RealSpaceMatrix&); | 
|---|
|  | 155 | const RealSpaceMatrix operator*(const RealSpaceMatrix&,const double); | 
|---|
| [325390] | 156 |  | 
|---|
| [cd82ec] | 157 | /** | 
|---|
|  | 158 | * Takes a symmetric matrix that stores the lower diagonal and produces a | 
|---|
|  | 159 | * full matrix. | 
|---|
|  | 160 | * | 
|---|
|  | 161 | * The array is laid out as follows: | 
|---|
|  | 162 | * | 
|---|
|  | 163 | * array -> matrix | 
|---|
|  | 164 | * 0 -> (0,0) | 
|---|
|  | 165 | * 1 -> (1,0);(0,1) | 
|---|
|  | 166 | * 2 -> (1,1) | 
|---|
|  | 167 | * 3 -> (2,0);(0,2) | 
|---|
|  | 168 | * 4 -> (2,1);(1,2) | 
|---|
|  | 169 | * 5 -> (2,2) | 
|---|
|  | 170 | */ | 
|---|
| [cca9ef] | 171 | RealSpaceMatrix ReturnFullMatrixforSymmetric(const double * const cell_size); | 
|---|
| [d10eb6] | 172 |  | 
|---|
| [cca9ef] | 173 | std::ostream &operator<<(std::ostream&,const RealSpaceMatrix&); | 
|---|
|  | 174 | Vector operator*(const RealSpaceMatrix&,const Vector&); | 
|---|
|  | 175 | Vector& operator*=(Vector&,const RealSpaceMatrix&); | 
|---|
| [c49c96] | 176 |  | 
|---|
| [cca9ef] | 177 | #endif /* REALSPACEMATRIX_HPP_ */ | 
|---|