Matrix storages

The largeMatrix library deals with large matrices stored in different ways: dense storages (all the matrix values are stored), skyline storages (values on a row or a column are stored from the first nonzero values) and compressed sparse storage (only nonzero values are stored). For each of these storages, different “access” methods, driving how the matrix values are stored, are proposed: by rows, by columns, by rows and columns (say dual access) and symmetric access which is like dual access but with same row and column access. The symmetric access is well suited for matrices having symmetry property (symmetric, skew-symmetric, self-adjoint or skew-adjoint), or matrix with symmetric storage property, currently met in finite element approximation. Note that the row dense storage corresponds to the storage used by the Matrix class defined in The Matrix class. As other access methods may be useful in finite element methods, it is the reason why the largeMatrix library provides dense storages.

Large matrices under consideration may be real or complex value matrices, but also matrices of real/complex matrices. This last case allows dealing with vector variational problems where the unknown and test function are vector functions, using the same storage structure as in the scalar case.

More precisely, the largeMatrix library is based on the template class LargeMatrix<T> which mainly handles a vector to store the matrix values and a pointer to a storage describing how the matrix values are stored. Different storage structures are organized in a hierarchy of classes, terminal classes carrying specific data and methods of storages, in particular the product of a matrix and a vector.

The LargeMatrix class

The LargeMatrix<T> class deals with large matrices with different types of storages. Its template parameter T is one of the following types: Real, Complex, RealMatrix or ComplexMatrix. Although it is possible to instantiate LargeMatrix<T> with other types, most functionalities do not work because some member functions are overloaded only for the previous types.

The class mainly manages a vector storing the matrix values and a pointer to a storage structure (MatrixStorage base class):

template <typename T>
class LargeMatrix
{
  public :
    ValueType valueType;      // type of values (real, complex)
    StrucType strucType;      // struc of values (scalar, matrix)
    Number nbRows;            // number of rows counted in T
    Number nbCols;            // number of columns counted in T
    SymType sym;              // type of symmetry
    Dimen nbRowsSub;          // number of rows of submatrices (1 if scalar values)
    Dimen nbColsSub;          // number of columns of submatrices (1 if scalar values)
    String name;              // optionnal name, useful for documentation
  protected :
    vector<T> values_;   // list of values ordered along storage method
    MatrixStorage* storage_p; // pointer to the matrix storage
    ....
}

When T=Matrix<Real> or T=Matrix<Complex>, the number of rows and columns are counted in Matrix not in scalar!

Hint

The first value (index 0) of the vector values_ is not used for storing a matrix value. It contains the default value (in general T) returned when trying to access to a matrix value outside the storage. Be care with default Matrix<T>(), it is a \(0 \times 0\) matrix!

There are various constructors : three from an 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.

LargeMatrix();
LargeMatrix(const LargeMatrix<T>&);
LargeMatrix(MatrixStorage*, const T&= T(), SymType= _noSymmetry);
LargeMatrix(MatrixStorage* ms, dimPair dims, SymType sy = _noSymmetry);
LargeMatrix(MatrixStorage*, SymType);
LargeMatrix(Number, Number, StorageType = _dense, AccessType = _row, const T& = T());
LargeMatrix(Number, Number, StorageType = _dense, SymType = _symmetric, const T& = T());
LargeMatrix(const String&, StorageType, Number , Number, StorageType = _dense, AccessType = _row,  Number nsr=1, Number nsc=1);
LargeMatrix(const String&, StorageType, Number, StorageType = _dense, SymType = _symmetric,  Number nsr=1);
~LargeMatrix();

void init(MatrixStorage*, const T&, SymType);
void allocateStorage(StorageType, AccessType, const T& = T());
void setType(const T&);

The StorageType , AccessType and SymType are enumerations :

enum StorageType { _noStorage =0 , _dense , _cs , _skyline , _coo } ;
enum AccessType { _noAccess =0 ,_sym , _row , _col , _dual } ;
enum SymType { _noSymmetry=0 , _symmetric , _skewSymmetric , _ selfAdjoint , _skewAdjoint , _diagonal } ;

The storage type _cs means compressed sparse. It refers to the well known CSR/CSC storage (see the next section for details on storage). The storage type _coo means coordinates storage \((i , j , a i j )\). 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:

pair <Dimen, Dimen> dimsOf ( const Real&) ;
pair <Dimen, Dimen> dimsOf ( const Complex&) ;
pair <Dimen, Dimen> dimsOf ( const Matrix<Real >&) ;
pair <Dimen, Dimen> dimsOf ( const Matrix<Complex>&) ;
bool isDiagonal ( ) const ;
bool i s I d ( const double & tol = 0 . ) const ;

Hint

Note that for the moment, it is not possible to load from a file a matrix as a matrix of matrices (use only nsr=1 and nsc=1 arguments in constructors from file).

There are few facilities to access to matrix values :

Number pos (Number i , Number j ) const ;
void positions ( const vector <Number>&, const vector <Number>&,
vector <Number>&, bool = f a l s e ) const ;
vector <pair <Number , Number>> getCol ( Number , Number =1 , Number=0) const ;
vector <pair <Number , Number>> getRow ( Number , Number =1 , Number =0) const ;
T operator() ( Number adr ) const;
T& operator() ( Number adr );
typename vector <T>::iteratorat( Number , Number );
typename vector <T>::const_iterator at (Number i, Number j ) const;
T operator() (Number i , Number j ) const;
T& operator() (Number i , Number j , bool =true );
Number nbNonZero() const;
void resetZero( );

The member function pos(i,j) returns the address in vector values_ of the matrix value \(M_{ij}\) for \(i,j\ge 1\). The returned address begins at 1 except when the value is outside the storage where the function returns 0. The access operator (i,j) returns the matrix value or a reference to it with the non-const version. In that case, it is possible to activate (errorOn=true) an error handler to prevent access to values outside the storage. It is the default behavior. When errorOn=false, the access operator returns a reference to values_[0], thus you can modify its value. This trick avoids performing tests like if(pos(i,j)!=0)... before writing matrix values. To be safe at end, reset the values_[0] to 0 with resetZero function. Note that the skyline or compressed storage structure can not be dynamically changed when “inserting” a new value.

Hint

pos(), getRow() and getCol() function should not be used in heavy computations because it is time-consuming. In heavy computation, write codes based on storage structure!

As LargeMatrix class provides a lot of functions to modify its contents, some being time expansive regarding the storage:

void deleteRows ( Number , Number ) ;
void deleteCols ( Number , Number ) ;
void setColToZero ( Number c1 =0 , Number c2 =0) ;
void setRowToZero ( Number r1 =0 , Number r2 =0) ;
void roundToZero ( r e a l _ t aszero =10 * theEpsilon ) ;

LargeMatrix<T>& operator * =( const K&) ;
LargeMatrix<T>& operator /=( const K&) ;
LargeMatrix<T>& operator +=( LargeMatrix<T>&) ;
LargeMatrix<T>& operator -=( LargeMatrix<T>&) ;
LargeMatrix<T>& add( i t e r a t o r &, i t e r a t o r &);
LargeMatrix<T>& add( const vector <T>&, const vector <Number>&,const vector <Number>&);
LargeMatrix<T>& add( const T& c , const vector <Number>&, const vector <Number>&);
LargeMatrix<T>& add( const LargeMatrix<K>&, const vector <Number>&,const vector <Number>&, C);
void addColToCol(Number ,Number ,Complex a, bool=false );
void addRowToRow(Number ,Number ,Complex a, bool=false );
void copyVal( const LargeMatrix<T>&, const vector <Number>&, const vector <Number>&);
LargeMatrix<T>& assign ( const LargeMatrix<K>&, const vector <Number>&, const vector <Number>&);

Some functions are provided to change the scalar/matrix representation of the matrix or the storage:

void toSkyline ( ) ;
void toStorage ( MatrixStorage * ) ;
LargeMatrix<K> * toScalar (K) ;
void toScalarEntries ( const vector <Matrix<K> >&, vector <K>&, const MatrixStorage&) ;
void toUnsymmetric ( AccessType at=_sym) ;
void extract2UmfPack ( vector <Int >& colPointer , vector <Int >& rowIndex , vector <T>& matA) const ;
bool getCsColUmfPack ( vector <Int >& colPointer , vector <Int >& rowIndex , const T*& matA) const ;

Matrices 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.

void ldltFactorize();
void ldlstarFactorize();
void luFactorize(bool withPermutation=true);
void umfpackFactorize();

void ldltSolve(vector<S1>& vec, vector<S2>& res) const;
void ldlstarSolve(vector<S>& vec, vector<T>& res) const;
void luSolve(vector<S1>& vec, vector<S2>& res) const;
void umfluSolve(vector<S1>& vec, vector<S2>& res) const;

