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