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{Matrix}
18
19The purpose of the {\class Matrix} class is mainly to deal with complex or real \textbf{dense}
20matrices. In particular, this class is used in the definition of the user functions (see the
21section Function).  This class is compliant with the {\class Vector} class. Although, it can
22deal with matrices of anything, it is only fully functional for real or complex matrices:
23\vspace{.2cm}
24\begin{lstlisting}[]{}
25Matrix<Real> rA;         // an empty matrix
26Matrix<Real> rB(3,2);    // a 3@$\times$@2 zeros matrix
27Matrix<Real> rC(3,2,1);  // a 3@$\times$@2 ones matrix
28Vector<Real> w(3,2.5);   // w=[2.5 2.5 2.5]
29Complex i(0,1);          // the complex i
30Matrix<Real> cA(3,2,i);  // a 3@$\times$@2 i matrix
31\end{lstlisting}
32\vspace{.2cm}
33
34It is possible to construct diagonal matrix from a {\class Vector} or a matrix from a {\class
35Vector} of {\class Vector}, to load (and save) a matrix from a file and to construct particular
36matrices ({\tt \_zeroMatrix, \_onesMatrix, \_idMatrix, \_hilbertMatrix}):
37\vspace{.2cm}
38\begin{lstlisting}[]{}
39Vector<Real> u(3,2.);              // vector [2. 2. 2.]
40Matrix<Real> rA(u);                // 3@$\times$@3 matrix with u as diagonal
41Matrix<Real> rB("mat.dat");        // matrix loaded from "mat.dat" file
42Matrix<Real> rO(3,_zeroMatrix);    // a 3@$\times$@3 zeros matrix
43Matrix<Real> r1(3,_onesMatrix);    // a 3@$\times$@3 ones matrix
44Matrix<Real> rI(3,_idMatrix);      // a 3@$\times$@3 identity matrix
45Matrix<Real> rH(3,_hilbertMatrix); // the 3@$\times$@3 Hilbert matrix
46\end{lstlisting}
47\vspace{.2cm}
48Construction of complex matrix from real data are allowed (automatic cast). But the contrary
49is not. \\
50
51There are some functions to access to the matrix properties:
52\vspace{.2cm}
53\begin{lstlisting}[]{}
54numberOfRows(), numberOfColumns()
55isSymmetric(), isSkewSymmetric(), isSelfAdjoint(), isSkewAdjoint()
56\end{lstlisting}
57\vspace{.2cm}
58ans some utilities to access to a coefficient, a row or a column or the diagonal of the matrix
59:
60\begin{lstlisting}[]{}
61Matrix<Real> A(2,2,1);   //a 2@$\times$@2 ones matrix
62A(1,1)=2.;                 //change the coefficient A11
63Vector<Real> r=A.row(1); //first row of A
64r=A.column(2);             //second column of A
65r=A.diag();                //diagonal of A
66A.column(1,r);             //assign a vector to the first column
67A.row(2,r);                //assign a vector to the second row
68A.diag(r);                 //assign a vector to the diagonal
69\end{lstlisting}
70\vspace{.2cm}
71All these functions support automatic cast from real to complex but not the contrary. \\
72Advanced users can use member functions {\cmd begin()} and {\cmd end()} returning respectively
73iterators (or const iterators) to the beginning and the end of the Matrix. The data values of the matrix
74are stored according to the C convention, i.e. row-wise. \\
75
76There are also generalized access tools either to extract submatrix ({\cmd get()} or {\cmd operator()}) or to set submatrix of matrix ({\cmd set()}):
77\vspace{.2cm}
78\begin{lstlisting}[]{}
79Matrix<Real> M(3,3);
80for(Number i=1;i<=3;i++)
81  for(Number j=1;j<=3;j++)
82	   M(i,j)=i+j;                //M=[2 3 4; 3 4 5; 4 5 6]
83Matrix<Real> N= M.get(2,3,2,3); //N=[4 5; 5 6]
84//N=M(2,3,2,3) gives the same
85Vector<Number> is(2);
86is(1)=1;is(2)=3;
87N= M.get(is,is);                //N=[2 4; 4 6]
88//N=M(is,is) gives the same
89Matrix<Real> Z(2,2,0);          //Z=[0 0; 0 0]
90M.set(1,2,1,2,Z);               //M=[0 0 4; 0 0 5; 4 5 6]
91Matrix<Real> U(2,2,1);          //U=[1 1; 1 1]
92M.set(is,is,Z);                 //M=[1 0 1; 0 0 5; 1 5 1]
93\end{lstlisting}
94\vspace{.2cm}
95Other syntaxes are proposed, see the developer's documentation.\\
96\vspace{.2cm}
97
98Besides, the {\class Matrix} class proposes some transformations either as internal functions
99or external functions:
100\vspace{.2cm}
101\begin{lstlisting}[]{}
102Matrix<Real> A(2,2,1),B;
103Matrix<Complex> C(2,2,i),D;
104A.transpose();     // self transposition of A
105B=transpose(A);    // transposition of A, A not changed
106C.adjoint();       // self transposition and conjugate C
107D=adjoint(A);      // transpose and conjugate, C not changed
108B=diag(A);         // from diagonal of A to a diagonal matrix
109A=real(C);         // real part of C
110B=imag(C);         // imaginary part of C
111D=conj(C);         // conjugate of C
112D=cmplx(A);        // forced casting from real to complex
113Real n2=norm2(A);      //Frobenius norm
114Real ninf=norminfty(A);//infinite norm
115\end{lstlisting}
116\vspace{.2cm}
117
118Standard algebraic operations (+=,-=,*=,/=,+,-,*,/) are supported by the {\class Matrix} class.
119Some shortcuts are also possible, for instance a matrix plus a scalar, a scalar plus a matrix,
120... Automatic cast from real to complex is supported. There is no comparison operator.\\
121\\
122The {\class Matrix} proposes some solvers:
123\begin{lstlisting}[deletekeywords={[3] row}]
124Matrix<Real> A
125RealVector B;
126...
127// solve AX=B or AXs=Bxs using Gauss reduction
128gaussSolver(A, B, piv, row);
129gaussMultipleSolver(A, B, nbrhs, piv, row);
130//inverse of a square matrix
131RealMatrix invA=inverse(A);
132// QR factorization
133RealMatrix Q,R;
134qr(A,Q,R);
135// SVD factorization A= U S V*, if A a (m,n)-matrix, U is a (m,r)-matrix, V a (n,r)-matrix and S a r-vector where r=min(m,n)
136RealMatrix U, V;
137RealVector S;   //singular values
138svdMat(A, U, S, V);
139\end{lstlisting}
140\vspace{.2cm}
141\begin{remark}
142	QR and SVD are available only if Eigen library is set on.
143\end{remark}
144\begin{remark}
145It is also possible to deal with matrix of matrices, for instance:
146\vspace{.2cm}
147\begin{lstlisting}[]{}
148Matrix<Real> ones(2,2,1);          //a 2@$\times$@2 ones matrix
149Matrix<Matrix<Real> > A(2,2,ones); //a 2@$\times$@2 matrix of 2@$\times$@2 ones matrix
150\end{lstlisting}
151\vspace{.2cm}
152but all operations are not supported for such matrices!\\
153\end{remark}
154
155To avoid explicit templates in user program, the following aliases are provided:
156\begin{itemize}[itemsep=-5pt]
157\item {\class RealMatrix} stands for {\class Matrix<Real>},
158\item {\class ComplexMatrix} stands for {\class Matrix<Complex>}.
159\item {\class RealMatrices} stands for {\class Matrix<Matrix<Real> >},
160\item {\class ComplexMatrices} stands for {\class Matrix<Matrix<Complex> >}.
161\end{itemize}
162
163