void sorLowerMatrixVector(const vector<S1>& vec, vector<S2>& res, Real) const;
void sorDiagonalMatrixVector(const vector<S1>& vec, vector<S2>& res,Real) const;
void sorUpperMatrixVector(const vector<S1>& vec, vector<S2>& res,Real) const;
void sorLowerSolve(const vector<S1>& vec, vector<S2>& res,Real) const;
void sorDiagonalSolve(const vector<S1>& vec, vector<S2>& res,Real) const;
void sorUpperSolve(const vector<S1>& vec, vector<S2>& res,Real ) const;

Real umfpackSolve(const vector<S>& vec,
vector<typename Conditional<NumTraits<ScalarType>::IsComplex, ScalarType, S>::type>& res, bool soa=true);
Real umfpackSolve(const vector<vector<S>*>&, vector<vector<typename Conditional<NumTraits<ScalarType>::IsComplex, ScalarType,S>::type>* >&, bool soa=true);
Real umfpackSolve(const vector<Int>& colPointer, const vector<Int>& rowIndex, const vector<T>& values, const vector<S>& vec, vector<typename Conditional<NumTraits<ScalarType>::IsComplex, ScalarType, S>::type>& res,
bool soa=true);

In a same way, some functions interfaces eigensolvers:

friend Number eigenDavidsonSolve(const LargeMatrix<K>* pA, const LargeMatrix<K>* pB, vector<pair<Complex, Vector<Complex> > >& res, Number nev, Real tol, String which, bool isInverted, FactorizationType fac, bool isShift);
friend Number eigenKrylovSchurSolve(const LargeMatrix<K>* pA, const LargeMatrix<K>* pB, vector<pair<Complex, Vector<Complex> > >& res, Number nev, Real tol, String which, bool isInverted, FactorizationType fac, bool isShift);

The Frobenius norm \(\sqrt{\sum{|a_{ij}|^2}}\) and infinity norm \(\max_{|a_{ij}|}\) are provided:

Real norm2() const;
Real norminfty() const;
Real partialNormOfCol(Number, Number, Number) const;

Some operations are available as external functions to the class.

The 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.

friend void addMatrixMatrix(const LargeMatrix<S>& matA, const LargeMatrix<S>& matB, LargeMatrix<S>& matC);
friend void addMatrixMatrix(const LargeMatrix<Complex>& matA, const LargeMatrix<Real>& matB, LargeMatrix<Complex>& matC);
friend void addMatrixMatrix(const LargeMatrix<Real>& matA, const LargeMatrix<Complex>& matB, LargeMatrix<Complex>& matC);

The operator + is overloaded for LargeMatrix in the same manner.

LargeMatrix<T> operator+(const LargeMatrix<T>& matA, const LargeMatrix<T>& matB);
LargeMatrix<Complex> operator+(const LargeMatrix<Real>& matA, const LargeMatrix<Complex>& matB);
LargeMatrix<Complex> operator+(const LargeMatrix<Complex>& matA, const LargeMatrix<Real>& matB);

Like addMatrixMatrix, the result of multMatrixScalar shares the same storage of the input large matrix.

The operator * can be also overloaded in the same manner.

LargeMatrix<T> operator*(const LargeMatrix<T>& mat, const T v);
LargeMatrix<T> operator*(const T v, const LargeMatrix<T>& mat);
LargeMatrix<Complex> operator*(const LargeMatrix<Complex>& mat, const Real v);
LargeMatrix<Complex> operator*(const Real v, const LargeMatrix<Complex>& mat);
LargeMatrix<Complex> operator*(const LargeMatrix<Real>& mat, const Complex v);
LargeMatrix<Complex> operator*(const Complex v, const LargeMatrix<Real>& mat);

The product of matrix and vector is one of the most important operation required. The LargeMatrix class interfaces matrix-vector product and the vector-matrix matrix; the product being really done by the storage classes. They are external template functions to the class (declared friend of class), with specializations to mixed real and complex types (only casting from real to complex is allowed) :

void multMatrixVector(const LargeMatrix<S>&, const vector<S>&, vector<S>&);
void multVectorMatrix(const LargeMatrix<S>&, const vector<S>&, vector<S>&);
void multMatrixVector(const LargeMatrix<Matrix<S> >&, const vector<Vector<S> >&, vector<Vector<S> >&);
void multVectorMatrix(const LargeMatrix<Matrix<S> >&, const vector<Vector<S> >&, vector<Vector<S> >&);

The operator * is overloaded for LargeMatrix and vector or vector and LargeMatrix in a same way:

friend vector<S> operator*(const LargeMatrix<S>&, const vector<S>&);
friend vector<S> operator*(const vector<S>&, const LargeMatrix<S>&);
friend vector<Vector<S> operator*(const LargeMatrix<Matrix<S> >&, const vector<Vector<S> >&);
friend vector<Vector<S> operator*(const vector<Vector<S> >&, const LargeMatrix<Matrix<S> >&);

The 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):

void multMatrixMatrix(const LargeMatrix<SA>&, const LargeMatrix<SB>&, LargeMatrix<SR>&);

The LargeMatrix also provides some functions to factorize a matrix and solve the factorized system. Like multMatrixVector, these functions are external (friend of class), template with hybrid complex-real specialization.

void multMatrixMatrix(const LargeMatrix<SA>&, const LargeMatrix<SB>&, LargeMatrix<SR>&);

The LargeMatrix also provides some functions to factorize a matrix and solve the factorized system. Like multMatrixVector, these functions are external (friend of class), template with hybrid complex-real specialization.

void ldltFactorize(LargeMatrix<S>& mat);
void ldlstarFactorize(LargeMatrix<S>& mat);
void luFactorize(LargeMatrix<S>& mat);

Hint

These functions are available for SkylineStorage and DenseStorage.

If a matrix is factorized, some other operation are available

void multFactMatrixVector(const LargeMatrix<S1>&, const vector<S2>&, vector<S3>&);
void multVectorFactMatrix(const LargeMatrix<S1>&, const vector<S2>&, vector<S3>&);
void multInverMatrixVector(const LargeMatrix<S1>& mat, const vector<S2>& vec, vector<typename Conditional<NumTraits<S1>::IsComplex, S1, S2>::type>& res, FactorizationType);

In order to generate a diagonal matrix which uses an existent matrix as prototype, LargeMatrix class provides function

LargeMatrix<S> diagonalMatrix(const LargeMatrix<S>&, const S);

A special version of this function is a function to generate an identity matrix from an existing one.

LargeMatrix<S> identityMatrix(const LargeMatrix<S>&);

Hint

These functions work with LargeMatrix class having all kinds of matrix storage except RowCsStorage class and ColCsStorage class.

Finally, the class has input/output member functions :

void print(ostream&) const;
void viewStorage(ostream&) const;
ostream& operator<<(ostream&, const LargeMatrix<S>&);
void saveToFile(const String &, StorageType, bool encode=false) const;
void loadFromFile(const String&, StorageType);
String encodeFileName(const String&, StorageType) const;

Only two formats are allowed for saving matrix to file or loading matrix from file :

• the dense format (_dense), where all the matrix values are written row by row (one row by line), space character separates each value :

  A11 A12 ... A1n                       re(A11) im(A11)  ... re(A1n) im(A1n)
  A21 A22 ... A2n   if complex values   re(A21) im(A21)  ... re(A2n) im(A2n)
  ...                                   ...
  Am1 Am2 ... Amn                       re(Am1) im(Am1)  ... re(Amn) im(Amn)
  

• the coordinates format (_coo)

  i1 j1 Ai1j1                           i1 j1 re(Ai1j1) im(Ai1j1)
  i2 j2 Ai2j2     if complex values     i2 j2 re(Ai2j2) im(Ai2j2)
  ...
  

Hint

For matrix of matrices (T=Matrix<Real> or T=Matrix<Complex>), submatrix structure is not preserved when writing matrix values!

To keep some general information about matrix, with the argument encode=true in the saveFile function, the name of file may be modified by inserting the string :(m_n_storageType_valueType_strucType_p_q)

where m and n are the “dimensions” of the matrix, storageType = dense or coo valueType = real or complex , strucType = scalar or matrix and p and q are the dimensions of sub-matrices. In case of scalar value, these parameters are omitted. This trick avoids including information in file in order that it is easily loadable by other software programs, in particular Matlab.

Some examples

• Manipulate a \(3\times 2\) matrix with row dense storage:

