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{Skyline storages} 18 19With dense storages, skyline storages are the only storages that are compatibles with factorisation methods like LU. Skyline storages consist in storage rows or columns from its first non zero value up to its last non zero value. When lower triangular part is addressed, the last non zero value on a row is the diagonal matrix entry. When upper triangular part is addressed, the last non zero value on a column is the diagonal matrix entry. Only dual and symmetric skyline storages are proposed here (resp. {\class DualSkylineStorage} and {\class SymSkylineStorage} classes). 20\begin{itemize} 21\item Dual skyline stores first the diagonal of the matrix, then the strict lower "triangular" part and the strict upper "triangular" part of the matrix using the following storage vectors ($A$ is a $m\times n$ matrix):\\ 22\mbox{}\hspace{10mm} values=(0, $A_{11}, ...,A_{ii}, ..., A_{il},..., A_{i,p_i}, ..., A_{i,l_i},..., A_{q_j,j}, ..., A_{q_j,c_j}$)\\ 23\mbox{}\hspace{10mm} rowPointer=($s_1,...,s_i,...,s_m,s_{m+1}$)\\ 24\mbox{}\hspace{10mm} colPointer=($t_1,...,t_j,...,t_n,t_{n+1}$)\\ 25where 26\[\begin{array}{l} 27l_i=\min(i-1,n-1),\ c_j=\min(j-1,m-1)\\ 28p_i=\min{0<j<l_i,\ A_{ij}\ne 0} \text{ and } q_j=\min{0<i<c_j,\ A_{ij}\ne 0}\\ 29s_1=0;\ s_i=s_{i-1}+i-s_{i-1}\ \forall i=2,m, \text{ and } t_1=0;\ t_j=s_{j-1}+j-t_{j-1}\ \forall j=2,n 30\end{array} 31\] 32The rowPointer (resp. colPointer) vector gives the position in lower (resp. upper) part of the first non zero entry of a row (resp. a column); $s_{m+1}$ (resp.$t_{n+1}$) gives the number of entries stored in lower part (resp. upper part). \\ \\ 33For instance, the following $5\times 6$ non symmetric matrix 34\[ 35A=\left[ 36\begin{array}{cccccc} 3711 & 12 & 0 & 14 & 0 & 16\\ 380 & 22 & 23 & 0 & 0 & 26 \\ 3931 & 32 & 33 & 0 & 0 & 0 \\ 400 & 0 & 43 & 44 & 0 & 0 \\ 4151 & 0 & 53 & 54 & 55 & 0 42\end{array} 43\right] 44\] 45has the following storage vectors: \\ \\ 46\mbox{}\hspace{10 mm}values = ( 0 11 22 33 44 45\ \ 22\ \ 31 32\ \ 43\ \ 51 0 53 54 \ \ 12\ \ 23\ \ 14 0 0 \ \ 16 26 0 0 ) \\ 47\mbox{}\hspace{10 mm}rowPointer = ( 0 0 0 2 3 7 )\\ 48\mbox{}\hspace{10 mm}colPointer = ( 0 0 1 2 5 5 8 ) \\ 49\\ 50The lengh of row $i$ (resp. column $j$) is given by $s_{i+1}-s_i$ (resp. $t_{t+1}-t_j$). The position in values vector of entry $(i,j)$ is given by the following relations ($1\le i\le m$, $1\le j\le n$): 51\[\begin{array}{l} 52\text {if } i=j\ \ \text{adr}(i,j)=i \\ 53\text {if } p_i\le j\le l_i \ \ \text{adr}(i,j)=\min(m,n)+s_{i+1}+j-i\\ 54\text {if } q_j\le i\le c_j \ \ \text{adr}(i,j)=\min(m,n)+s_{m+1}+t_{j+1}+i-j\\ 55\end{array} 56\] 57\end{itemize} 58 59\subsection{The {\classtitle SkylineStorage} class} 60 61The {\class SkylineStorage} class has no member attribute. It has some basic constructors just calling {\class MatrixStorage} constructors: 62\vspace{.1cm} 63\begin{lstlisting}[] 64class SkylineStorage : public MatrixStorage 65{public : 66 SkylineStorage(AccessType = _dual); 67 SkylineStorage(Number, AccessType = _dual); 68 SkylineStorage(Number, Number, AccessType = _dual); 69 virtual ~SkylineStorage() {} 70\end{lstlisting} 71\vspace{.2cm} 72{\class SkylineStorage} objects do not have to be instantiated. This abstract class acts as an interface to particular skyline storages and gathers all common functionalities of skyline storages, some being virtuals (some template line declarations are omitted): 73\vspace{.1cm} 74\begin{lstlisting}[] 75virtual Number size() const = 0; 76virtual Number lowerPartSize() const = 0; 77virtual Number upperPartSize() const {return 0;} 78template<typename Iterator> 79 void printEntriesTriangularPart (StrucType, Iterator&, Iterator&, const std::vector<Number>&, Number, Number, Number, const String&, Number, std::ostream&) const; 80 void printCooTriangularPart(std::ostream&, Iterator&, const std::vector<Number>&, Number, Number, bool, SymType sym = _noSymmetry) const; 81 void printCooDiagonalPart(std::ostream&, Iterator&, Number) const; 82 83template <typename T> 84 void loadSkylineFromFileDense(std::istream&, std::vector<T>&, std::vector<Number>&, std::vector<Number>&, SymType, bool); 85 void loadSkylineFromFileCoo(std::istream&, std::vector<T>&, std::vector<Number>&, std::vector<Number>&, SymType, bool); 86 87template<typename MatIterator, typename VecIterator, typename ResIterator> 88 void diagonalMatrixVector(MatIterator&, VecIterator&, ResIterator&, ResIterator&) const; 89 void lowerMatrixVector(const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&, SymType sym) const; 90 void upperMatrixVector(const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&, SymType) const; 91 void diagonalVectorMatrix(MatIterator&, VecIterator&, ResIterator&, ResIterator&) const; 92 void lowerVectorMatrix(const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&, SymType sym) const; 93 void upperVectorMatrix(const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&, SymType) const; 94 void sumMatrixMatrix(Mat1Iterator&, Mat2Iterator&, ResIterator&, ResIterator&) const; 95\end{lstlisting} 96\vspace{.1cm} 97The {\class SkylineStorage} class provides the real code of product of matrix and vector, using iterators as arguments to by-pass the template type of matrix values. \\ 98In order to deal with some specific solvers, this class also comes with several following functions. For simplicity, we only declare template line once. 99\vspace{.1cm} 100\begin{lstlisting}[deletekeywords={[3] size\_t}] 101template<typename MatIterator, typename VecIterator, typename XIterator> 102void bzLowerSolver(const MatIterator&, const MatIterator&, const VecIterator&, const XIterator&, const XIterator&, const std::vector<size_t>::const_iterator) const; 103void bzLowerD1Solver(const MatIterator&, const VecIterator&, const XIterator&, const XIterator&, const std::vector<size_t>::const_iterator) const; 104void bzLowerConjD1Solver(const MatIterator&, const VecIterator&, const XIterator&, const XIterator&, const std::vector<size_t>::const_iterator) const; 105void bzDiagonalSolver(const MatIterator&, const VecIterator&, const XIterator&, const XIterator&) const; 106void bzUpperSolver(const MatRevIterator&, const MatRevIterator&, const VecRevIterator&, const XRevIterator&, const XRevIterator&, const std::vector<size_t>::const_reverse_iterator) const; 107void bzUpperD1Solver(const MatRevIterator&, const VecRevIterator&, const XRevIterator&, const XRevIterator&, const std::vector<size_t>::const_reverse_iterator) const; 108void bzUpperConjD1Solver(const MatRevIterator&, const VecRevIterator&, const XRevIterator&, const XRevIterator&, const std::vector<size_t>::const_reverse_iterator) const; 109\end{lstlisting} 110Like the product of matrix and vector, these functions take advantage of iterators as arguments to by-pass the template type of matrix values. 111\vspace{.2cm} 112\displayInfos{library=largeMatrix, header=SkylineStorage.hpp, implementation=SkylineStorage.cpp, test={test\_LargeMatrixSkylinetorage.cpp}, 113header dep={MatrixStorage.hpp, config.h, utils.h}} 114 115\subsection{{\classtitle DualSkylineStorage} and {\classtitle SymSkylineStorage} classes} 116 117Dual skyline storage travels the matrix as diagonal, strict lower and upper parts. Row pointers vectors are attached to strict lower part and column pointers vectors are attached to strict upper part: 118\vspace{.1cm} 119\begin{lstlisting} 120class DualSkylineStorage : public SkylineStorage 121{protected : 122 std::vector<Number> rowPointer_; // vector of positions of begining of rows 123 std::vector<Number> colPointer_; // vector of positions of begining of cols 124 ... 125\end{lstlisting} 126\vspace{.2cm} 127{\class DualSkylineStorage} class provides basic constructor (no construction of storage vectors), constructor from a pair of global numbering vectors or list of column indices by rows. 128Both of them create the compressed storage vector using the auxiliary function \verb?buildStorage?: 129\vspace{.1cm} 130\begin{lstlisting}[] 131DualSkylineStorage(Number nr = 0, Number nc = 0); 132DualSkylineStorage(Number, Number, const std::vector< std::vector<Number> >&, 133 const std::vector< std::vector<Number> >&); 134DualSkylineStorage(Number, Number, const std::vector< std::vector<Number> >&); 135~DualSkylineStorage() {}; 136\end{lstlisting} 137\vspace{.2cm} 138The class provides most of the virtual methods declared in {\class MatrixStorage} class. For sake of simplicity, we report here only function with \verb?Real? type but versions with \verb?Complex?, \verb?Matrix<Real>?, \verb?Matrix<Complex>? and mixed types are also provided. Some template line declarations are also omitted. 139\vspace{.1cm} 140\begin{lstlisting}[] 141Number size() const; 142Number lowerPartSize() const ; 143Number upperPartSize() const; 144Number pos(Number i, Number j, SymType s=_noSymmetry) const; 145void positions(const std::vector<Number>&, const std::vector<Number>&, 146 std::vector<Number>&, bool errorOn=true, 147 SymType s=_noSymmetry) const; 148void printEntries(std::ostream&, const std::vector<Real>&, Number, SymType) const; 149void printCooMatrix(std::ostream&, const std::vector<Real>&, 150 SymType s=_noSymmetry) const; 151void loadFromFileDense(std::istream&, std::vector<Real>&, SymType, bool ); 152void loadFromFileCoo(std::istream& , std::vector<Real>&, SymType, bool); 153template<typename M, typename V, typename R> 154 void multMatrixVector(const std::vector<M>&, const std::vector<V>&, 155 std::vector<R>&) const; 156 void multVectorMatrix(const std::vector<M>&, const std::vector<V>&, 157 std::vector<R>&) const; 158 void multMatrixVector(const std::vector<Real>&, const std::vector<Real>&, 159 std::vector<Real>&, SymType) const; 160 void multVectorMatrix(const std::vector<Real>&, const std::vector<Real>&, 161 std::vector<Real>&, SymType) const; 162template<typename M1, typename M2, typename R> 163 void addMatrixMatrix(const std::vector<M1>&, const std::vector<M2>&, std::vector<R>&) const; 164template<typename T> 165 void setDiagValueDualSkyline(std::vector<T>&, const T); 166template<typename M> 167void lu(std::vector<M>& m, std::vector<M>& fa) const; 168template<typename M, typename V, typename X> 169 void lowerD1Solver(const std::vector<M>&, std::vector<V>&, std::vector<X>&) const; 170 void diagonalSolver(const std::vector<M>&, std::vector<V>&, std::vector<X>&) const; 171 void upperD1Solver(const std::vector<M>&, std::vector<V>&, std::vector<X>&, const SymType) const; 172 void upperSolver(const std::vector<M>&, std::vector<V>&, std::vector<X>&) const; 173\end{lstlisting} 174\vspace{.2cm} 175 176The {\class SymSkylineStorage} class is very similar to the {\class DualSkylineStorage} class, except that it adresses only square matrix with symmetric storage and, thus, does not manage column pointers vectors because the upper part has same storage as lower part (with transposition) : 177\vspace{.1cm} 178\begin{lstlisting} 179class SymSkylineStorage : public SkylineStorage 180{protected : 181 std::vector<Number> rowPointer_; //!< vector of positions of begining of cols 182 ... 183\end{lstlisting} 184\vspace{.1cm} 185Note that the matrix may be not symmetric. In that case its upper part is stored.\\ 186Besides some similar functions to {\class DualSkylineStorage} class', the {\class SymSkylineStorage} provides more several functions to factorize matrix: LDLt and LDL* factorization 187\vspace{.1cm} 188\begin{lstlisting}[] 189template<typename M> 190 void ldlt(std::vector<M>& m, std::vector<M>& fa, const SymType sym = _symmetric) const; 191 void ldlstar(std::vector<M>& m, std::vector<M>& fa) const; 192\end{lstlisting} 193\vspace{.2cm} 194XXX=Dual or Sym 195 196\displayInfos{library=largeMatrix, header=XXXSkylineStorage.hpp, implementation=XXXSkylineStorage.cpp, test={test\_LargeMatrixSkylineStorage.cpp}, header dep={SkylineStorage.hpp, MatrixStorage.hpp, config.h, utils.h}} 197