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 MatrixStorage} class} 18 19The {\class MatrixStorage} is the abstract base class of all supported matrix storages. Vectors carrying storage pointers are defined in child classes and the vector storing the matrix entries is defined outside the storage (see {\class LargeMatrix} class for instance). Various matrix storages are defined by two characteristics : 20\begin{itemize} 21\item a StorageType : \verb?_dense? for dense storage, \verb?_cs? for compressed sparse storage (CSR/CSC) and \verb?_skyline? for skyline strorage 22\item an AccessType : \verb?_row? for a row wise access, \verb?_col? for a column wise access,\verb?_dual? for a row wise access for lower triangular part of the matrix and a column wise access for the upper triangular part of the matrix, and \verb?_sym? if the storage is symmetric (only the lower triangular part is described with a row wise access). 23\end{itemize} 24\begin{lstlisting} 25class MatrixStorage 26{protected: 27 StorageType storageType_; // storage type 28 AccessType accessType_; // access type 29 Number nbRows_; // number of rows 30 Number nbCols_; // number of columns 31 Number nbObjectsSharingThis_; // number of matrices sharing "this" storage 32 ... 33\end{lstlisting} 34\vspace{.1cm} 35The number of rows and columns are counted in value type of the matrix.\\ \\ 36The class provides some constructors which only initialize members of the class. 37\vspace{.1cm} 38\begin{lstlisting}[] 39MatrixStorage(); 40MatrixStorage(StorageType, AccessType at = _noAccess); 41MatrixStorage(StorageType, AccessType, Number, Number); 42virtual ~MatrixStorage(); 43\end{lstlisting} 44\vspace{.1cm} 45Protected members may be read by accessors : 46\vspace{.1cm} 47\begin{lstlisting}[] 48String name() const; 49AccessType accessType() const; 50Number nbOfRows() const; 51Number nbOfColumns() const; 52Number numberOfObjects() const 53void objectPlus() ; 54void objectMinus(); 55\end{lstlisting} 56\vspace{.1cm} 57The class provides virtual member functions returning the sizes of storage (number of matrix entries stored) : 58\vspace{.1cm} 59\begin{lstlisting}[] 60virtual Number size() const = 0; 61Number diagonalSize() const; 62virtual Number lowerPartSize() const=0; 63virtual Number upperPartSize() const=0; 64\end{lstlisting} 65\vspace{.2cm} 66Most of member functions depend on type of storage, so they are virtual member functions. Besides, as storages do not depend on value type of matrix, the storage classes are not templated by the matrix value type. Because virtual template functions are not allowed in C++, all the value type specializations are explicit in base class and in child classes. Up to now, the supported specialisations are \verb?Real?, \verb?Complex?, \verb?Matrix<Real>? and \verb?Matrix<Complex>?. We describe here only the \verb?Real? specialisation. \\ \\ 67There are two functions related to matrix entry position in storage, the first one gives the position of entry $(i,j)$ with $i,j\ge 1$ and the second one compute the positions of the entries of the submatrix given by its row indices and column indices. These functions give \textbf{positions from index 1} if the entry is inside the storage else the returned position is 0. In the second function, it is possible to activate an error when entry index is outside the storage. 68In both functions, the symmetry has to be specified if the matrix has a symmetry property, in that case the position is always in lower triangular part. Be care with the fact that a matrix with no symmetry may have a symmetric storage! 69\vspace{.1cm} 70\begin{lstlisting}[] 71virtual Number pos(Number i, Number j, SymType sym=_noSymmetry) const; 72virtual void positions(const std::vector<Number>&, const std::vector<Number>&, 73 std::vector<Number>&, bool errorOn=true, SymType=_noSymmetry) const; 74\end{lstlisting} 75\vspace{.2cm} 76To print, to save to file or load from file large matrices, the following functions are proposed : 77\vspace{.1cm} 78\begin{lstlisting}[] 79virtual void printEntries(std::ostream&, const std::vector<Real>&, 80 Number vb = 0, SymType sym = _noSymmetry) const; 81template<typename T> 82 void printDenseMatrix(std::ostream&, const std::vector<T>&, SymType s=_noSymmetry) const; 83template<typename T> 84 void printCooMatrix(std::ostream&, const std::vector<T>&, SymType s=_noSymmetry) const; 85virtual void loadFromFileDense(std::istream&, std::vector<Real>&, SymType, bool); 86virtual void loadFromFileCoo(std::istream&, std::vector<Real>&, SymType, bool); 87virtual void visual(std::ostream&) const; 88\end{lstlisting} 89\vspace{.1cm} 90The {\cmd printEntries} function print on output stream the stored matrix entries in a well suited presentation, the {\cmd printDenseMatrix} function print the matrix entries in a dense format as it was dense and the {\cmd printCooMatrix} function print the matrix entries in a coordinates format ($i,j,A_{ij}$. A $m\times n$ matrix saved as dense produced a file like : 91\begin{verbatim} 92 A11 A12 ... A1n re(A11) im(A11) ... re(A1n) im(A1n) 93 A21 A22 ... A2n if complex values re(A21) im(A21) ... re(A2n) im(A2n) 94 ... ... 95 Am1 Am2 ... Amn re(Am1) im(Am1) ... re(Amn) im(Amn) 96\end{verbatim} 97and when saved as coordinates format : 98\begin{verbatim} 99 i1 j1 Ai1j1 i1 j1 re(Ai1j1) im(Ai1j1) 100 i2 j2 Ai2j2 if complex values i2 j2 re(Ai2j2) im(Ai2j2) 101 ... ... 102\end{verbatim} 103Note that {\cmd printDenseMatrix} and {\cmd printCooMatrix} functions are templated functions. Their coding is not optimal because they use the virtual \verb?pos? function 104and other extern utilities which "interpret" the type of the template : 105\vspace{.1cm} 106\begin{lstlisting}[] 107void printDense(std::ostream&, const Real&, Number); 108void printDense(std::ostream&, const Complex&, Number); 109template<typename T> void printDense(std::ostream&, const Matrix<T>&, Number); 110void printCoo(std::ostream&, const Real&, Number); 111void printCoo(std::ostream&, const Complex&, Number); 112template<typename T> void printCoo(std::ostream&, const Matrix<T>&, Number); 113Number numberOfRows(const Real&); 114Number numberOfRows(const Complex&); 115Number numberOfCols(const Real&); 116Number numberOfCols(const Complex&); 117template<typename T> Number numberOfRows(const Matrix<T>&); 118template<typename T> Number numberOfCols(const Matrix<T>& m); 119\end{lstlisting} 120\vspace{.2cm} 121The \verb?visual? function prints only the matrix entries inside the storage as 'x' characters (d for diagonal entries): 122\begin{lstlisting}[deletekeywords={[3] size}] 123matrix storage :row_compressed sparse (csr,csc), size = 49, 9 rows, 9 columns (shared by 0 objects). 124 123456789 125 1 dx.xx.... 126 2 xdxxxx... 127 3 .xd.xx... 128 4 xx.dx.xx. 129 5 xxxxdxxxx 130 6 .xx.xd.xx 131 7 ...xx.dx. 132 8 ...xxxxdx 133 9 ....xx.xd 134 123456789 135\end{lstlisting} 136\vspace{.1cm} 137The {\class MatrixStorage} class provides the product matrix $\times$ vector and the product vector $\times$ matrix: 138\vspace{.1cm} 139\begin{lstlisting}[] 140virtual void multMatrixVector(const std::vector<Real>&, const std::vector<Real>&, std::vector<Real>&, SymType) const; 141virtual void multMatrixVector(const std::vector< Matrix<Real> >&, const std::vector< Vector<Real> >&, std::vector< Vector<Real> >&, SymType) const; 142virtual void multVectorMatrix(const std::vector<Real>&, const std::vector<Real>&, std::vector<Real>&, SymType) const; 143virtual void multVectorMatrix(const std::vector< Matrix<Real> >&, const std::vector< Vector<Real> >&, std::vector< Vector<Real> >&, SymType) const; 144\end{lstlisting} 145\vspace{.1cm} 146Note that mixing real vector/matrix and complex vector/matrix is also provided. Each child class implements the best algorithm for the product.\\ 147The addition of two matrices are supplied with: 148\vspace{.1cm} 149\begin{lstlisting}[] 150virtual void addMatrixMatrix(const std::vector<Real>&, const std::vector<Real>&, 151 std::vector<Real>&) const; 152\end{lstlisting} 153\vspace{.1cm} 154Similar to multiplication of matrix-vector, this function is implemented in each child class. 155The {\class MatrixStorage} class also provides the matrix factorizations 156\vspace{.1cm} 157\begin{lstlisting}[] 158virtual void ldlt(std::vector<Real>&, std::vector<Real>&, 159 const SymType sym = _symmetric) const; 160virtual void ldlt(std::vector<Complex>&, std::vector<Complex>&, 161 const SymType sym = _symmetric) const; 162virtual void lu(std::vector<Real>&, std::vector<Real>&, 163 const SymType sym = _noSymmetry) const; 164virtual void lu(std::vector<Complex>&, std::vector<Complex>&, 165 const SymType sym = _noSymmetry) const; 166virtual void ldlstar(std::vector<Real>&, std::vector<Real>&) const; 167virtual void ldlstar(std::vector<Complex>&, std::vector<Complex>&) const; 168\end{lstlisting} 169\vspace{.1cm} 170as well as some special solvers 171\vspace{.1cm} 172\begin{lstlisting}[] 173virtual void lowerD1Solver(const std::vector<Real>&, std::vector<Real>&, 174 std::vector<Real>&) const; 175virtual void diagonalSolver(const std::vector<Real>&, std::vector<Real>&, 176 std::vector<Real>&) const 177virtual void upperD1Solver(const std::vector<Real>&, std::vector<Real>&, 178 std::vector<Real>&, const SymType sym = _noSymmetry) const; 179virtual void upperSolver(const std::vector<Real>&, std::vector<Real>&, 180 std::vector<Real>&, const SymType sym = _noSymmetry) const; 181virtual void sorDiagonalMatrixVector(const std::vector<Real>&, 182 const std::vector<Real>&, std::vector<Real>&, const Real) const; 183virtual void sorLowerMatrixVector(const std::vector<Real>&, 184 const std::vector<Real>&, std::vector<Real>&, 185 const Real, const SymType = _noSymmetry) const; 186virtual void sorUpperMatrixVector(const std::vector<Real>&, 187 const std::vector<Real>&, std::vector<Real>&, 188 const Real, const SymType = _noSymmetry) const; 189virtual void sorDiagonalSolve(const std::vector<Real>&, 190 const std::vector<Real>&, std::vector<Real>&, const Real) const; 191virtual void sorLowerSolve(const std::vector<Real>&, 192 const std::vector<Real>&, std::vector<Real>&, const Real) const; 193virtual void sorUpperSolve(const std::vector<Real>&, 194 const std::vector<Real>&, std::vector<Real>&, 195 const Real, const SymType sym = _noSymmetry) const; 196\end{lstlisting} 197\vspace{.1cm} 198Each descendant will implements the specifiic algorithm for these virtual functions. 199Values on the diagonal of a matrix is set with the function 200\vspace{.1cm} 201\begin{lstlisting}[] 202virtual void setDiagValue(std::vector<Real>&, const Real); 203\end{lstlisting} 204\vspace{.1cm} 205\displayInfos{library=largeMatrix, header=MatrixStorage.hpp, implementation=MatrixStorage.cpp, test={test\_LargeMatrix...Storage.cpp}, header dep={config.h, utils.h}} 206 207