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