1%%%%%%%%%%%%%%%%%%% 2% XLiFE++ is an extended library of finite elements written in C++ 3% Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin 4% 5% This program is free software: you can redistribute it and/or modify 6% it under the terms of the GNU General Public License as published by 7% the Free Software Foundation, either version 3 of the License, or 8% (at your option) any later version. 9% This program is distributed in the hope that it will be useful, 10% but WITHOUT ANY WARRANTY; without even the implied warranty of 11% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12% GNU General Public License for more details. 13% You should have received a copy of the GNU General Public License 14% along with this program. If not, see <http://www.gnu.org/licenses/>. 15%%%%%%%%%%%%%%%%%%% 16 17\section{The {\classtitle Matrix} class} 18 19The purpose of the {\class Matrix} class is to deal with complex or real dense matrices. It 20is derived from the \emph{std:vector<K>} and templated by the coefficient type; the coefficients 21are stored row by row. 22\vspace{.2cm} 23\begin{lstlisting} 24template<typename K> 25class Matrix : public std::vector<K> 26{private: 27 dimen_t rows_; //!<number of rows 28 public: 29 typedef K type_t; //!<alias on matrix coefficient type 30 typedef typename std::vector<K>::iterator itvk_t; 31 typedef typename std::vector<K>::const_iterator citvk_t; 32 ... 33\end{lstlisting} 34\vspace{.2cm} 35As it is templated, it is also possible to deal with matrix of matrices. But all functions 36are not supported, in particular those involving autocast.\\ 37 38It offers some basic constructors and assignment function: 39\vspace{.2cm} 40\begin{lstlisting}[]{} 41Matrix(); //by default 42Matrix(const dimen_t, const dimen_t); //by dimensions 43Matrix(const dimen_t r, const dimen_t, const K&); //by dimensions and constant 44Matrix(const Matrix<K>&); //by copy 45Matrix(const vector<K>&); //diagonal square matrix from vector 46Matrix(const vector<vector<K> >&); //from a vector of vectors 47Matrix(const vector<Vector<K> >&); //from a vector of Vectors 48Matrix(const dimen_t, const SpecialMatrix); //for special square matrix 49Matrix(const String&); //from file 50Matrix(const char *); //from file 51\end{lstlisting} 52\vspace{.2cm} 53As the {\class Vector} class inherited from {\class std::vector}, the constructor from {\class 54std::vector} of {\class std::vector} works also for {\class Vector} of {\class std::vector}, 55for {\class std::vector} of {\class Vector} and for {\class Vector} of {\class Vector}. The 56special matrices are enumerated by the \emph{enum}: 57\vspace{.2cm} 58\begin{lstlisting}[]{} 59enum SpecialMatrix{ _zeroMatrix=0, _idMatrix, _onesMatrix, _hilbertMatrix }; 60\end{lstlisting} 61\vspace{.2cm} 62Besides the constructors, the following assignment operators are allowed: 63\vspace{.2cm} 64\begin{lstlisting}[]{} 65Matrix<K>& operator=(const Matrix<K>&); 66Matrix<K>& operator=(const Matrix<KK>&); 67Matrix<K>& operator=(const vector<vector<K> >&); 68Matrix<K>& operator=(const vector<Vector<K> >&); 69private: void assign(const KK&); 70\end{lstlisting} 71Be cautious with the generalized templated version \lstinline{Matrix<K>\& operator=(const Matrix<KK>\&)}.\\ 72 73One can access to the dimensions, to the coefficients (both read and read/write access) either 74by iterators ({\cmd begin()} and {\cmd end()}) or by index from 1 to the length, to the rows 75or the columns (both read and read/write access), to the main diagonal: 76\vspace{.2cm} 77\begin{lstlisting}[]{} 78size_t size() const; //overloaded size() 79size_t vsize() const; //column length (number of rows) 80dimen_t numberOfRows() const; //number of rows 81size_t hsize() const; //row length (number of columns) 82dimen_t numberOfColumns() const; //number of columns 83itvk_t begin(); //overloaded iterator begin 84citvk_t begin() const; //overloaded const_iterator begin 85itvk_t end(); //overloaded iterator end 86citvk_t end() const; //overloaded const_iterator begin 87const K& operator()(const dimen_t, const dimen_t) const;//protected access 88K& operator()(const dimen_t, const dimen_t); //read/write access 89Vector<K> row(const dimen_t) const; //return r-th row (r > 0) as Vector 90void row(const dimen_t, const vector<K>&): //set r-th row from a stl vector 91Vector<K> column(const dimen_t) const; //return c-th column as a Vector 92void column(const dimen_t, const vector<K>&); //set c-th column from a stl vector 93Vector<K> diag() const; //return diagonal as a Vector 94Matrix<K>& diag(const K); //reset matrix as a diag. matrix 95Matrix<K>& diag(constvector<K>&); //reset matrix as a diag. matrix 96\end{lstlisting} 97\vspace{.2cm} 98Be care with functions on the diagonal, they change the current matrix to a diagonal matrix.\\ 99 100There are also generalized accesors that allow to extract or set a part of matrix. The row indices and column indices may be given either by lower and upper indices or by vector of indices: 101\vspace{.2cm} 102\begin{lstlisting}[]{} 103Matrix<K> get(const std::vector<Dimen>&, const std::vector<Dimen>&) const; 104Matrix<K> get(const std::pair<Dimen,Dimen>&, const std::pair<Dimen, Dimen>&) const; 105Matrix<K> get(Dimen, const std::pair<Dimen, Dimen>&) const; 106Matrix<K> get(const std::pair<Dimen,Dimen>&, Dimen) const; 107Matrix<K> get(Dimen, const std::vector<Dimen>&) const; 108Matrix<K> get(const std::vector<Dimen>&, Dimen) const; 109Matrix<K> get(Dimen, Dimen, Dimen, Dimen) const; 110Matrix<K> operator()(const std::vector<Dimen>&, const std::vector<Dimen>&) const; 111Matrix<K> operator()(const std::pair<Dimen,Dimen>&, const std::pair<Dimen, Dimen>&) const; 112Matrix<K> operator()(Dimen, const std::pair<Dimen, Dimen>&) const; 113Matrix<K> operator()(const std::pair<Dimen,Dimen>&, Dimen) const; 114Matrix<K> operator()(Dimen, const std::vector<Dimen>&) const; 115Matrix<K> operator()(const std::vector<Dimen>&, Dimen) const; 116Matrix<K> operator()(Dimen, Dimen, Dimen, Dimen) const; 117void set(const std::vector<Dimen>&, const std::vector<Dimen>&, const std::vector<K>&); 118void set(const std::pair<Dimen,Dimen>&, const std::pair<Dimen, Dimen>&, const std::vector<K>&); 119void set(Dimen, const std::pair<Dimen, Dimen>&, const std::vector<K>&); 120void set(const std::pair<Dimen,Dimen>&, Dimen, const std::vector<K>&); 121void set(Dimen, const std::vector<Dimen>&, const std::vector<K>&); 122void set(const std::vector<Dimen>&, Dimen, const std::vector<K>&); 123void set(Dimen, Dimen, Dimen, Dimen, const std::vector<K>&); 124\end{lstlisting} 125\vspace{.2cm} 126 127Through overloaded operators, the {\class Matrix} class provides the following (internal) algebraic 128operations: 129\vspace{.2cm} 130\begin{lstlisting}[]{} 131Matrix<K>& operator+=(const Matrix<KK>&);//add a matrix to the matrix 132Matrix<K>& operator+=(const K&); //add a scalar to the matrix 133Matrix<K>& operator-=(const Matrix<KK>&);//substract a matrix to the matrix 134Matrix<K>& operator-=(const K&); //substract a scalar to the matrix 135Matrix<K>& operator*=(const KK&); //multiply by a scalar the matrix 136Matrix<K>& operator/=(const KK&); //divide by a scalar the matrix 137\end{lstlisting} 138\vspace{.2cm} 139Be careful, the operations such as \lstinline{Matrix<Matrix<K> > += Matrix<K>} does not work.\\ 140 141Moreover, the {\class Matrix} class provide a lot of functionalities: 142\vspace{.2cm} 143\begin{lstlisting}[]{} 144bool shareDims(const Matrix<KK>&) const; //true if matrices bear same dimensions 145bool isSymmetric() const; //symmetric matrix test 146bool isSkewSymmetric() const; //skew-symmetric matrix test 147bool isSelfAdjoint() const; //self adjoint matrix test 148bool isSkewAdjoint() const; //skew adjoint matrix test 149Matrix<K>& transpose(); //matrix self-transpostion 150Matrix<K>& adjoint(); //matrix self-adjoint 151Matrix<Real> real() const; //real part of a matrix 152Matrix<Real> imag() const; //imaginary part of a matrix 153\end{lstlisting} 154\vspace{.2cm} 155and input/output and error utilities: 156\vspace{.2cm} 157\begin{lstlisting}[]{} 158void loadFromFile(const char*); 159void saveToFile(const char*); 160istream& operator >> (istream&, Matrix<K>&); 161ostream& operator<<(ostream&, const Matrix<Real>&); 162ostream& operator<<(ostream&, const Matrix<Complex>&); 163void print(ofstream&); 164void print() const; 165 166void nonSquare(const String&, const size_t, const size_t) const; 167void mismatchDims(const String& ,const size_t, const size_t) const; 168void divideByZero(const String&) const; 169void isSingular() const; 170void complexCastWarning(const String&, const size_t, const size_t) const; 171\end{lstlisting} 172\vspace{.2cm} 173 174Finally, this class provides almost algebraic operations as external functions, with specialized 175versions allowing automatic cast from real to complex type: 176\begin{lstlisting} 177bool operator==(const Matrix<K>&,const Matrix<K>&); //matrix comparison 178bool operator!=(const Matrix<K>&,const Matrix<K>&); //matrix comparison 179 180Matrix<K> operator+(const Matrix<K> &); //unary operator+ 181Matrix<K> operator-(const Matrix<K> &); //unary operator- 182 183Matrix<K> operator+(const Matrix<K>&, const Matrix<K>&);//matrix + matrix 184Matrix<K> operator+(const KK&, const Matrix<K>&); //scalar + matrix 185Matrix<K> operator+(const Matrix<K>&, const KK&) //matrix + scalar 186Matrix<Complex> operator+(const Complex&, const Matrix<Real>&); 187Matrix<Complex> operator+(const Matrix<Real>&, const Complex&); 188Matrix<Complex> operator+(const Matrix<Real>&, const Matrix<Complex>&); 189Matrix<Complex> operator+(const Matrix<Complex>&, const Matrix<Real>&); 190 191Matrix<K> operator-(const Matrix<K>&, const Matrix<K>&); //matrix - matrix 192Matrix<K> operator-(const KK&, const Matrix<K>&); //scalar - matrix 193Matrix<K> operator-(const Matrix<K>&, const KK&) //matrix - scalar 194Matrix<Complex> operator-(const Complex&, const Matrix<Real>&); 195Matrix<Complex> operator-(const Matrix<Real>&, const Complex&); 196Matrix<Complex> operator-(const Matrix<Real>&, const Matrix<Complex>&); 197Matrix<Complex> operator-(const Matrix<Complex>&, const Matrix<Real>&); 198 199Matrix<K> operator*(const KK&, const Matrix<K>&); //scalar * matrix 200Matrix<K> operator*(const Matrix<K>& ,const KK&); //matrix * scalar 201Matrix<Complex> operator*(const Complex&, const Matrix<Real>&); 202Matrix<Complex> operator*(const Matrix<Real>&, const Complex&); 203Matrix<K> operator/(const Matrix<K>&, const KK&); //matrix / scalar 204Matrix<Complex> operator/(const Matrix<Real>&, const Complex&); 205 206Vector<K> operator*(const Matrix<K>&, const Vector<V>&); //matrix * vector 207Vector<Complex> operator*(const Matrix<Real>&, const Vector<Complex>&); 208Vector<Complex> operator*(const Matrix<Complex>&, const Vector<Real>&); 209 210Matrix<K> operator*(const Matrix<K>&, const Matrix<K>&); //matrix * matrix 211Matrix<Complex> operator*(const Matrix<Real>&, const Matrix<Complex>&); 212Matrix<Complex> operator*(const Matrix<Complex>&, const Matrix<Real>&); 213 214Matrix<K> transpose(const Matrix<K>&); //transpose matrix 215void diag(Matrix<K>&); //change matrix into its diagonal 216Matrix<Real> real(const Matrix<Complex>&); //real part 217Matrix<Real> imag(const Matrix<Complex>&); //imaginary part 218Matrix<Real> conj(const Matrix<Real>&); //conjugate complex matrix 219Matrix<Complex> cmplx(const Matrix<Real>&); //cast a real matrix to complex one 220Matrix<Complex> conj(const Matrix<Complex>&);//conjugate complex matrix 221\end{lstlisting} 222\vspace{.2cm} 223Some of algebraic operations are based on general templated functions using iterators on vector 224or matrix: 225\begin{lstlisting}[]{} 226void transpose(const dimen_t, const dimen_t, M_it, MM_it); //transposition 227void adjoint(const dimen_t, const dimen_t, M_it, MM_it); //adjoint 228void matvec(M_it, const V1_it, const V1_it, V2_it, V2_it); //matrix * vector 229V2_it matvec(const Matrix<K>&, const V1_it, const V2_it); //matrix * vector 230void vecmat(M_it, const V1_it, const V1_it, V2_it, V2_it); //vector * matrix 231V2_it vecmat(const Matrix<K>&, const V1_it, const V2_it); //vector * matrix 232void matmat(A_it, const dimen_t, B_it, const dimen_t, const dimen_t, R_it); //matrix * matrix 233\end{lstlisting} 234\vspace{.2cm} 235\emph{Use these algebraic operations only with \emph{Matrix<real\_t>} and \emph{Matrix<complex\_t>}. 236They do not all work for more complicated structure such that \emph{Matrix<Matrix<K> >}.}\\ 237\\ 238If the Eigen library is available (it should be), QR and SVD decomposition are available: 239\begin{lstlisting}[]{} 240void qr(const Matrix<T>&, Matrix<T>&, Matrix<T>&); 241void svdMat(const Matrix<T>&, Matrix<T>&, Vector<T>&, Matrix<T>&, real_t eps=0); 242\end{lstlisting} 243\vspace{.2cm} 244 245Useful aliases hiding template syntax, are available for end users: 246\begin{lstlisting} 247typedef Matrix<Real> RealMatrix, 248typedef Matrix<Complex> ComplexMatrix; 249typedef Matrix<Matrix<Real> > RealMatrices; 250typedef Matrix<Matrix<Complex> > ComplexMatrices; 251\end{lstlisting} 252They are defined in \emph{users\_typedef.hpp} ({\lib init} library). 253 254\displayInfos{library=utils, header=Matrix.hpp, implementation=Matrix.cpp, test=test\_Matrix.cpp, 255header dep={config.h, Vector.hpp, Trace.hpp, Messages.hpp, String.hpp, Algorithms.hpp, complexUtils.hpp}} 256 257\vspace{5mm} 258