Vector<Real> x(2,2);
LargeMatrix<Real> A(3,2,_dense,_row,1.);
out<<"matrix A : "<<A;
out<<"access A(1,1)= "<<A(1,1);
out<<"product A*x : "<<A*x;
A.saveToFile("A_dense.mat",_dense);
LargeMatrix<Real> Areload("A_dense.mat",_dense,3,2,_dense,_row);

• Manipulate a \(3\times 2\) matrix of \(2\times 2\) real matrices with dual dense storage:

Matrix<Real> I1(2,_idMatrix);
Vector<Real> x(6,3);
LargeMatrix<Matrix<Real> > B(6,3,_dense,_dual,2*I1);
out<<"product X*B : "<<X*B;

• Manipulate a \(9\times 9\) symmetric matrix with symmetric compressed storage:

//construct storage for a Q1 mesh
Number n=3;
vector<vector<Number> > elts((n-1)*(n-1),vector<Number>(4));
for(Number j=1; j<=n-1; j++)
  for(Number i=1; i<=n-1; i++) {
     Number e=(j-1)*(n-1)+i-1, p=(j-1)*n+i;
     elts[e][0]=p; elts[e][1]=p+1; elts[e][2]=p+n; elts[e][3]=p+n+1;
  }
SymCsStorage* cssym= new SymCsStorage(n*n,elts,elts);
//construct large matrix with symmetric compressed storage
LargeMatrix<Real> C(cssym,1.,_symmetric);
C.viewStorage(out);
Vector<Real> x(n*n,0.);for(Number i=0;i<n*n;i++) x1[i]=i;
out<<"product C*x : "<<C*x;
C.saveToFile("C_coo.mat"),_coo);
//construct large matrix with symmetric compressed storage
LargeMatrix<Real> C1(cssym,1.,_symmetric);
out<<"addition C+C1 "<< C+C1;
// multiply a matrix with a scalar
out<<"product C*10.0 "<< C*10.0;
// generate a diagonal matrix
out<<"diagonal matrix from matrix C " << diagonalMatrix(C, 10.0);
// generate an identity matrix
out<<"identity matrix from matrix C " << identityMatrix(C);

Note that compressed storage has to be built before to instantiate a large matrix with compressed storage. The function viewStorage produced the output:

symmetric_compressed sparse (csr,csc), size = 29, 9 rows, 9 columns
 (shared by 1 objects).
        123456789
     1  d........
     2  xd.......
     3  .xd......
     4  xx.d.....
     5  xxxxd....
     6  .xx.xd...
     7  ...xx.d..
     8  ...xxxxd.
     9  ....xx.xd
        123456789

LDLt-Factorize a positive definite 3x3 matrix with symmetric skyline storage:

LargeMatrix<Real> rResSymSkylineLdlt(inputsPathTo("matrix3x3SymPosDef.data"), _dense, rowNum, _skyline, _symmetric);
ldltFactorize(rResSymSkylineLdlt);
out << "The result of LDLt factorizing sym skyline matrix is " << rResSymSkylineLdlt << endl;

See also

• library=largeMatrix,

• header=LargeMatrix.hpp,

• implementation=LargeMatrix.cpp,

• test=test_LargeMatrixDenseStorage.cpp, test_LargeMatrixCsStorage.cpp,

• header dep= config.h, utils.h, MatrixStorage.hpp, DenseStorage.hpp, CsStorage.hpp, SkylineStorage.hpp, config.h, utils.h

The MatrixStorage class

The 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 LargeMatrix<T> class for instance).

Various matrix storages are defined by three characteristics :

• A StorageType: _dense for dense storage, _cs for compressed sparse storage (CSR/CSC) and _skyline for skyline storage;

