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 LargeMatrix} class} 18 19The {\class LargeMatrix} class deals with large matrices with different type of storages. It is a template class \verb?LargeMatrix<T>? where \verb?T? is either the type \verb?Real?, \verb?Complex?, \verb?Matrix<Real>? or \verb?Matrix<Complex>?. Although it is possible to instantiate \verb?LargeMatrix<T>? with other types, most of functionalities do not work because some member functions are overloaded only for the previous types.\\ 20 21The class mainly manages a vector storing the matrix values and a pointer to a storage structure ({\class MatrixStorage} base class) : 22\vspace{.1cm} 23\begin{lstlisting}[deletekeywords={[3] name, valueType, strucType}] 24template <typename T> 25class LargeMatrix{ 26public : 27 ValueType valueType; // type of values (real, complex) 28 StrucType strucType; // struc of values (scalar, matrix) 29 Number nbRows; // number of rows counted in T 30 Number nbCols; // number of columns counted in T 31 SymType sym; // type of symmetry 32 Dimen nbRowsSub; // number of rows of submatrices (1 if scalar values) 33 Dimen nbColsSub; // number of columns of submatrices (1 if scalar values) 34 String name; // optionnal name, useful for documentation 35protected : 36 vector<T> values_; // list of values ordered along storage method 37 MatrixStorage* storage_p; // pointer to the matrix storage 38 ... 39\end{lstlisting} 40\vspace{.1cm} 41When \verb?T=Matrix<Real>? or \verb?T=Matrix<Complex>?, the number of rows and columns are counted in \verb?Matrix? not in scalar! \\ \\ 42\textdbend \ \ The first value (index 0) of the vector \verb?values_? is not used for storing a matrix value. It contains the default value (in general \verb?T()?) returned when trying to access to a matrix value outside the storage. Be care with default Matrix<T>(), its a $0 \times 0$ matrix!\\ 43\\ 44There are various constructors : three from a existing storage, some from file where the large matrix is loaded from the file and some from matrix dimensions and a default value. A large matrix which is created by the copy constructor (shallow copy) will share a same storage with the input one. Be care with these last constructors, they allocate the vector of values only for dense storage. 45 46\vspace{.1cm} 47\begin{lstlisting}[deletekeywords={[3] dims}] 48LargeMatrix(); 49LargeMatrix(const LargeMatrix<T>&); 50LargeMatrix(MatrixStorage*, const T&= T(), SymType= _noSymmetry); 51LargeMatrix(MatrixStorage* ms, dimPair dims, SymType sy = _noSymmetry); 52LargeMatrix(MatrixStorage*, SymType); 53LargeMatrix(Number, Number, StorageType = _dense, AccessType = _row, const T& = T()); 54LargeMatrix(Number, Number, StorageType = _dense, SymType = _symmetric, const T& = T()); 55LargeMatrix(const String&, StorageType, Number , Number, StorageType = _dense, 56 AccessType = _row, Number nsr=1, Number nsc=1); 57LargeMatrix(const String&, StorageType, Number, StorageType = _dense, 58 SymType = _symmetric, Number nsr=1); 59 ~LargeMatrix(); 60 61void init(MatrixStorage*, const T&, SymType); 62void allocateStorage(StorageType, AccessType, const T& = T()); 63void setType(const T&); 64\end{lstlisting} 65\vspace{.1cm} 66The \verb?StorageType?, \verb?AccessType? and \verb?SymType? are enumerations : 67\vspace{.1cm} 68\begin{lstlisting}[] 69enum StorageType{_noStorage=0,_dense,_cs,_skyline,_coo}; 70enum AccessType {_noAccess=0,_sym,_row,_col,_dual}; 71enum SymType {_noSymmetry=0,_symmetric,_skewSymmetric,_selfAdjoint,_skewAdjoint, 72 _diagonal}; 73\end{lstlisting} 74\vspace{.1cm} 75The storage type \verb?_cs? means compressed sparse. It refers to the well known CSR/CSC storage (see the next section for details on storage). The storage type \verb?_coo? means coordinates storage ($i,j,a_{ij}$). So far, it is not proposed as a storage method for large matrix but large matrix may be saved to file in this format. Constructors use the utilities: init, allocateStorage and setType. The SetType function uses some utilities (external functions) to get the dimensions of sub-matrices: 76\vspace{.1cm} 77\begin{lstlisting}[] 78pair<Dimen,Dimen> dimsOf(const Real&); 79pair<Dimen,Dimen> dimsOf(const Complex&); 80pair<Dimen,Dimen> dimsOf(const Matrix<Real>&); 81pair<Dimen,Dimen> dimsOf(const Matrix<Complex>&); 82bool isDiagonal() const; 83bool isId(const double & tol=0.) const; 84\end{lstlisting} 85\vspace{.1cm} 86\textdbend \ \ Note that for the moment, it is not possible to load from a file a matrix as a matrix of matrices (use only \verb?nsr=1? and \verb?nsc=1? arguments in constructors from file).\\ 87\\ 88There are few facilities to access to matrix values : 89\vspace{.1cm} 90\begin{lstlisting}[] 91Number pos(Number i, Number j) const; 92void positions(const vector<number_t>&, const vector<number_t>&, 93 vector<number_t>&, bool = false) const; 94vector<pair<number_t, number_t>> getCol(number_t, number_t =1, number_t=0) const; 95vector<pair<number_t, number_t>> getRow(number_t, number_t =1, number_t =0) const; 96T operator()(number_t adr) const; 97T& operator()(number_t adr); 98typename vector<T>::iterator at(number_t, number_t); 99typename vector<T>::const_iterator at(number_t i, number_t j) const; 100T operator() (Number i, Number j) const; 101T& operator()(Number i, Number j, bool =true); 102Number nbNonZero() const; 103void resetZero(); 104\end{lstlisting} 105\vspace{.1cm} 106The member function \verb?pos(i,j)? returns the adress in vector \verb?values_? of the matrix value $M_{ij}$ for $i,j\ge 1$. The returned adress begins at 1 except when the value is outside the storage where the function returns 0. The access operator \verb?(i,j)? returns the matrix value or a reference to it with the {\em non const} version. In that case, it is possible to activate (\verb?errorOn=true?) an error handler to prevent access to values outside the storage. It is the default behaviour. When \verb?errorOn=false?, the access operator returns a reference to \verb?values_[0]?, thus you can modify its value. This trick avoids to perform test like \verb?if(pos(i,j)!=0)...? before writing matrix values. To be safe at end, reset the \verb?values_[0]? to 0 with {\cmd resetZero} function. Note that the skyline or compressed storage structure can not be dynamically changed when "inserting" a new value. \\ \\ 107\textdbend \ \ \verb?pos(i,j)?, \verb?getRow? and \verb?getCol? function should not be used in heavy computations because it is time consuming. In heavy computation, write codes based on storage structure! \\ 108\\ 109As {\class LargeMatrix} class provides a lot of functions to modify its contents, some being time expansive regarding the storage : 110\vspace{.1cm} 111\begin{lstlisting}[] 112void deleteRows(number_t, number_t); 113void deleteCols(number_t, number_t); 114void setColToZero(number_t c1=0, number_t c2=0); 115void setRowToZero(number_t r1=0, number_t r2=0); 116void roundToZero(real_t aszero=10*theEpsilon); 117 118LargeMatrix<T>& operator*=(const K&); 119LargeMatrix<T>& operator/=(const K&); 120LargeMatrix<T>& operator+=(LargeMatrix<T>&); 121LargeMatrix<T>& operator-=(LargeMatrix<T>&); 122 123LargeMatrix<T>& add(iterator&,iterator&); 124LargeMatrix<T>& add(const vector<T>&, const vector<number_t>&,const vector<number_t>&); 125LargeMatrix<T>& add(const T& c, const vector<number_t>&, const vector<number_t>&); 126LargeMatrix<T>& add(const LargeMatrix<K>&, const vector<number_t>&,const vector<number_t>&, C); 127void addColToCol(number_t, number_t, complex_t a, bool =false); 128void addRowToRow(number_t, number_t, complex_t a, bool =false); 129void copyVal(const LargeMatrix<T>&, const vector<number_t>&, const vector<number_t>&); 130LargeMatrix<T>& assign(const LargeMatrix<K>&, const vector<number_t>&, const vector<number_t>&); 131\end{lstlisting} 132\vspace{.2cm} 133Some functions are provided to change the scalar/matrix representation of the matrix or the storage: 134\vspace{.1cm} 135\begin{lstlisting}[] 136void toSkyline(); 137void toStorage(MatrixStorage*); 138LargeMatrix<K>* toScalar(K); 139void toScalarEntries(const vector<Matrix<K> >&, vector<K>&, const MatrixStorage&); 140void toUnsymmetric(AccessType at=_sym); 141void extract2UmfPack(vector<int_t>& colPointer,vector<int_t>& rowIndex,vector<T>& matA) const; 142bool getCsColUmfPack(vector<int_t>& colPointer,vector<int_t>& rowIndex, const T*& matA) const; 143\end{lstlisting} 144\vspace{.2cm} 145Matrices that are stored as Dense or Skyline can be factorized as LU, LDLt, LDL* and matrices that are stored as Skyline or Sparse can be factorized using UMFpack if it is available. All the stuff related to solve linear system is also available. 146\vspace{.1cm} 147\begin{lstlisting}[] 148void ldltFactorize(); 149void ldlstarFactorize(); 150void luFactorize(bool withPermutation=true); 151void umfpackFactorize(); 152 153void ldltSolve(vector<S1>& vec, vector<S2>& res) const; 154void ldlstarSolve(vector<S>& vec, vector<T>& res) const; 155void luSolve(vector<S1>& vec, vector<S2>& res) const; 156void umfluSolve(vector<S1>& vec, vector<S2>& res) const; 157 158void sorLowerMatrixVector(const vector<S1>& vec, vector<S2>& res, real_t) const; 159void sorDiagonalMatrixVector(const vector<S1>& vec, vector<S2>& res,real_t) const; 160void sorUpperMatrixVector(const vector<S1>& vec, vector<S2>& res,real_t) const; 161void sorLowerSolve(const vector<S1>& vec, vector<S2>& res,real_t) const; 162void sorDiagonalSolve(const vector<S1>& vec, vector<S2>& res,real_t) const; 163void sorUpperSolve(const vector<S1>& vec, vector<S2>& res,real_t ) const; 164 165real_t umfpackSolve(const vector<S>& vec, 166vector<typename Conditional<NumTraits<ScalarType>::IsComplex, ScalarType, S>::type>& res, 167bool soa=true); 168real_t umfpackSolve(const vector<vector<S>*>&, 169vector<vector<typename Conditional<NumTraits<ScalarType>::IsComplex, ScalarType,S>::type>* >&, 170bool soa=true); 171real_t umfpackSolve(const vector<int_t>& colPointer, const vector<int_t>& rowIndex, 172const vector<T>& values, const vector<S>& vec, 173vector<typename Conditional<NumTraits<ScalarType>::IsComplex, ScalarType, S>::type>& res, 174bool soa=true); 175\end{lstlisting} 176\vspace{.2cm} 177In a same way, some functions interfaces eigen solvers: 178\vspace{.1cm} 179\begin{lstlisting}[] 180 181friend number_t eigenDavidsonSolve(const LargeMatrix<K>* pA, const LargeMatrix<K>* pB, 182 vector<pair<complex_t,Vector<complex_t> > >& res, number_t nev, real_t tol, 183 string_t which, bool isInverted, FactorizationType fac, bool isShift); 184friend number_t eigenKrylovSchurSolve(const LargeMatrix<K>* pA, const LargeMatrix<K>* pB, 185 vector<pair<complex_t,Vector<complex_t> > >& res,number_t nev, real_t tol, 186 string_t which, bool isInverted, FactorizationType fac, bool isShift); 187 188\end{lstlisting} 189\vspace{.2cm} 190The Frobenius norm $\sqrt{\sum{|a_{ij}|^2}}$ and infinity norm $\max_{|a_{ij}|}$ are provided: 191\vspace{.1cm} 192\begin{lstlisting}[] 193real_t norm2() const; 194real_t norminfty() const; 195real_t partialNormOfCol(number_t, number_t, number_t) const; 196\end{lstlisting} 197\vspace{.2cm} 198Some operations are available as external functions to the class. \\ 199 200The {\class LargeMatrix} supplies methods to add two large matrices, which share a same storage. The matrix result will use the same storage as the two added matrices. If the result doesn't point to the storage of addend, after the addition, it will be forced to use this storage, and the number of object sharing this storage increases by one. Similar to multiplication of matrix-vector, these functions are also external templates, with specializations to mixed real and complex types. 201\vspace{.1cm} 202\begin{lstlisting}[] 203friend void addMatrixMatrix(const LargeMatrix<S>& matA, const LargeMatrix<S>& matB, LargeMatrix<S>& matC); 204friend void addMatrixMatrix(const LargeMatrix<Complex>& matA, const LargeMatrix<Real>& matB, LargeMatrix<Complex>& matC); 205friend void addMatrixMatrix(const LargeMatrix<Real>& matA, const LargeMatrix<Complex>& matB, LargeMatrix<Complex>& matC); 206\end{lstlisting} 207\vspace{.1cm} 208The operator \verb?+? is overloaded for {\class LargeMatrix} in the same manner. 209\vspace{.1cm} 210\begin{lstlisting}[] 211LargeMatrix<T> operator+(const LargeMatrix<T>& matA, const LargeMatrix<T>& matB); 212LargeMatrix<Complex> operator+(const LargeMatrix<Real>& matA, const LargeMatrix<Complex>& matB); 213LargeMatrix<Complex> operator+(const LargeMatrix<Complex>& matA, const LargeMatrix<Real>& matB); 214\end{lstlisting} 215\vspace{.1cm} 216\textdbend \ \ The symmetry of the result matrix depends on the ones of both addends. If one of these has \verb?_noSymmetry? as \verb?SymType?, the result will have the same \verb?SymType? (i.e: \verb?_noSymmetry?). \\ 217A large matrix can multiply by a scalar with the external templated functions 218\vspace{.1cm} 219\begin{lstlisting}[] 220friend LargeMatrix<S> multMatrixScalar(const LargeMatrix<S>&, const S); 221friend LargeMatrix<Complex> multMatrixScalar(const LargeMatrix<Complex>&, const Real); 222friend LargeMatrix<Complex> multMatrixScalar(const LargeMatrix<Real>&, const Complex); 223\end{lstlisting} 224\vspace{.1cm} 225Like {\cmd addMatrixMatrix}, the result of {\cmd multMatrixScalar} shares the same storage of the input large matrix.\\ 226The operator \verb?*? can be also overloaded in the same manner. 227\vspace{.1cm} 228\begin{lstlisting}[] 229LargeMatrix<T> operator*(const LargeMatrix<T>& mat, const T v); 230LargeMatrix<T> operator*(const T v, const LargeMatrix<T>& mat); 231LargeMatrix<Complex> operator*(const LargeMatrix<Complex>& mat, const Real v); 232LargeMatrix<Complex> operator*(const Real v, const LargeMatrix<Complex>& mat); 233LargeMatrix<Complex> operator*(const LargeMatrix<Real>& mat, const Complex v); 234LargeMatrix<Complex> operator*(const Complex v, const LargeMatrix<Real>& mat); 235\end{lstlisting} 236\vspace{.2cm} 237The product of matrix and vector is one of the most important operation required. The {\class LargeMatrix} class interfaces matrix-vector product and the vector-matrix matrix; the product being really done by the storage classes. They are external templated functions to the class (declared friend of class), with specializations to mixed real and complex types (only casting from real to complex is allowed) : 238\vspace{.1cm} 239\begin{lstlisting}[] 240void multMatrixVector(const LargeMatrix<S>&, const vector<S>&, vector<S>&); 241void multVectorMatrix(const LargeMatrix<S>&, const vector<S>&, vector<S>&); 242void multMatrixVector(const LargeMatrix<Matrix<S> >&, 243 const vector<Vector<S> >&, vector<Vector<S> >&); 244void multVectorMatrix(const LargeMatrix<Matrix<S> >&, 245 const vector<Vector<S> >&, vector<Vector<S> >&); 246\end{lstlisting} 247\vspace{.1cm} 248The operator \verb?*? is overloaded for {\class LargeMatrix} and vector or vector and {\class LargeMatrix} in a same way : 249\goodbreak 250\vspace{.1cm} 251\begin{lstlisting}[] 252friend vector<S> operator*(const LargeMatrix<S>&, const vector<S>&); 253friend vector<S> operator*(const vector<S>&, const LargeMatrix<S>&); 254friend vector<Vector<S> operator*(const LargeMatrix<Matrix<S> >&, 255 const vector<Vector<S> >&); 256friend vector<Vector<S> operator*(const vector<Vector<S> >&, 257 const LargeMatrix<Matrix<S> >&); 258... 259\end{lstlisting} 260\vspace{.1cm} 261The product of two matrices is provided, but be cautious, up to now the result matrix is a dense matrix (except if one matrix is a diagonal one): 262\vspace{.1cm} 263\begin{lstlisting}[] 264void multMatrixMatrix(const LargeMatrix<SA>&, const LargeMatrix<SB>&, LargeMatrix<SR>&); 265\end{lstlisting} 266\vspace{.1cm} 267 268The {\class LargeMatrix} also provides some functions to factorize a matrix and solve factorized syste. Like {\cmd multMatrixVector}, these functions are external (friend of class), template with hybrid complex-real specialization. 269\vspace{.1cm} 270\begin{lstlisting}[] 271void ldltFactorize(LargeMatrix<S>& mat); 272void ldlstarFactorize(LargeMatrix<S>& mat); 273void luFactorize(LargeMatrix<S>& mat); 274\end{lstlisting} 275\vspace{.1cm} 276\begin{remark} 277These functions are available for {\class SkylineStorage} and {\class DenseStorage}. 278\end{remark} 279If a matrix is factorized, some other operation are available 280\vspace{.1cm} 281\begin{lstlisting}[] 282void multFactMatrixVector(const LargeMatrix<S1>&, const vector<S2>&, vector<S3>&); 283void multVectorFactMatrix(const LargeMatrix<S1>&, const vector<S2>&, vector<S3>&); 284void multInverMatrixVector(const LargeMatrix<S1>& mat, const vector<S2>& vec, 285 vector<typename Conditional<NumTraits<S1>::IsComplex, S1, S2>::type>& res, 286 FactorizationType); 287\end{lstlisting} 288\vspace{.1cm} 289 290In order to generate a diagonal matrix which uses an existent matrix as prototype, {\class LargeMatrix} class provides function 291\vspace{.1cm} 292\begin{lstlisting}[] 293LargeMatrix<S> diagonalMatrix(const LargeMatrix<S>&, const S); 294\end{lstlisting} 295\vspace{.1cm} 296A special version of this function is a function to generate an identity matrix from an existing one. 297\vspace{.1cm} 298\begin{lstlisting}[] 299LargeMatrix<S> identityMatrix(const LargeMatrix<S>&); 300\end{lstlisting} 301\vspace{.1cm} 302\textdbend \ \ These functions work with {\class LargeMatrix} class having all kinds of matrix storage except {\class RowCsStorage} class and {\class ColCsStorage} class. \\ 303Finally, the class has input/output member functions : 304\vspace{.1cm} 305\begin{lstlisting}[] 306void print(ostream&) const; 307void viewStorage(ostream&) const; 308ostream& operator<<(ostream&, const LargeMatrix<S>&); 309void saveToFile(const String &, StorageType, bool encode=false) const; 310void loadFromFile(const String&, StorageType); 311String encodeFileName(const String&, StorageType) const; 312\end{lstlisting} 313\vspace{.1cm} 314Only two formats are allowed for saving matrix to file or loading matrix from file : 315\begin{itemize} 316\item the dense format (\verb?_dense?), where all the matrix values are written row by row (one row by line), space character separates each value : 317\begin{verbatim} 318A11 A12 ... A1n re(A11) im(A11) ... re(A1n) im(A1n) 319A21 A22 ... A2n if complex values re(A21) im(A21) ... re(A2n) im(A2n) 320... ... 321Am1 Am2 ... Amn re(Am1) im(Am1) ... re(Amn) im(Amn) 322\end{verbatim} 323\item the coordinates format (\verb?_coo?) 324\begin{verbatim} 325i1 j1 Ai1j1 i1 j1 re(Ai1j1) im(Ai1j1) 326i2 j2 Ai2j2 if complex values i2 j2 re(Ai2j2) im(Ai2j2) 327... 328\end{verbatim} 329\end{itemize} 330\textdbend \ \ For matrix of matrices (\verb?T=Matrix<Real>? or \verb?T=Matrix<Complex>?), submatrix structure is not preserved when writing matrix values!\\ \\ 331To keep some general informations about matrix, with the argument \verb?encode=true? in the \verb?saveFile? function, the name of file may be modified by inserting the string : 332\begin{center} 333\verb?(m_n_storageType_valueType_strucType_p_q)? 334\end{center} 335where \verb?m? and \verb?n? are the "dimensions" of the matrix, \verb?storageType = dense or coo?, \verb?valueType = real or complex? , \verb?strucType = scalar or matrix? and \verb?p? and \verb?q? are the dimensions of sub-matrices. In case of scalar value, this parameters are omitted. This trick avoids to include informations in file in order that it is easily loadable by other softwares, in particular Matlab. 336 337\subsection*{Some examples} 338Manipulate a $3\times 2$ matrix with row dense storage: 339\vspace{.1cm} 340\begin{lstlisting}[] 341Vector<Real> x(2,2); 342LargeMatrix<Real> A(3,2,_dense,_row,1.); 343out<<"matrix A : "<<A; 344out<<"access A(1,1)= "<<A(1,1); 345out<<"product A*x : "<<A*x; 346A.saveToFile("A_dense.mat",_dense); 347LargeMatrix<Real> Areload("A_dense.mat",_dense,3,2,_dense,_row); 348\end{lstlisting} 349\vspace{.2cm} 350Manipulate a $3\times 2$ matrix of $2\times 2$ real matrices with dual dense storage: 351\vspace{.1cm} 352\begin{lstlisting}[] 353Matrix<Real> I1(2,_idMatrix); 354Vector<Real> x(6,3); 355LargeMatrix<Matrix<Real> > B(6,3,_dense,_dual,2*I1); 356out<<"product X*B : "<<X*B; 357\end{lstlisting} 358\vspace{.2cm} 359Manipulate a $9\times 9$ symmetric matrix with symmetric compressed storage: 360\vspace{.1cm} 361\begin{lstlisting}[deletekeywords={[3] x}] 362//construct storage for a Q1 mesh 363Number n=3; 364vector<vector<Number> > elts((n-1)*(n-1),vector<Number>(4)); 365for(Number j=1; j<=n-1; j++) 366 for(Number i=1; i<=n-1; i++) { 367 Number e=(j-1)*(n-1)+i-1, p=(j-1)*n+i; 368 elts[e][0]=p; elts[e][1]=p+1; elts[e][2]=p+n; elts[e][3]=p+n+1; 369 } 370SymCsStorage* cssym= new SymCsStorage(n*n,elts,elts); 371//construct large matrix with symmetric compressed storage 372LargeMatrix<Real> C(cssym,1.,_symmetric); 373C.viewStorage(out); 374Vector<Real> x(n*n,0.);for(Number i=0;i<n*n;i++) x1[i]=i; 375out<<"product C*x : "<<C*x; 376C.saveToFile("C_coo.mat"),_coo); 377//construct large matrix with symmetric compressed storage 378LargeMatrix<Real> C1(cssym,1.,_symmetric); 379out<<"addition C+C1 "<< C+C1; 380// multiply a matrix with a scalar 381out<<"product C*10.0 "<< C*10.0; 382// generate a diagonal matrix 383out<<"diagonal matrix from matrix C " << diagonalMatrix(C, 10.0); 384// generate an identity matrix 385out<<"identity matrix from matrix C " << identityMatrix(C); 386\end{lstlisting} 387\vspace{.1cm} 388Note that compressed storage has to be built before to instantiate a large matrix with compressed storage. The function {\cmd viewStorage} produced the output: 389\begin{verbatim} 390symmetric_compressed sparse (csr,csc), size = 29, 9 rows, 9 columns 391 (shared by 1 objects). 392 123456789 393 1 d........ 394 2 xd....... 395 3 .xd...... 396 4 xx.d..... 397 5 xxxxd.... 398 6 .xx.xd... 399 7 ...xx.d.. 400 8 ...xxxxd. 401 9 ....xx.xd 402 123456789 403\end{verbatim} 404\vspace{.1cm} 405LDLt-Factorize a positive definite 3x3 matrix with symmetric skyline storage: 406\vspace{.1cm} 407\begin{lstlisting}[] 408LargeMatrix<Real> rResSymSkylineLdlt(inputsPathTo("matrix3x3SymPosDef.data"), _dense, rowNum, _skyline, _symmetric); 409ldltFactorize(rResSymSkylineLdlt); 410out << "The result of LDLt factorizing sym skyline matrix is " << rResSymSkylineLdlt << endl; 411\end{lstlisting} 412\vspace{.1cm} 413\displayInfos{library=largeMatrix, header=LargeMatrix.hpp, implementation=LargeMatrix.cpp, test={test\_LargeMatrixDenseStorage.cpp, test\_LargeMatrixCsStorage.cpp}, header dep={MatrixStorage.hpp, DenseStorage.hpp, CsStorage.hpp, SkylineStorage.hpp, config.h, utils.h}} 414