| 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 <iosfwd>
 | 
|---|
| 12 | #include "Helpers/defs.hpp"
 | 
|---|
| 13 | 
 | 
|---|
| 14 | /**
 | 
|---|
| 15 |  * Simple class to store matrices with a few extra functions
 | 
|---|
| 16 |  */
 | 
|---|
| 17 | 
 | 
|---|
| 18 | class Vector;
 | 
|---|
| 19 | class MatrixContent;
 | 
|---|
| 20 | 
 | 
|---|
| 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
 | 
|---|
| 26 | {
 | 
|---|
| 27 |   friend Vector operator*(const RealSpaceMatrix&,const Vector&);
 | 
|---|
| 28 | public:
 | 
|---|
| 29 |   RealSpaceMatrix();
 | 
|---|
| 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 |    */
 | 
|---|
| 47 |   RealSpaceMatrix(const double*);
 | 
|---|
| 48 |   RealSpaceMatrix(const RealSpaceMatrix&);
 | 
|---|
| 49 |   RealSpaceMatrix(const MatrixContent&);
 | 
|---|
| 50 |   virtual ~RealSpaceMatrix();
 | 
|---|
| 51 | 
 | 
|---|
| 52 |   /**
 | 
|---|
| 53 |    * Set the matrix to a unit matrix.
 | 
|---|
| 54 |    */
 | 
|---|
| 55 |   void setIdentity();
 | 
|---|
| 56 | 
 | 
|---|
| 57 |   /**
 | 
|---|
| 58 |    * Set all matrix entries to zero.
 | 
|---|
| 59 |    */
 | 
|---|
| 60 |   void setZero();
 | 
|---|
| 61 | 
 | 
|---|
| 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 |    */
 | 
|---|
| 67 |   void setRotation(const double x, const double y, const double z);
 | 
|---|
| 68 | 
 | 
|---|
| 69 |   /**
 | 
|---|
| 70 |    * Access the matrix at index (i,j)
 | 
|---|
| 71 |    */
 | 
|---|
| 72 |   double &at(size_t i, size_t j);
 | 
|---|
| 73 |   /**
 | 
|---|
| 74 |    * Access the matrix at index (i,j)
 | 
|---|
| 75 |    */
 | 
|---|
| 76 |   const double at(size_t i, size_t j) const;
 | 
|---|
| 77 | 
 | 
|---|
| 78 |   /**
 | 
|---|
| 79 |    * Set the matrix at index (i,j).
 | 
|---|
| 80 |    *
 | 
|---|
| 81 |    * Slightly faster than at(i,j)=x
 | 
|---|
| 82 |    */
 | 
|---|
| 83 |   void set(size_t i, size_t j, const double value);
 | 
|---|
| 84 | 
 | 
|---|
| 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 | 
 | 
|---|
| 103 |   /**
 | 
|---|
| 104 |    * Calculate the determinant of the matrix
 | 
|---|
| 105 |    */
 | 
|---|
| 106 |   double determinant() const;
 | 
|---|
| 107 | 
 | 
|---|
| 108 |   /**
 | 
|---|
| 109 |    * Calculate the inverse of the matrix.
 | 
|---|
| 110 |    *
 | 
|---|
| 111 |    * Rather costly, so use precomputation as often as possible.
 | 
|---|
| 112 |    */
 | 
|---|
| 113 |   RealSpaceMatrix invert() const;
 | 
|---|
| 114 | 
 | 
|---|
| 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 |    */
 | 
|---|
| 126 |   RealSpaceMatrix transpose() const;
 | 
|---|
| 127 |   RealSpaceMatrix &transpose();
 | 
|---|
| 128 | 
 | 
|---|
| 129 |   // operators
 | 
|---|
| 130 |   RealSpaceMatrix &operator=(const RealSpaceMatrix&);
 | 
|---|
| 131 | 
 | 
|---|
| 132 |   const RealSpaceMatrix &operator+=(const RealSpaceMatrix&);
 | 
|---|
| 133 |   const RealSpaceMatrix &operator-=(const RealSpaceMatrix&);
 | 
|---|
| 134 |   const RealSpaceMatrix &operator*=(const RealSpaceMatrix&);
 | 
|---|
| 135 | 
 | 
|---|
| 136 |   const RealSpaceMatrix &operator*=(const double);
 | 
|---|
| 137 | 
 | 
|---|
| 138 |   const RealSpaceMatrix operator+(const RealSpaceMatrix&) const;
 | 
|---|
| 139 |   const RealSpaceMatrix operator-(const RealSpaceMatrix&) const;
 | 
|---|
| 140 |   const RealSpaceMatrix operator*(const RealSpaceMatrix&) const;
 | 
|---|
| 141 | 
 | 
|---|
| 142 |   bool operator==(const RealSpaceMatrix&) const;
 | 
|---|
| 143 | 
 | 
|---|
| 144 | private:
 | 
|---|
| 145 |   RealSpaceMatrix(MatrixContent*);
 | 
|---|
| 146 |   void createViews();
 | 
|---|
| 147 |   MatrixContent *content;
 | 
|---|
| 148 |   // we keep around some Vector views of the matrix, to return references
 | 
|---|
| 149 |   Vector* rows_ptr[NDIM];
 | 
|---|
| 150 |   Vector* columns_ptr[NDIM];
 | 
|---|
| 151 |   Vector* diagonal_ptr;
 | 
|---|
| 152 | };
 | 
|---|
| 153 | 
 | 
|---|
| 154 | const RealSpaceMatrix operator*(const double,const RealSpaceMatrix&);
 | 
|---|
| 155 | const RealSpaceMatrix operator*(const RealSpaceMatrix&,const double);
 | 
|---|
| 156 | 
 | 
|---|
| 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 |  */
 | 
|---|
| 171 | RealSpaceMatrix ReturnFullMatrixforSymmetric(const double * const cell_size);
 | 
|---|
| 172 | 
 | 
|---|
| 173 | std::ostream &operator<<(std::ostream&,const RealSpaceMatrix&);
 | 
|---|
| 174 | Vector operator*(const RealSpaceMatrix&,const Vector&);
 | 
|---|
| 175 | Vector& operator*=(Vector&,const RealSpaceMatrix&);
 | 
|---|
| 176 | 
 | 
|---|
| 177 | #endif /* REALSPACEMATRIX_HPP_ */
 | 
|---|