• An AccessType: _row for a row wise access, _col for a column wise access, _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 _sym if the storage is symmetric (only the lower triangular part is described with a row wise access;

• A BuildStorageType: describing the way the storage is constructed :_undefBuild, _feBuild, _dgBuild,_diagBuild, _ecBuild, _globalBuild, _otherBuild;

• A flag telling if the storage comes from a vector to scalar conversion.

class MatrixStorage
{
  protected:
  StorageType storageType_;     // storage type
  AccessType accessType_;       // access type
  StorageBuildType buildType_;  // storage build type
  bool scalarFlag_;             // scalar conversion flag
  Number nbRows_;               // number of rows
  Number nbCols_;               // number of columns
  Number nbObjectsSharingThis_; // number of matrices sharing "this" storage
  ...
}

The number of rows and columns are counted in value type of the matrix.

It manages also a string index to identify the storage (generally encoding the row/col space pointers)

public:
 String stringId; // string identifier based on characteristic pointers

The class provides some constructors which only initialize members of the class.

MatrixStorage();
MatrixStorage(StorageType, AccessType at = _noAccess);
MatrixStorage(StorageType, AccessType, Number, Number);
virtual ~MatrixStorage();

Protected members may be read by accessors:

String name() const;
StorageType storageType() const;
AccessType accessType() const;
StorageBuildType buildType() const;
Number nbOfRows() const;
Number nbOfColumns() const;
Number numberOfObjects() const;
void objectPlus();
void objectMinus();

All storage object pointers are stored in a static vector

static vector<MatrixStorage*> theMatrixStorages; // list of all matrix storages
static void clearGlobalVector();                 // delete all MatrixStorage objects
static void printAllStorages(ostream&);          // print the list of MatrixStorage in memory
static void printAllStorages(PrintStream& os);

and some static functions dealing with are available.

When building a new storage, it is possible to look for an existing storage by using the external function:

MatrixStorage* findMatrixStorage(const String& id, StorageType st, AccessType at, StorageBuildType sb, bool scal=false, Number nbr=0, Number nbc=0);

Travelling theMatrixStorages vector this function checks if it exists a storage having the same characteristics (id, storage type, access type, build type) and may check the scalar flag and the number of rows or columns if nonzero values are given.

The class provides virtual member functions returning the sizes of storage (number of matrix entries stored) :

virtual Number size() const = 0;
Number diagonalSize() const;
virtual Number lowerPartSize() const=0;
virtual Number upperPartSize() const=0;

Most of the 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 specializations are Real, Complex, RealMatrix and ComplexMatrix.

We describe here only the :cpp:type:~xlifepp::Real specialization. There 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 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.

In 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!

virtual Number pos(Number i, Number j, SymType sym=_noSymmetry) const;
virtual void positions(const std::vector<Number>&, const std::vector<Number>&, std::vector<Number>&, bool errorOn=true, SymType=_noSymmetry) const;

To print, to save to file or load from file large matrices, the following functions are proposed :

virtual void printEntries(std::ostream&, const std::vector<Real>&, Number vb = 0, SymType sym = _noSymmetry) const;
template<typename T>
void printDenseMatrix(std::ostream&, const std::vector<T>&, SymType s=_noSymmetry) const;
void printCooMatrix(std::ostream&, const std::vector<T>&, SymType s=_noSymmetry) const;
virtual void loadFromFileDense(std::istream&, std::vector<Real>&, SymType, bool);
virtual void loadFromFileCoo(std::istream&, std::vector<Real>&, SymType, bool);
virtual void visual(std::ostream&) const;

The printEntries function print on output stream the stored matrix entries in a well suited presentation, the printDenseMatrix function print the matrix entries in a dense format as it was dense and the printCooMatrix function print the matrix entries in a coordinate format (\(i,j,A_{ij}\)). An \(m\times n\) matrix saved as dense produced a file like:

A11 A12 .... A1n                       re(A11) im(A11)  .... re(A1n) im(A1n)
A21 A22 .... A2n   if complex values   re(A21) im(A21)  .... re(A2n) im(A2n)
....                                     ....
Am1 Am2 .... Amn                       re(Am1) im(Am1)  .... re(Amn) im(Amn)

and when saved as coordinates format :

i1 j1 Ai1j1                           i1 j1 re(Ai1j1) im(Ai1j1)
i2 j2 Ai2j2     if complex values     i2 j2 re(Ai2j2) im(Ai2j2)
...                                   ....

Note that printDenseMatrix and printCooMatrix functions are template functions. Their coding is not optimal because they use the virtual pos() function and other external utilities which “interpret” the type of the template:

void printDense(ostream&, const Real&, Number);
void printDense(ostream&, const Complex&, Number);
template<typename T>
  void printDense(ostream&, const Matrix<T>&, Number);
  void printCoo(ostream&, const Matrix<T>&, Number);
void printCoo(ostream&, const Real&, Number);
void printCoo(ostream&, const Complex&, Number);
Number numberOfRows(const Real&);
Number numberOfRows(const Complex&);
Number numberOfCols(const Real&);
Number numberOfCols(const Complex&);
template<typename T> Number numberOfRows(const Matrix<T>&);
template<typename T> Number numberOfCols(const Matrix<T>& m);

The visual() function prints only the matrix entries inside the storage as ‘x’ characters (d for diagonal entries):

matrix storage :row_compressed sparse (csr,csc), size = 49, 9 rows, 9 columns (shared by 0 objects).
        123456789
     1  dx.xx....
     2  xdxxxx...
     3  .xd.xx...
     4  xx.dx.xx.
     5  xxxxdxxxx
     6  .xx.xd.xx
     7  ...xx.dx.
     8  ...xxxxdx
     9  ....xx.xd
        123456789

The MatrixStorage class provides the product matrix \(\times\) vector and the product vector \(\times\) matrix:

virtual
  void multMatrixVector(const vector<Real>&, const vector<Real>&, vector<Real>&, SymType) const;
  void multMatrixVector(const vector< Matrix<Real> >&, const vector<Vector<Real>>&, vector Vector<Real>>&, SymType) const;
  void multVectorMatrix(const vector<Real>&, const vector<Real>&, vector<Real>&, SymType) const;
  void multVectorMatrix(const vector<Matrix<Real>>&, const vector<Vector<Real>>&, vector<Vector<Real>>&, SymType) const;

Note that mixing real vector/matrix and complex vector/matrix is also provided. Each child class implements the best algorithm for the product.

The addition of two matrices is supplied with:

virtual void addMatrixMatrix(const vector<Real>&, const vector<Real>&, vector<Real>&) const;

Similar to multiplication of matrix-vector, this function is implemented in each child class. The MatrixStorage class also provides the matrix factorizations:

virtual
 void ldlt(vector<Real>&, vector<Real>&, const SymType sym = _symmetric) const;
 void ldlt(vector<Complex>&, vector<Complex>&, const SymType sym = _symmetric) const;
 void lu(vector<Real>&, vector<Real>&, const SymType sym = _noSymmetry) const;
 void lu(vector<Complex>&, vector<Complex>&, const SymType sym = _noSymmetry) const;
 void ldlstar(vector<Real>&, vector<Real>&) const;
 void ldlstar(vector<Complex>&, vector<Complex>&) const;

as well as some special solvers

virtual
 void lowerD1Solver(const vector<Real>&, vector<Real>&, vector<Real>&) const;
 void diagonalSolver(const vector<Real>&, vector<Real>&, vector<Real>&) const
 void upperD1Solver(const vector<Real>&, vector<Real>&, vector<Real>&, const SymType sym = _noSymmetry) const;
 void upperSolver(const vector<Real>&, vector<Real>&, vector<Real>&, const SymType sym = _noSymmetry) const;
 void sorDiagonalMatrixVector(const vector<Real>&, const vector<Real>&, vector<Real>&, const Real) const;
 void sorLowerMatrixVector(const vector<Real>&, const vector<Real>&, vector<Real>& ,Real, const SymType = _noSymmetry) const;
 void sorUpperMatrixVector(const vector<Real>&, const vector<Real>&, vector<Real>&, Real, const SymType = _noSymmetry) const;
 void sorDiagonalSolve(const vector<Real>&, const vector<Real>&, vector<Real>&, Real) const;
 void sorLowerSolve(const vector<Real>&, const vector<Real>&, vector<Real>&, Real) const;
 void sorUpperSolve(const vector<Real>&, const vector<Real>&, vector<Real>&, Real, const SymType sym = _noSymmetry) const;

Each descendant will implement the specific algorithm for these virtual functions.

Values on the diagonal of a matrix is set with the function

virtual void setDiagValue(vector<Real>&, Real);

See also

• library=largeMatrix,

• header=MatrixStorage.hpp,

• implementation=MatrixStorage.cpp,

• test=test_LargeMatrixDenseStorage.cpp, test_LargeMatrixCsStorage.cpp,

• header dep= config.h, utils.h

Dense storages

The DenseStorage class

Dense storages are proposed as children of MatrixStorage class in order to be consistent with other storages. Contrary to the Matrix class (see utils lib) which is a row dense storage, several dense storages are offered for an \(m\times n\) matrix (say \(A\)) :

• the row dense storage (RowDenseStorage class)

  \((A_{11}, A_{12}, \cdots, A_{1n},A_{21}, A_{22}, \cdots, A_{2n}, \cdots, A_{m1}, A_{m2}, \cdots, A_{mn})\)

• The column dense storage (ColDenseStorage class)

  \((A_{11}, A_{21}, \cdots, A_{2m},A_{12}, A_{22}, \cdots, A_{m2}, \cdots, A_{1n}, A_{2n}, \cdots, A_{mn})\)

• The dual dense storage (DualDenseStorage class), diagonal part is first stored, then strict lower triangular part and strict upper triangular part:

  \((A_{11}, A_{22}, \cdots, A_{mm}, A_{21}, \cdots, A_{i1},\cdots, A_{i,i-1}, \cdots, A_{12}, \cdots, A_{1j}, \cdots, A_{j-1,j}, \cdots) \ \) \(\text{case }m=n\)

• The symmetric dense storage (SymDenseStorage class), diagonal part is first stored, then strict lower triangular part :

  \((A_{11}, A_{22}, \cdots, A_{mm}, A_{21}, \cdots, A_{i1},\cdots, A_{i,i-1}, \cdots) \ \text{case }m=n\)

Dual storage works also for non-square matrices. Obviously, symmetric storage works only for square matrices!

All dense storages inherit from the DenseStorage abstract class which inherits from MatrixStorage abstract class.

The DenseStorage class

The DenseStorage class has no member attribute. There are some basic constructors just calling MatrixStorage constructors:

class DenseStorage : public MatrixStorage
{
  public :
    DenseStorage();
    DenseStorage(AccessType);
    DenseStorage(AccessType, Number);
    DenseStorage(AccessType, Number, Number);
    virtual ~DenseStorage() {}
    ....
}

DenseStorage objects do not have to be instantiated. This class acts as an interface to particular dense storages and gathers all common functionalities of dense storages, some being virtual:

• Public size functions:

  Number size() const;
  virtual Number lowerPartSize() const;
  virtual Number upperPartSize() const;
  

• Protected input/output functions (template line is written once):

  template<typename Iterator>
  void printScalarEntries(Iterator&, Number, Number, Number, Number, Number, const String&, Number, std::ostream&) const;
  void printScalarEntriesTriangularPart(Iterator& , Iterator&, Number, Number, Number, Number, Number, const String&, Number, std::ostream&) const;
  void printMatrixEntries(Iterator&, Number, Number, const String&, Number, std::ostream&) const;
  void printMatrixEntriesTriangularPart(Iterator&, Iterator&, Number, Number, const String&, Number, std::ostream&) const;
  void loadFromFileDense(std::istream&, std::vector<Real>&, SymType, bool);
  void loadFromFileDense(std::istream&, std::vector<Complex>&, SymType, bool);
  void loadFromFileCoo(std::istream&, std::vector<Real>&, SymType, bool);
  void loadFromFileCoo(std::istream&, std::vector<Complex>&, SymType, bool);
  

• product of matrix and vector (template line is written once):

  template<typename MatIterator, typename VecIterator, typename ResIterator>
  void diagonalMatrixVector(MatIterator&, VecIterator&, ResIterator&, ResIterator&) const;
  void diagonalVectorMatrix(MatIterator&, VecIterator&, ResIterator&, ResIterator&) const;
  void lowerMatrixVector(MatIterator&, VecIterator&, VecIterator&, ResIterator&, ResIterator&, SymType) const;
  void lowerVectorMatrix(MatIterator&, VecIterator&, VecIterator&, ResIterator&, ResIterator&, SymType) const;
  void upperMatrixVector(MatIterator&, VecIterator&, VecIterator&, ResIterator&, ResIterator&, SymType) const;
  void upperVectorMatrix(MatIterator&, VecIterator&, VecIterator&, ResIterator&, ResIterator&, SymType) const;
  void rowMatrixVector(MatIterator&, VecIterator&, VecIterator&, ResIterator&, ResIterator&)  const;
  void rowVectorMatrix(MatIterator&, VecIterator&, VecIterator&, ResIterator&, ResIterator&)  const;
  void columnMatrixVector(MatIterator&, VecIterator&, VecIterator&, ResIterator&, ResIterator&)  const;
  void columnVectorMatrix(MatIterator&, VecIterator&, VecIterator&, ResIterator&, ResIterator&)  const;
  

Note that the product of matrix and vector are performed either by row or column access, for diagonal part, lower part or upper part. They use iterators to by-pass the templated matrix values vector! With this trick, product algorithms are all located in DenseStorage class.

• Addition of two matrices:

template<typename Mat1Iterator, typename Mat2Iterator, typename ResIterator>
void sumMatrixMatrix(Mat1Iterator&, Mat2Iterator&, ResIterator&, ResIterator&) const;

See also

• library=largeMatrix,

• header=DenseStorage.hpp,

• implementation=DenseStorage.cpp,

• test=,test_LargeMatrixDenseStorage.cpp,

• header dep= MatrixStorage.hpp, config.h, utils.h

The RowDenseStorage/ColDenseStorage/DualDenseStorage/SymDenseStorage classes

All child classes (RowDenseStorage, ColDenseStorage, DualDenseStorage and SymDenseStorage) of DenseStorage class have the same structure. We give here only the description of RowDenseStorage class.

RowDenseStorage class manages dense storage where matrix values are stored row by row in values vector. It has no member attribute and has constructors just calling DenseStorage constructors:

class RowDenseStorage : public DenseStorage
{
  public :
   RowDenseStorage();
   RowDenseStorage(Number);
   RowDenseStorage(Number, Number);
   virtual ~RowDenseStorage() {}
   ....
}

The class provides all the virtual methods declared in MatrixStorage class :

Number pos(Number i, Number j, SymType s=_noSymmetry) const;
void positions(const std::vector<Number>&, const std::vector<Number>&, std::vector<Number>&, bool errorOn=true, SymType=_noSymmetry ) const;
void printEntries(std::ostream&, const std::vector<Real>&, Number vb = 0, SymType sym = _noSymmetry) const;
template<typename M, typename V, typename R>
void multMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const;
void multVectorMatrix(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const;
template<typename T>
void setDiagValueColDense(std::vector<T>& m, const T k);
template<typename M1, typename M2, typename R>
void addMatrixMatrix(const std::vector<M1>&, const std::vector<M2>&, std::vector<R>&) const;
...

Only version of multMatrixVector and multVectorMatrix are written here for Real type but versions for Complex, RealMatrix, ComplexMatrix and mixed types are also provided.

DualDenseStorage and SymDenseStorage classes provides also the size of strict lower and strict upper triangular part:

Number lowerPartSize() const;
Number upperPartSize() const;

See also

XXX=Row, Col, Dual or Sym

• library=largeMatrix,

• header=XXXDenseStorage.hpp,

• implementation=XXXDenseStorage.cpp,

• test=test_LargeMatrixDenseStorage.cpp, test_LargeMatrixCsStorage.cpp,

• header dep= DenseStorage.hpp, MatrixStorage.hpp, config.h, utils.h

Compressed sparse storages

Compressed sparse storages are designed to store only nonzero values of a matrix. It exists several methods to do it. The well known CSR/CSC storage being the most common one, we chose it. It consists in storing nonzero values of a row (resp. column) in a compressed vector and their column indices (resp. row indices) in a compressed vector. A vector of addresses where begin rows (resp. columns) is also required. More precisely, for an \(m\times n\) matrix

• Row compressed storage is defined by the three vectors:

  values = \((0, A_{1\mathcal{C}_1}, \ldots, A_{i\mathcal{C}_i}, \ldots, A_{m\mathcal{C}_m})\)

  colIndex = \((\mathcal{C}_1, \ldots, \mathcal{C}_i, \ldots, \mathcal{C}_m)\)

  rowPointer = \((p_1, \ldots, p_i, \ldots, p_m, p_{m+1})\)

  where \(\mathcal{C}_i\) is the set of index \(j\) where \(A_{ij}\ne 0\) and \(p_i\) is the position in vector colIndex of the first entry of row \(i\le m\); \(p_{m+1}\) gives the number of nonzero values.

  By convention, all the indices begins at 0 and first entry of values vector is always 0!

  For instance, the following \(5\times 6\) non-symmetric matrix

  \[\begin{split}A=\left[ \begin{array}{cccccc} 11 & 12 & 0 & 14 & 0 & 16\\ 0 & 22 & 23 & 0 & 0 & 26 \\ 31 & 32 & 33 & 0 & 0 & 0 \\ 0 & 0 & 43 & 44 & 0 & 0 \\ 51 & 0 & 53 & 54 & 55 & 0 \end{array} \right]\end{split}\]

  has the row compressed storage vectors:

  values = \((0\ 11\ 12\ 14\ 16\ \ 22\ 23\ 26\ \ 31\ 32\ 33\ \ 43\ 44\ \ 51\ 53\ 54\ 55)\)

  colIndex = \((0\ 1\ 3\ 5\ \ 1\ 2\ 5\ \ 0\ 1\ 2\ \ 2\ 3\ \ 0\ 2\ 3\ 4)\)

  rowPointer = \((0\ 4\ 7\ 10\ 12\ 16)\)

• The Col compressed storage is defined by the three vectors:

  values = \((0,\ A_{\mathcal{R}_11}, \ldots , \ A_{\mathcal{R}_jj}, \ldots, \ A_{\mathcal{R}_nn})\)

  rowIndex = \((\mathcal{R}_1, \ldots , \ \mathcal{R}_j, \ldots,\ \mathcal{R}_n)\)

  colPointer = \((p_1, \ldots ,\ p_j, \ldots ,\ p_n,\ p_{n+1})\)

  where \(\mathcal{R}_j\) is the set of index \(i\) where \(A_{ij}\ne 0\) and \(p_j\) is the position in rowIndex vector of the first entry of column \(j\le n\); \(p_{n+1}\) gives the number of nonzero values.

  By convention, all the indices begins at 0 and first entry of values vector is always 0!

  For instance, the previous \(5\times 6\) non-symmetric matrix \(A\) has the column compress storage vectors:

  values = \((0\ 11\ 31\ 51\ \ 12\ 22\ 32\ \ 23\ 33\ 43\ 53\ \ 55\ \ 16\ 26)\)

  rowIndex = \((0\ 2\ 4\ \ 0\ 1\ 2\ \ 1\ 2\ 3\ 4\ \ 0\ 3\ 4\ \ 4\ \ 0\ 1)\)

  colPointer = \((0\ 4\ 7\ 10\ 12\ 16)\)

• The Dual compressed storage is defined by the five vectors:

  values = \((0,\text{diag}(A), A_{\mathcal{C}_11}, \ldots, A_{\mathcal{C}_ii}, \ldots, A_{\mathcal{C}_mm}, A_{\mathcal{R}_11}, \ldots, A_{\mathcal{R}_jj}, \ldots, A_{\mathcal{R}_nn})\)

  colIndex = \((\mathcal{C}_1, \ldots, \mathcal{C}_i, \ldots, \mathcal{C}_m)\)

  rowPointer = \((p_1, \ldots ,\ p_i,\ldots, p_m, p_{m+1})\)

  rowIndex = \((\mathcal{R}_1, \ldots, \mathcal{R}_j, \ldots, \mathcal{R}_n)\)

  colPointer = \((q_1, \ldots, q_j, \ldots, ,q_n, q_{n+1})\)

  where \(\mathcal{C}_i\) is the set of index \(j\) where \(A_{ij}\ne 0\) and \(j\le\min(i-1,n)\), \(\mathcal{R}_j\) is the set of index \(i\) where \(A_{ij}\ne 0\) and \(i\le\min(j-1,m)\), \(p_i\) (resp. \(q_j\)) is the position in colIndex (resp. rowIndex) vector of the first entry of row \(i\le m\) (resp. column \(j\le n\)). The colIndex and rowPointer vectors describes the row compressed storage of the strict lower part of matrix while the rowIndex and colPointer vectors describes the column compressed storage of the strict upper part of matrix.

  In this storage, all the diagonal entries of matrix are always first stored even the null entries!

  For instance, the previous \(5\times 6\) non-symmetric matrix \(A\) has the column compress storage vectors:

  values = \((0\ 11\ 22\ 33\ 44\ 55\ \ 31\ 32\ \ 43\ \ 51\ 53\ 54\ \ \ 23\ \ \ 14\ 44\ \ \ 16\ 26)\)

  colIndex = \((0\ 1\ \ 2\ \ 0\ 2\ 3)\)

  rowPointer = \((0\ 0\ 0\ 2\ 3\ 6)\)

  rowIndex = \((0\ \ 1\ \ 0\ \ 0\ 1)\)

  colPointer = \((0\ 0\ 1\ 2\ 2\ 3\ 5)\)

  Note

  For an empty row \(i\) (resp. column \(j\)), rowPointer(\(i\))=rowPointer(\(i-1\)) (resp. colPointer(\(j\))=colPointer(\(j-1\))).

  By the way, rowPointer(\(i\))-rowPointer(\(i-1\)) gives the number of stored entries on row \(i-1\).

• The Symmetric compressed storage is devoted to square matrix with symmetric storage, the matrix may be not symmetric.

  In that case, only the storage of lower triangular part of matrix is described by the three vectors (those of row part of dual storage):

  values = \((0,\text{diag}(A), A_{\mathcal{C}_11}, \ldots, A_{\mathcal{C}_ii}, \ldots, A_{\mathcal{C}_mm})\) if \(A\) has a symmetry property

  values = \((0,\text{diag}(A), A_{\mathcal{C}_11}, \ldots, A_{\mathcal{C}_ii}, \ldots, A_{\mathcal{C}_mm}, A_{\mathcal{C}_11}, \ldots, A_{\mathcal{C}_jj}, \ldots, A_{\mathcal{C}_mm})\) if \(A\) has no symmetry property

  colIndex = \((\mathcal{C}_1, \ldots, \mathcal{C}_i, \ldots, \mathcal{C}_m)\)

  rowPointer = \((p_1, \ldots, p_i, \ldots ,p_m, p_{m+1})\)

  where the meaning of \(\mathcal{C}_i\) and \(p_i\) is the same as dual storage.

  For instance, the \(5\times 5\) non-symmetric matrix with symmetric storage

  \[\begin{split}B=\left[ \begin{array}{ccccc} 11 & 0 & 13 & 0 & 15 \\ 0 & 22 & 23 & 0 & 0 \\ 31 & 32 & 33 & 34 & 35\\ 0 & 0 & 43 & 44 & 45 \\ 51 & 0 & 53 & 54 & 55 \end{array} \right]\end{split}\]

  has the symmetric compressed storage vectors:

  values = \((0\ 11\ 22\ 33\ 44\ 55\ \ 31\ 32\ \ 43\ \ 51\ 53\ 54\ \ 13\ 23\ \ 34\ \ 15\ 35\ 45)\)

  colIndex = \((0\ 1\ \ 2\ \ 0\ 2\ 3)\)

  rowPointer = \((0\ 0\ 0\ 2\ 3\ 6)\)

  Warning

  Be careful with the fact that actual matrix values begins at address 1 in values vector defined outside the storage stuff.

The CsStorage class

The CsStorage class has no member attribute. There are some basic constructors just calling MatrixStorage constructors:

class CsStorage : public MatrixStorage
{
  public :
    CsStorage(AccessType = _dual);
    CsStorage(Number, AccessType = _dual);
    CsStorage(Number, Number, AccessType = _dual);
    virtual ~CsStorage() {}
  protected :
    void buildCsStorage(const std::vector<std::vector<Number> >&, std::vector<Number>&, std::vector<Number>&);
  ...
}

The buildCsStorage() is a general member function building storage vectors from lists of column indices by row. It is used by any child class of CsStorage class.

CsStorage objects do not have to be instantiated. This class acts as an interface to particular compressed storages and gathers all common functionalities of dense storages, some being virtual (some template line declaration are omitted):

virtual Number size() const = 0;
virtual Number lowerPartSize() const;
virtual Number upperPartSize() const;
template<typename Iterator>
  void printEntriesAll(StrucType, Iterator& , const std::vector<Number>&, const std::vector<Number>&, Number, Number, Number, const String&, Number, std::ostream&) const;
  void printEntriesTriangularPart(StrucType, Iterator&, Iterator&, const std::vector<Number>&, const std::vector<Number>&, Number, Number, Number, const String&, Number, std::ostream&) const;
  void printCooTriangularPart(std::ostream&, Iterator&, const std::vector<Number>&, const std::vector<Number>&, bool, SymType sym = _noSymmetry) const;
  void printCooDiagonalPart(std::ostream&, Iterator&, Number) const;

template <typename T>
  void loadCsFromFileDense(std::istream&, std::vector<T>&, std::vector<Number>&, std::vector<Number>&, SymType, bool);
  void loadCsFromFileCoo(std::istream&, std::vector<T>&, std::vector<Number>&, std::vector<Number>&, SymType, bool);
  void loadCsFromFileDense(std::istream&, std::vector<T>&, std::vector<Number>&, std::vector<Number>&, std::vector<Number>&, std::vector<Number>&, SymType, bool);
  void loadCsFromFileCoo(std::istream&, std::vector<T>&, std::vector<Number>&, std::vector<Number>&, std::vector<Number>&, std::vector<Number>&, SymType, bool);

The CsStorage class provides the real code of product of matrix and vector, using iterators as arguments to by-pass the template type of matrix values. For dual and symmetric storage, the computation is performed by splitting the matrix in its diagonal, lower and upper part (template line are omitted):

template<typename MatIterator, typename VecIterator, typename ResIterator>
  void diagonalMatrixVector(MatIterator&, VecIterator&, ResIterator&, ResIterator&) const;
  void lowerMatrixVector(const std::vector<Number>&, const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&, SymType sym) const;
  void upperMatrixVector(const std::vector<Number>&, const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&, SymType) const;
  void rowMatrixVector(const std::vector<Number>&, const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&)  const;
  void columnMatrixVector(const std::vector<Number>&, const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&)  const;
  void diagonalVectorMatrix(MatIterator&, VecIterator&, ResIterator&, ResIterator&) const;
  void lowerVectorMatrix(const std::vector<Number>&, const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&, SymType sym) const;
  void upperVectorMatrix(const std::vector<Number>&, const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&, SymType) const;
  void rowVectorMatrix(const std::vector<Number>&, const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&)  const;
  void columnVectorMatrix(const std::vector<Number>&, const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&)  const;

Besides the product of matrix and vector, the CsStorage provides solvers for matrix in its diagonal, lower and upper part

template<typename MatIterator, typename VecIterator, typename ResIterator>
  void bzSorDiagonalMatrixVector(MatIterator& itd, const VecIterator& itvb, ResIterator& itrb, const Real w) const;
template<typename MatIterator, typename VecIterator, typename XIterator>
  void bzSorLowerSolver(const MatIterator&, const MatIterator&, VecIterator&, XIterator&, XIterator&, const std::vector<Number>&, const std::vector<Number>&, const Real) const;
template<typename MatIterator, typename VecIterator, typename XIterator>
  void bzSorDiagonalSolver(const MatIterator& itdb, VecIterator& itbb, XIterator& itxb, XIterator& itxe, const Real w) const;
template<typename MatRevIterator, typename VecRevIterator, typename XRevIterator>
  void bzSorUpperSolver(const MatRevIterator&, const MatRevIterator&, VecRevIterator&, XRevIterator&, XRevIterator&, const std::vector<Number>&, const std::vector<Number>&, const Real) const;

Addition of two matrices are implemented by

template<typename Mat1Iterator, typename Mat2Iterator, typename ResIterator>
  void sumMatrixMatrix(Mat1Iterator&, Mat2Iterator&, ResIterator&, ResIterator&)

See also

• library=largeMatrix,

• header=CsStorage.hpp,

• implementation=CsStorage.cpp,

• test=test_LargeMatrixCsStorage.cpp,

• header dep= MatrixStorage.hpp, config.h, utils.h

The RowCsStorage / ColCsStorage classes

As mentioned before, row compressed storage are based on the column indices vector and the row pointers vector. So the RowCsStorage class carries these two vectors:

class RowCsStorage : public CsStorage
{
  protected :
    std::vector<Number> colIndex_; // vector of column index of nonzero value
    std::vector<Number> rowPointer_; // vector of positions of begining of rows
  ...
}

The class provides a default constructor with no construction of storage vectors and constructors from a list of \((i,j)\) index (coordinates) or a list of column indices by row, constructing the compressed storage vectors:

RowCsStorage(Number nr=0, Number nc=0);
RowCsStorage(Number, Number, const std::vector< std::vector<Number> >&, const std::vector< std::vector<Number> >&);
RowCsStorage(Number, Number, const std::vector< std::vector<Number> >&);
~RowCsStorage() {};

The class provides most of the virtual methods declared in MatrixStorage class. For sake of simplicity, we report here only function with Real type but versions with Complex, Matrix<Real>, Matrix<Complex> and mixed types are also provided. Some template line declarations are also omitted.

Number size() const;
Number pos(Number i, Number j, SymType s=_noSymmetry) const;
void positions(const std::vector<Number>&, const std::vector<Number>&, std::vector<Number>&, bool errorOn=true, SymType s=_noSymmetry) const;
void printEntries(std::ostream&, const std::vector<Real>&, Number, SymType) const;
void printCooMatrix(std::ostream&, const std::vector<Real>&, SymType s=_noSymmetry) const;
void loadFromFileDense(std::istream&, std::vector<Real>&, SymType, bool );
void loadFromFileCoo(std::istream& , std::vector<Real>&, SymType, bool);
template<typename M, typename V, typename R>
  void multMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const;
  void multVectorMatrix(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const;
void multMatrixVector(const std::vector<Real>&, const std::vector<Real>&, std::vector<Real>&, SymType) const;
void multVectorMatrix(const std::vector<Real>&, const std::vector<Real>&, std::vector<Real>&, SymType) const;
template<typename M1, typename M2, typename R>
  void addMatrixMatrix(const std::vector<M1>&, const std::vector<M2>&, std::vector<R>&) const; //!< templated row dense Matrix + Matrix
void addMatrixMatrix(const std::vector<Real>& m, const std::vector<Real>& v, std::vector<Real>& rv) const
  { addMatrixMatrix<>(m, v, rv);}

There is no printDenseMatrix() function because the general one defined in MatrixStorage class is sufficiently efficient.

The ColCsStorage class is very similar to the RowCsStorage class, except that it manages a row indices vector and a column pointers vector:

class ColCsStorage : public CsStorage
{
 protected :
  std::vector<Number> rowIndex_; // vector of row index of nonzero value
  std::vector<Number> colPointer_; // vector of positions of begining of cols
  ...
}

Its constructors have a similar structure (same arguments) and all prototypes of member functions are the same as ColCsStorage member functions. We do not repeat them.

See also

XXX=Row or Col

• library=largeMatrix,

• header=XXXCsStorage.hpp,

• implementation=XXXCsStorage.cpp,

• test=test_LargeMatrixCsStorage.cpp,

• header dep= CsStorage.hpp, MatrixStorage.hpp, config.h, utils.h

The DualCsStorage / SymCsStorage classes

Dual compressed storage travels the matrix as diagonal, strict lower and upper parts. Column indices and row pointers vectors are attached to strict lower part and row indices and column pointers vectors are attached to strict upper part:

class DualCsStorage : public CsStorage
{
  protected :
   std::vector<Number> colIndex_;   //vector of column index of nonzero value
   std::vector<Number> rowPointer_; //vector of positions of begining of rows
   std::vector<Number> rowIndex_;   //vector of row index of nonzero value
   std::vector<Number> colPointer_; //vector of positions of begining of cols
 ...
}

DualCsStorage class provides basic constructor (no construction of storage vectors), constructor from a pair of global numbering vectors or list of column indices by rows.

Both of them create the compressed storage vector using the auxiliary function buildStorage():

DualCsStorage(Number nr=0, Number nc=0);
DualCsStorage(Number, Number, const std::vector< std::vector<Number> >&, const std::vector< std::vector<Number> >&);
DualCsStorage(Number, Number, const std::vector<std::vector<Number> >&);
void buildStorage(const std::vector<std::vector<Number> >&);
~DualCsStorage() {};

The class provides most of the virtual methods declared in MatrixStorage class. For sake of simplicity, we report here only function with Real type but versions with Complex, RealMatrix, ComplexMatrix and mixed types are also provided. Some template line declarations are also omitted.

Number size() const;
Number lowerPartSize() const ;
Number upperPartSize() const;
Number pos(Number i, Number j, SymType s=_noSymmetry) const;
void positions(const std::vector<Number>&, const std::vector<Number>&, std::vector<Number>&, bool errorOn=true, SymType s=_noSymmetry) const;
void printEntries(std::ostream&, const std::vector<Real>&, Number, SymType) const;
void printCooMatrix(std::ostream&, const std::vector<Real>&, SymType s=_noSymmetry) const;
void loadFromFileDense(std::istream&, std::vector<Real>&, SymType, bool );
void loadFromFileCoo(std::istream& , std::vector<Real>&, SymType, bool);
template<typename M, typename V, typename R>
  void multMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const;
void multVectorMatrix(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const;
void multMatrixVector(const std::vector<Real>&, const std::vector<Real>&, std::vector<Real>&, SymType) const;
void multVectorMatrix(const std::vector<Real>&, const std::vector<Real>&, std::vector<Real>&, SymType) const;
template<typename M1, typename M2, typename R>
  void addMatrixMatrix(const std::vector<M1>&, const std::vector<M2>&, std::vector<R>&) const;
void addMatrixMatrix(const std::vector<Real>& m, const std::vector<Real>& v, std::vector<Real>& rv) const
  { addMatrixMatrix<>(m, v, rv);}

The SymCsStorage class is very similar to the DualCsStorage class, except that it addresses only square matrix with symmetric storage and, thus, does not manage row indices vector and column pointers vectors because the upper part has same storage as lower part (transposed).

Note that matrix may have no symmetry property but a symmetric storage. In that case, the upper part is stored.

class SymCsStorage : public CsStorage
{
  protected :
   std::vector<Number> colndex_;   //!< vector of row index of nonzero value
   std::vector<Number> rowPointer_; //!< vector of positions of begining of cols
  ...
}

Contrary to the row compressed storage, in the SymCsStorage, the diagonal is first stored then the lower part of the matrix.

Its constructors have a similar structure (same arguments) to the DualCsStorage constructors and all prototypes of member functions are the same as DualCsStorage member functions. We do not repeat them.

Different from other CsStorage classes, the SymCsStorage class provides some functions to deal with the SorSolver class and SsorSolver class. Only specialized functions with Real are listed in the following

void sorDiagonalMatrixVector(const std::vector<Real>& m, const std::vector<Real>& v, std::vector<Real>& rv, const Real w) const
{ sorDiagonalMatrixVector<>(m, v, rv, w); }
void sorLowerMatrixVector(const std::vector<Real>& m, const std::vector<Real>& v, std::vector<Real>& rv, const Real w, const SymType sym) const
{ sorLowerMatrixVector<>(m, v, rv, w, sym); }
void sorUpperMatrixVector(const std::vector<Real>& m, const std::vector<Real>& v, std::vector<Real>& rv, const Real w, const SymType sym) const
{ sorUpperMatrixVector<>(m, v, rv, w, sym); }
void sorDiagonalSolve(const std::vector<Real>& m, const std::vector<Real>& b, std::vector<Real>& x, const Real w) const
{ sorDiagonalSolve<>(m, b, x, w); }
void sorLowerSolve(const std::vector<Real>& m, const std::vector<Real>& b, std::vector<Real>& x, const Real w) const
{ sorLowerSolve<>(m, b, x, w); }
void sorUpperSolve(const std::vector<Real>& m, const std::vector<Real>& b, std::vector<Real>& x, const Real w, const SymType sym) const
{ sorUpperSolve<>(m, b, x, w, sym); }

See also

XXX=Dual or Sym

• library=largeMatrix,

• header=XXXCsStorage.hpp,

• implementation=XXXCsStorage.cpp,

• test=test_LargeMatrixCsStorage.cpp,

• header dep= CsStorage.hpp, MatrixStorage.hpp, config.h, utils.h

Skyline storages

With dense storages, skyline storages are the only storages that are compatibles with factorization methods like LU. Skyline storages consist in storage rows or columns from its first nonzero value up to its last nonzero value. When lower triangular part is addressed, the last nonzero value on a row is the diagonal matrix entry. When upper triangular part is addressed, the last nonzero value on a column is the diagonal matrix entry. Only dual and symmetric skyline storages are proposed here (resp. DualSkylineStorage and SymSkylineStorage classes).

• 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):

  values = \((0, A_{11}, \ldots,A_{ii}, \ldots, A_{il},\ldots, A_{i,p_i}, \ldots, A_{i,l_i},\ldots, A_{q_j,j}, \ldots, A_{q_j,c_j})\)

  rowPointer = \((s_1,\ldots,s_i,\ldots,s_m,s_{m+1})\)

  colPointer = \((t_1,\ldots,t_j,\ldots,t_n,t_{n+1})\)

  where

  \[\begin{split}\begin{array}{l} l_i=\min(i-1,n-1),\ c_j=\min(j-1,m-1)\\ p_i=\min{0<j<l_i,\ A_{ij}\ne 0} \mathrm{\;and\;} q_j=\min{0<i<c_j,\ A_{ij}\ne 0}\\ s_1=0;\ s_i=s_{i-1}+i-s_{i-1}\ \forall i=2,m, \mathrm{\;and\;} t_1=0;\ t_j=s_{j-1}+j-t_{j-1}\ \forall j=2,n \end{array}\end{split}\]

  The rowPointer (resp. colPointer) vector gives the position in lower (resp. upper) part of the first nonzero 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).

  For instance, the following \(5\times 6\) non-symmetric matrix

  \[\begin{split}A=\left[ \begin{array}{cccccc} 11 & 12 & 0 & 14 & 0 & 16\\ 0 & 22 & 23 & 0 & 0 & 26 \\ 31 & 32 & 33 & 0 & 0 & 0 \\ 0 & 0 & 43 & 44 & 0 & 0 \\ 51 & 0 & 53 & 54 & 55 & 0 \end{array} \right]\end{split}\]

  has the following storage vectors:

  values = \((0\ 11\ 22\ 33\ 44\ 45\ \ 22\ \ 31\ 32\ \ 43\ \ 51\ 0\ 53\ 54\ \ \ 12\ \ 23\ \ 14\ 0\ 0\ \ \ 16\ 26\ 0\ 0)\)

  rowPointer = \((0\ 0\ 0\ 2\ 3\ 7)\)

  colPointer = \((0\ 0\ 1\ 2\ 5\ 5\ 8)\)

The length 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\)):

\[\begin{split}\begin{array}{l} \text{if } i=j\ \ \text{adr}(i,j)=i \\ \text{if } p_i\le j\le l_i \ \ \text{adr}(i,j)=\min(m,n)+s_{i+1}+j-i\\ \text{if } q_j\le i\le c_j \ \ \text{adr}(i,j)=\min(m,n)+s_{m+1}+t_{j+1}+i-j\\ \end{array}\end{split}\]

The SkylineStorage class

The SkylineStorage class has no member attribute. It has some basic constructors just calling MatrixStorage constructors:

class SkylineStorage : public MatrixStorage
{
  public :
   SkylineStorage(AccessType = _dual);
   SkylineStorage(Number, AccessType = _dual);
   SkylineStorage(Number, Number, AccessType = _dual);
   virtual ~SkylineStorage() {}
  ...
}

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):

virtual Number size() const = 0;
virtual Number lowerPartSize() const = 0;
virtual Number upperPartSize() const {return 0;}
template<typename Iterator>
 void printEntriesTriangularPart (StrucType, Iterator&, Iterator&, const std::vector<Number>&, Number, Number, Number, const String&, Number, std::ostream&) const;
 void printCooTriangularPart(std::ostream&, Iterator&, const std::vector<Number>&, Number, Number, bool, SymType sym = _noSymmetry) const;
 void printCooDiagonalPart(std::ostream&, Iterator&, Number) const;

template <typename T>
 void loadSkylineFromFileDense(std::istream&, std::vector<T>&, std::vector<Number>&, std::vector<Number>&, SymType, bool);
 void loadSkylineFromFileCoo(std::istream&, std::vector<T>&, std::vector<Number>&, std::vector<Number>&, SymType, bool);

template<typename MatIterator, typename VecIterator, typename ResIterator>
 void diagonalMatrixVector(MatIterator&, VecIterator&, ResIterator&, ResIterator&) const;
 void lowerMatrixVector(const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&, SymType sym) const;
 void upperMatrixVector(const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&, SymType) const;
 void diagonalVectorMatrix(MatIterator&, VecIterator&, ResIterator&, ResIterator&) const;
 void lowerVectorMatrix(const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&, SymType sym) const;
 void upperVectorMatrix(const std::vector<Number>&, MatIterator&, VecIterator&, ResIterator&, SymType) const;
 void sumMatrixMatrix(Mat1Iterator&, Mat2Iterator&, ResIterator&, ResIterator&) const;

The 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.

In order to deal with some specific solvers, this class also comes with several following functions. For simplicity, we only declare template line once.

template<typename MatIterator, typename VecIterator, typename XIterator>
 void bzLowerSolver(const MatIterator&, const MatIterator&, const VecIterator&,  const XIterator&, const XIterator&, const std::vector<size_t>::const_iterator) const;
 void bzLowerD1Solver(const MatIterator&, const VecIterator&, const XIterator&, const XIterator&, const std::vector<size_t>::const_iterator) const;
 void bzLowerConjD1Solver(const MatIterator&, const VecIterator&,  const XIterator&, const XIterator&,  const std::vector<size_t>::const_iterator) const;
 void bzDiagonalSolver(const MatIterator&, const VecIterator&,  const XIterator&, const XIterator&) const;
 void bzUpperSolver(const MatRevIterator&, const MatRevIterator&, const VecRevIterator&,  const XRevIterator&, const XRevIterator&, const std::vector<size_t>::const_reverse_iterator) const;
 void bzUpperD1Solver(const MatRevIterator&, const VecRevIterator&, const XRevIterator&, const XRevIterator&, const std::vector<size_t>::const_reverse_iterator) const;
 void bzUpperConjD1Solver(const MatRevIterator&, const VecRevIterator&, const XRevIterator&, const XRevIterator&, const std::vector<size_t>::const_reverse_iterator) const;

Like the product of matrix and vector, these functions take advantage of iterators as arguments to by-pass the template type of matrix values.

See also

• library=largeMatrix,

• header=SkylineStorage.hpp,

• implementation=SkylineStorage.cpp,

• test=test_LargeMatrixSkylineStorage.cpp,

• header dep= MatrixStorage.hpp, config.h, utils.h

The DualSkylineStorage/SymSkylineStorage classes

Dual 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:

class DualSkylineStorage : public SkylineStorage
{
  protected :
   std::vector<Number> rowPointer_; // vector of positions of begining of rows
   std::vector<Number> colPointer_; // vector of positions of begining of cols
  ...
}

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.

Both of them create the compressed storage vector using the auxiliary function buildStorage():

DualSkylineStorage(Number nr = 0, Number nc = 0);
DualSkylineStorage(Number, Number, const std::vector< std::vector<Number> >&, const std::vector< std::vector<Number> >&);
DualSkylineStorage(Number, Number, const std::vector< std::vector<Number> >&);
~DualSkylineStorage() {};

The class provides most of the virtual methods declared in MatrixStorage class. For sake of simplicity, we report here only function with Real type but versions with Complex, RealMatrix, ComplexMatrix and mixed types are also provided. Some template line declarations are also omitted.

Number size() const;
Number lowerPartSize() const ;
Number upperPartSize() const;
Number pos(Number i, Number j, SymType s=_noSymmetry) const;
void positions(const std::vector<Number>&, const std::vector<Number>&, std::vector<Number>&, bool errorOn=true, SymType s=_noSymmetry) const;
void printEntries(std::ostream&, const std::vector<Real>&, Number, SymType) const;
void printCooMatrix(std::ostream&, const std::vector<Real>&, SymType s=_noSymmetry) const;
void loadFromFileDense(std::istream&, std::vector<Real>&, SymType, bool );
void loadFromFileCoo(std::istream& , std::vector<Real>&, SymType, bool);
template<typename M, typename V, typename R>
  void multMatrixVector(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const;
void multVectorMatrix(const std::vector<M>&, const std::vector<V>&, std::vector<R>&) const;
void multMatrixVector(const std::vector<Real>&, const std::vector<Real>&, std::vector<Real>&, SymType) const;
void multVectorMatrix(const std::vector<Real>&, const std::vector<Real>&, std::vector<Real>&, SymType) const;
template<typename M1, typename M2, typename R>
  void addMatrixMatrix(const std::vector<M1>&, const std::vector<M2>&, std::vector<R>&) const;
template<typename T>
  void setDiagValueDualSkyline(std::vector<T>&, const T);
template<typename M>
  void lu(std::vector<M>& m, std::vector<M>& fa) const;
template<typename M, typename V, typename X>
  void lowerD1Solver(const std::vector<M>&, std::vector<V>&, std::vector<X>&) const;
  void diagonalSolver(const std::vector<M>&, std::vector<V>&, std::vector<X>&) const;
  void upperD1Solver(const std::vector<M>&, std::vector<V>&, std::vector<X>&, const SymType) const;
  void upperSolver(const std::vector<M>&, std::vector<V>&, std::vector<X>&) const;

The SymSkylineStorage class is very similar to the DualSkylineStorage class, except that it addresses 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) :

class SymSkylineStorage : public SkylineStorage
{
  protected :
   std::vector<Number> rowPointer_; //!< vector of positions of begining of cols
  ...
}

Note that the matrix may be not symmetric. In that case its upper part is stored.

Besides some similar functions to DualSkylineStorage class, the SymSkylineStorage provides more several functions to factorize matrix: LDLt and LDL* factorization

template<typename M>
  void ldlt(std::vector<M>& m, std::vector<M>& fa, const SymType sym = _symmetric) const;
  void ldlstar(std::vector<M>& m, std::vector<M>& fa) const;

See also

XXX=Dual or Sym

• library=largeMatrix,

• header=XXXSkylineStorage.hpp,

• implementation=XXXSkylineStorage.cpp,

• test=test_LargeMatrixSkylineStorage.cpp,

• header dep= SkylineStorage.hpp, MatrixStorage.hpp, config.h, utils.h