# Block-diagonal and banded matrices

Some basic classes are available for block-diagonal matrices (symmetric or not). Two implementations of those matrices are completed in files `src/Algebra/MatrixBlockDiagonal.cxx` and `src/Algebra/OptMatrixBlockDiagonal.cxx`. This last implementation is faster especially for the inversion of block-diagonal matrices, since it uses the class `TinyMatrix`. Otherwise the efficiency is similar. You can check in `MontjoieHeader.hxx` which implementation is included. An example of use of thoses matrices is available in file src/Program/Test/MatrixBlockDiagTest.cc.

An implementation of banded matrices and arrow matrices is proposed. Arrow matrices are banded matrices with additional dense rows and columns located at the end. These two type of matrices are particularly popular for 1-D problems, since the LU factorisation without pivoting does keep the same profile. Only non-symmetric versions of these matrices are implemented.

An implementation of diagonal matrices is provided in file DiagonalMatrix.hxx, it can be used to compute the diagonal of a finite element matrix, since the functions common to other matrices are implemented (e.g. AddInteraction, AddInteractionRow). All values that are outside the diagonal are ignored when these functions are invoked. Detailed documentation is available in this section.

## Declaration

```// block-diagonal matrix
Matrix<double, General, BlockDiagRow> A;
// symmetric block-diagonal matrix
Matrix<double, Symmetric, BlockDiagRowSym> B;

// banded matrix
Matrix<double, General, BandedCol> C;
// arrow matrix
Matrix<double, General, ArrowCol> D;

// diagonal matrix
Matrix<double, Symmetric, DiagonalRow> E;
```

## Use

These matrices can be constructed either by specifying the initial pattern, either by converting a sparse matrix into a block-diagonal matrix (respectively a banded matrix). Access to elements is achieved through the `operator(int)`, and indices start from 0 :

```// matrix is constructed through a conversion procedure
Matrix<double, Symmetric, BlockDiagRowSym> A;
Matrix<double, General, ArrayRowSparse> Asp(5, 5);
// you fill Asp as you want then you convert it to a block-diagonal matrix :
Copy(Asp, A);

// alternative way:
// the pattern is directly provided by specifying row numbers for each block
Vector<IVect> RowNum(2);
RowNum(0).Reallocate(2); RowNum(0)(0) = 0; RowNum(0)(1) = 2;
RowNum(1).Reallocate(3); RowNum(1)(0) = 1; RowNum(1)(1) = 3; RowNum(1)(1) = 5;
// previous values of A will be removed
A.SetPattern(RowNum);
// to fill the matrix, use AddInteraction

// you can display an entry
cout<<"A(3, 5)" << A(3,5) <lt; endl;

// for banded/arrow matrices, you can either convert a sparse matrix
// either fill directly the matrix with AddInteraction
// in that last case, you must give kl and ku in Reallocate
Matrix<double, General, BandedCol> Ab;
// m : number of rows, n : number of columns, kl : number of diagonals
// below the main diagonal, ku : number of diagonals above
Ab.Reallocate(m, n, kl, ku);

// then you an fill the matrix by using AddInteraction
// the abort signal is called if the value is outside the band
```

The methods AddInteraction/AddInteractionRow are quite different from those in sparse matrices, since if the asked entry does not belong to the sparsity pattern of the matrix, the value will be ignored, whereas for sparse matrices, this method creates a new non-zero entry in the case it does not exist. Unsymmetric block-diagonal matrices are handled but the sparsity pattern is assumed to be symmetric, only values may be not symmetric. For banded and arrow matrices, if the value is outside the pattern, the computation is aborted.

## Methods for block-diagonal matrices :

 Operators for block-diagonal matrices GetM returns the number of rows GetN returns the number of columns Clear clears the matrix ClearRow clears a row of the matrix ClearColumn clears a column of the matrix GetMemorySize returns the size used by the object in bytes GetRowSize returns the number of non-zero entries on row i GetSizeMaxBlock returns the maximal size of a block of the matrix GetBlockNumber returns the number of the block associated with row i SetPattern clears the matrix and specifies a new sparsity pattern Copy copies another block-diagonal matrix Reallocate changes the number of rows/columns of the matrix Resize changes the number of rows/columns of the matrix, previous rows are kept GetDataSize returns the number of elements in the matrix GetNbBlocks returns the number of blocks of the matrix GetBlockSize returns the size of a block of the matrix GetNbRowNumbers returns the number of row numbers (including distant rows in parallel) GetRowNumData returns pointer containing row numbers GetRowPtr returns pointer containing start indices of blocks GetReverseNum returns global to local array Index returns the global row number of a local row of a block Value returns the local entry of a block GetData returns the pointer to all the values of the matrix GetPtr returns start indices of blocks GetInd returns row numbers GetRev returns global to local array GetNumBlock returns the array containing block numbers to which each row belongs to GetOffset returns offset array to use for accessing a given block in data array. Fill sets all non-zero entries to a given value FillRand fills randomly the non-zero entries Zero sets all non-zero entries to zero Get access to element (i, j) of the matrix Set sets element (i, j) of the matrix AddInteraction adds a value to a non-zero entry (if present) AddInteractionRow adds values to non-zero entries of a given row WriteText prints the matrix in a file

## Functions for block-diagonal matrices :

 Mlt multiplies the elements of the matrix by a scalar MltAdd performs a matrix vector product Add adds two matrices GetInverse computes the inverse of a matrix Copy conversion to/from a sparse matrix ConvertToSparse conversion to a sparse matrix ConvertToBlockDiagonal conversion to a block-diagonal matrix GetCholesky computes the cholesky factorisation of the block-diagonal matrix SolveCholesky solves linear system L x = b or LT x = b, once GetCholesky has been called MltCholesky computes y = L x or t = LT x, once GetCholesky has been called

## Methods for Banded/Arrow matrices :

 Operators for banded or arrow matrices GetM returns the number of rows GetN returns the number of columns GetDataSize returns the number of elements in the matrix GetMemorySize returns the number of elements in the matrix GetKL returns the number of non-zero diagonals strictly below the main diagonal GetKU returns the number of non-zero diagonals strictly above the main diagonal GetNbLastRow returns the number of dense rows at the end of the matrix GetNbLastCol returns the number of dense columns at the end of the matrix GetData returns the pointer to all the values of the matrix (banded matrices only) Clear clears the matrix Zero sets to zero the values of non-zero entries Reallocate sets the size of the matrix the bandwidth, the number of last rows and columns AddInteraction adds a value to a non-zero entry (if present) AddInteractionRow adds values to non-zero entries of a given row ClearRow sets to zero the values on a row of the matrix Val returns direct access to an element of the matrix Get returns access to an element of the matrix Set sets an element of the matrix SetIdentity sets the matrix to the identity matrix Fill all non-zero entries are affected a same value FillRand each non-zero entry is affected a random value Copy copies the contents of a matrix into the current matrix Factorize the matrix is replaced by its LU factors Solve solves A x = b, assuming that Factorize has been called previously Write the matrix is written in binary format WriteText the matrix is written in text format

## Functions for banded matrices or arrow matrices :

 Mlt multiplies the elements of the matrix by a scalar or performs a matrix vector product Copy converts a sparse matrix into a banded matrix Add adds two matrices MltAdd performs a matrix vector product GetLU computes LU factors (by using Lapack or not) SolveLU solves the linear system A x = b, assuming that GetLU has been previously called

### Operators for block-diagonal or banded matrices

#### Syntax

``` T operator()(int i, int j) const;
operator = (const Matrix)
operator *= (const T& alpha)
```

For block-diagonal matrices, the operator () can be used to modify the matrix, whereas for banded and arrow matrices, the methods Get or Val have to be used.

#### Example :

```Matrix<float, General, BlockDiagRow> A, B;

// for block-diagonal matrices, the operator () can be used to modify the matrix
A(2, 0) = 1.5;

// operator = to copy the contents of a matrix
B = A;

// operator *= to multiply by a scalar
B *= 0.8;

// for banded and arrow matrix, operator() is used to know the value
Matrix<double, General, BandedCol> Aband, Bband;

cout << "value of A(2, 0) = " << Aband(2, 0);

// operator = can be used
Bband = Aband;

// and operator *=
Bband *= 2.5;
```

### GetM, GetN

#### Syntax

```  int GetM() const;
int GetN() const;
```

Those methods are identic and return the number of rows/columns of the matrix. Since block-diagonal unsymmetric matrices have a symmetric sparsity pattern, the number of rows is equal to the number of columns. For banded and arrow matrices, the matrices are assumed to be square.

#### Example :

```Matrix<float, General, BlockDiagRow> A;
// GetM() should return 0
cout << "Number of rows of A " << A.GetM() << endl;
```

### GetDataSize

#### Syntax

```  int GetDataSize() const;
```

This method returns the number of elements effectively stored.

#### Example :

```Matrix<float, Symmetric, BlockDiagRowSym> A;
Vector<IVect> BlockNum(2);
BlockNum(0).Reallocate(2); BlockNum(0)(0) = 0; BlockNum(0)(1) = 1;
BlockNum(1).Reallocate(2); BlockNum(1)(0) = 2; BlockNum(1)(1) = 3;
A.SetPattern(BlockNum);
// GetDataSize() should return 6
cout << "Number of elements stored in A " << A.GetDataSize() << endl;

// exists also in banded and arrow matrices
Matrix<double, General, BandedCol> B;

// 5x5 matrix with kl = 1, ku = 2
B.Reallocate(5, 5, 1, 2);

// GetDataSize() should return (2*kl+ku+1)*n = 25
cout << "Number of elements stored in B " << B.GetDataSize() << endl;

```

### GetMemorySize

#### Syntax

```  size_t GetMemorySize() const;
```

This method returns the memory used by the object in bytes.

#### Example :

```Matrix<double, Symmetric, BlockDiagRowSym> A;
Vector<IVect> BlockNum(2);
BlockNum(0).Reallocate(2); BlockNum(0)(0) = 0; BlockNum(0)(1) = 1;
BlockNum(1).Reallocate(2); BlockNum(1)(0) = 2; BlockNum(1)(1) = 3;
A.SetPattern(BlockNum);

cout << "Memory used to store A " << A.GetMemorySize() << endl;

// exists also in banded and arrow matrices
Matrix<double, General, BandedCol> B;

// 5x5 matrix with kl = 1, ku = 2
B.Reallocate(5, 5, 1, 2);

cout << "Memory used to store B " << B.GetMemorySize() << endl;

```

### GetRowSize

#### Syntax

```  int GetRowSize(int i) const;
```

This method returns the number of non-zero entries stored for row i.

#### Example :

```Matrix<double, Symmetric, BlockDiagRowSym> A;
Vector<IVect> BlockNum(2);
BlockNum(0).Reallocate(2); BlockNum(0)(0) = 0; BlockNum(0)(1) = 1;
BlockNum(1).Reallocate(2); BlockNum(1)(0) = 2; BlockNum(1)(1) = 3;
A.SetPattern(BlockNum);

// should return 2
cout << "Size of row 0 : " << B.GetRowSize(0) << endl;

// should return 1 (because of symmetry)
cout << "Size of row 1 : " << B.GetRowSize(1) << endl;

```

### GetSizeMaxBlock

#### Syntax

```  int GetSizeMaxBlock() const;
```

This method returns the maximal size for all the blocks stored in the matrix.

#### Example :

```Matrix<double, Symmetric, BlockDiagRowSym> A;
Vector<IVect> BlockNum(2);
BlockNum(0).Reallocate(3); BlockNum(0)(0) = 0; BlockNum(0)(1) = 1; BlockNum(0)(2) = 4;
BlockNum(1).Reallocate(2); BlockNum(1)(0) = 2; BlockNum(1)(1) = 3;
A.SetPattern(BlockNum);

// should return 3
cout << "Maximal size of a block : " << B.GetSizeMaxBlock() << endl;

```

### GetBlockNumber

#### Syntax

```  int GetBlockNumber(int i) const;
```

This method returns the block number associated with the row i.

#### Example :

```Matrix<double, Symmetric, BlockDiagRowSym> A;
Vector<IVect> BlockNum(2);
BlockNum(0).Reallocate(3); BlockNum(0)(0) = 0; BlockNum(0)(1) = 1; BlockNum(0)(2) = 4;
BlockNum(1).Reallocate(2); BlockNum(1)(0) = 2; BlockNum(1)(1) = 3;
A.SetPattern(BlockNum);

// should return 0
cout << "Block number for row 4 : " << B.GetBlockNumber(4) << endl;

```

### Copy

#### Syntax

```  void Copy(const Matrix& B);
```

This method copies the contents of matrix B in the current matrix. B has the same type as the current type.

#### Example :

```Matrix<double, Symmetric, BlockDiagRowSym> A, B;
Vector<IVect> BlockNum(2);
BlockNum(0).Reallocate(3); BlockNum(0)(0) = 0; BlockNum(0)(1) = 1; BlockNum(0)(2) = 4;
BlockNum(1).Reallocate(2); BlockNum(1)(0) = 2; BlockNum(1)(1) = 3;
A.SetPattern(BlockNum);

// A is filled
A(0, 0) = 2.0;

// then it can be copied into B
B.Copy(A);

```

### GetNbBlocks

#### Syntax

```  int GetNbBlocks() const;
```

This method returns the number of blocks of the matrix.

#### Example :

```Matrix<float, Symmetric, BlockDiagRowSym> A;
Vector<IVect> BlockNum(2);
BlockNum(0).Reallocate(2); BlockNum(0)(0) = 0; BlockNum(0)(1) = 1;
BlockNum(1).Reallocate(2); BlockNum(1)(0) = 2; BlockNum(1)(1) = 3;
A.SetPattern(BlockNum);
// GetNbBlocks() should return 2
cout << "Number of blocks in A " << A.GetNbBlocks() << endl;
```

### GetBlockSize

#### Syntax

```  int GetBlockSize(int n) const;
```

This method returns the size of the block i of the block-diagonal matrix.

#### Example :

```Matrix<float, General, BlockDiagRow> A;
// loop over the blocks
for (int i = 0; i < A.GetNbBlocks(); i++)
cout << "size of the block " << i << " = " << A.GetBlockSize(i) << endl;
```

### Index

#### Syntax

```  int Index(int i, int j) const;
```

This method returns the global row number of the row j of block i.

#### Example :

```Matrix<float, General, BlockDiagRow> A;

// global row numbers are provided through a call to SetPattern
Vector<IVect> num_dofs;
A.SetPattern(num_dofs);

// then when a loop over the blocks is performed
// global row numbers of the block are given by method Index
for (int i = 0; i < A.GetNbBlocks(); i++)
for (int j = 0; j < A.GetBlockSize(i); j++)
cout << "global row number = " << A.Index(i, j) <<endl;

```

### Value

#### Syntax

```  T Value(int i, int j, int k) const;
```

This method returns the value (j, k) of block i of the matrix.

#### Example :

```Matrix<float, General, BlockDiagRow> A;

// global row numbers are provided through a call to SetPattern
Vector<IVect> num_dofs;
A.SetPattern(num_dofs);

// then when a loop over the blocks is performed
// global row numbers of the block are given by method Index
// and values of local matrices are accessed through method Value
for (int i = 0; i < A.GetNbBlocks(); i++)
{
for (int j = 0; j < A.GetBlockSize(i); j++)
cout << "global row number = " << A.Index(i, j) <<endl;

// local matrix associated with block i
int n = A.GetBlockSize(i);
Matrix<float> local_mat(n, n);
for (int j = 0; j < n; j++)
for (int k = 0; k < n; k++)
local_mat(j, k) = A.Value(i, j, k);

// Value can be used to modify the matrix
for (int j = 0; j < n; j++)
for (int k = 0; k < n; k++)
A.Value(i, j, k) = exp(local_mat(j, k));
}
```

### GetData

#### Syntax

```  T* GetData();
```

This method returns the pointer to the array containing all values of non-zero entries. This method is available only for the first implementation of block-diagonal matrices (file MatrixBlockDiagonal.cxx) since in that case the elements of the matrix are stored in a contiguous manner. For banded matrices, the method GetData is also available for the same reasons.

#### Example :

```Matrix<double, General, BlockDiagRow> A;

// global row numbers are provided through a call to SetPattern
Vector<IVect> num_dofs;
A.SetPattern(num_dofs);

// values can be modified directly with the method GetData() (not recommended)
double* data = A.GetData();

data = 1.5;
// ...

// GetData() exists also for banded matrices
Matrix<double, General, BandedCol> B;
B.Reallocate(10, 10, 2, 2);

data = B.GetData();
```

### GetRowPtr, GetPtr

#### Syntax

``` int* GetRowPtr();
int* GetPtr();
```

This method returns the pointer to the array containing the beginning index of each block in array GetInd(). This method is a low level function, that should not be used in a normal use.

#### Example :

```Matrix<double, General, BlockDiagRow> A;

// global row numbers are provided through a call to SetPattern
Vector<IVect> num_dofs;
A.SetPattern(num_dofs);

// global row numbers are stored in an array called ind
int* ind = A.GetInd();
int* ptr = A.GetPtr();

// row numbers of block i are between ptr[i] and ptr[i+1]
int i = 2;
for (int j = ptr[i]; j < ptr[i+1]; j++)
cout  << "row number of block 2 : " << ind[j] << endl;
```

### GetInd

#### Syntax

```  int* GetInd();
```

This method returns the pointer containing row numbers of all the blocks. This method is a low level function, that should not be used in a normal use.

#### Example :

```Matrix<double, General, BlockDiagRow> A;

// global row numbers are provided through a call to SetPattern
Vector<IVect> num_dofs;
A.SetPattern(num_dofs);

// global row numbers are stored in an array called ind
int* ind = A.GetInd();
int* ptr = A.GetPtr();

// row numbers of block i are between ptr[i] and ptr[i+1]
int i = 2;
for (int j = ptr[i]; j < ptr[i+1]; j++)
cout  << "row number of block 2 : " << ind[j] << endl;
```

### GetRev

#### Syntax

```   int* GetRev();
int* GetReverseNum();
```

This method returns the pointer containing local positions of rows in each block. This method is a low level function, that should not be used in a normal use.

#### Example :

```Matrix<double, General, BlockDiagRow> A;

// global row numbers are provided through a call to SetPattern
Vector<IVect> num_dofs;
A.SetPattern(num_dofs);

// global row numbers are stored in an array called ind
int* ind = A.GetInd();

// for a given row j, you can access to the block number where it belongs :
int j = 13;
int* num_blok_ = A.GetNumBlock();
int blok_num = num_blok_[j];

// and the local position of the row inside this block
int* rev = A.GetRev();
int pos = rev[j];
```

### GetNumBlock

#### Syntax

```  int* GetNumBlock();
```

This method returns the pointer containing block numbers for all the rows. This method is a low level function, that should not be used in a normal use.

#### Example :

```Matrix<double, General, BlockDiagRow> A;

// global row numbers are provided through a call to SetPattern
Vector<IVect> num_dofs;
A.SetPattern(num_dofs);

// global row numbers are stored in an array called ind
int* ind = A.GetInd();

// for a given row j, you can access to the block number where it belongs :
int j = 13;
int* num_blok_ = A.GetNumBlock();
int blok_num = num_blok_[j];

// and the local position of the row inside this block
int* rev = A.GetRev();
int pos = rev[j];
```

### GetOffset

#### Syntax

```  int* GetOffset();
```

This method returns the pointer to the array containing offsets in order to access to the data contained in each block. This method is a low level function, that should not be used in a normal use.

#### Example :

```Matrix<double, General, BlockDiagRow> A;

// global row numbers are provided through a call to SetPattern
Vector<IVect> num_dofs;
A.SetPattern(num_dofs);

// values of non-zero entries are stored in a contiguous array data
double* data = A.GetData();

// if you want to modify directly values on a given block, use offset
int* offset = A.GetOffset();

// first value of the block n
int n = 10;
double first_val = data[offset[n]];
```

### Zero

#### Syntax

```  void Zero();
```

This method sets all non-zero entries to 0.

#### Example :

```Matrix<float, Symmetric, BlockDiagRowSym> A;
Vector<IVect> BlockNum(2);
BlockNum(0).Reallocate(2); BlockNum(0)(0) = 0; BlockNum(0)(1) = 1;
BlockNum(1).Reallocate(2); BlockNum(1)(0) = 2; BlockNum(1)(1) = 3;
A.SetPattern(BlockNum);
// you can first fill A with 1, A = [1 1 0 0; 1 1 0 0; 0 0 1 1; 0 0 1 1]
A.Fill(1);
// then you can set all values to 0
A.Zero();
```

### Fill

#### Syntax

```  void Fill(const T& x);
```

This method sets all non-zero entries to a given value.

#### Example :

```Matrix<float, Symmetric, BlockDiagRowSym> A;
Vector<IVect> BlockNum(2);
BlockNum(0).Reallocate(2); BlockNum(0)(0) = 0; BlockNum(0)(1) = 1;
BlockNum(1).Reallocate(2); BlockNum(1)(0) = 2; BlockNum(1)(1) = 3;
A.SetPattern(BlockNum);
// you can first fill A with 1, A = [1 1 0 0; 1 1 0 0; 0 0 1 1; 0 0 1 1]
A.Fill(1);
```

### FillRand

#### Syntax

```  void FillRand();
```

This method sets all non-zero entries to random values.

#### Example :

```Matrix<float, Symmetric, BlockDiagRowSym> A;
Vector<IVect> BlockNum(2);
BlockNum(0).Reallocate(2); BlockNum(0)(0) = 0; BlockNum(0)(1) = 1;
BlockNum(1).Reallocate(2); BlockNum(1)(0) = 2; BlockNum(1)(1) = 3;
A.SetPattern(BlockNum);
// you can first fill A with random values
//  A = [x x 0 0; x x 0 0; 0 0 x x; 0 0 x x]
A.FillRand();
```

### Clear

#### Syntax

```  void Clear();
```

This method clears the matrix. After this call, the matrix is an empty 0x0 matrix.

#### Example :

```Matrix<float, Symmetric, BlockDiagRowSym> A;
// after some modifications on A, you can clear it :
A.Clear();
```

### SetPattern

#### Syntax

```  void SetPattern(const Vector<IVect>& num);
```

This method initializes the sparsity pattern of the matrix. You specify the row numbers for each block.

#### Example :

```Matrix<float, Symmetric, BlockDiagRowSym> A;
// for example you can specify A with three blocks, one with rows 0, 2
// one with rows 1, 3, 5 and one with a single row 4
Vector<IVect> BlockNum(3);
BlockNum(0).Reallocate(2); BlockNum(0)(0) = 0; BlockNum(0)(1) = 2;
BlockNum(1).Reallocate(3); BlockNum(1)(0) = 1; BlockNum(1)(1) = 3; BlockNum(1)(2) = 5;
BlockNum(2).Reallocate(1); BlockNum(2)(0) = 4;
A.SetPattern(BlockNum);
// A is set to 0, you can fill it with AddInteraction
// ...
```

### Reallocate

#### Syntax

```  void Reallocate(int m, int n);

// banded matrix
void Reallocate(int m, int n, int kl, int ku);

// arrow matrix
void Reallocate(int m, int n, int kl, int ku, int nrow, int ncol);
```

This method is not needed for block-diagonal matrices, since the arrays are allocated when the method SetPattern is called. This method is present in order to ensure a compatibility with generic functions. For banded-matrix and arrow matrix, you can provide the number of diagonals in the lower part of A (kl) and the number of upper diagonals in the upper part of A (ku). If not provided, these numbers are assumed to be equal to 0. For an arrow matrix, you can provide the number of dense rows at the end of the matrix (nrow), and the number of dense columns (ncol).

#### Example :

```Matrix<double, Symmetric, BlockDiagRowSym> A;

// you can call Reallocate
A.Reallocate(6, 6);

// but the real allocation is performed with SetPattern
// for example you can specify A with three blocks, one with rows 0, 2
// one with rows 1, 3, 5 and one with a single row 4
Vector<IVect> BlockNum(3);
BlockNum(0).Reallocate(2); BlockNum(0)(0) = 0; BlockNum(0)(1) = 2;
BlockNum(1).Reallocate(3); BlockNum(1)(0) = 1; BlockNum(1)(1) = 3; BlockNum(1)(2) = 5;
BlockNum(2).Reallocate(1); BlockNum(2)(0) = 4;
A.SetPattern(BlockNum);
// A is set to 0, you can fill it with AddInteraction
// ...

// for banded-matrices, the allocation is effective when kl and ku are provided
Matrix<double, General, BandedCol> B;
B.Reallocate(6, 6, 2, 1); // kl = 2, ku = 1

// for arrow matrices, nrow and ncol are denoting the number of last rows and last columns
Matrix<double, General, ArrowCol> C;
C.Reallocate(6, 6, 2, 1, 3, 2); // kl = 2, ku = 1, nrow = 3, ncol = 2

```

### Resize

#### Syntax

```  void Resize(int m, int n);
```

This method is present for compatibility but does nothing. The matrix is allocated/resized by the method SetPattern.

### GetNbRowNumbers

#### Syntax

```  int GetNbRowNumbers() const
```

This method returns the number of rows (including distant rows for distributed block-diagonal matrices). For sequential matrices, the method is equivalent to GetM().

#### Example :

```Matrix<double, Symmetric, BlockDiagRowSym> A;

// for example you can specify A with three blocks, one with rows 0, 2
// one with rows 1, 3, 5 and one with a single row 4
Vector<IVect> BlockNum(3);
BlockNum(0).Reallocate(2); BlockNum(0)(0) = 0; BlockNum(0)(1) = 2;
BlockNum(1).Reallocate(3); BlockNum(1)(0) = 1; BlockNum(1)(1) = 3; BlockNum(1)(2) = 5;
BlockNum(2).Reallocate(1); BlockNum(2)(0) = 4;
A.SetPattern(BlockNum);

// A.GetM() and A.GetNbRowNumbers() should return 6
cout << "Number of rows : " << A.GetNbRowNumbers() << endl;
```

### GetRowNumData

#### Syntax

```  int* GetRowNumData() const
```

This method returns the pointer storing all the row numbers (including distant row numbers for distributed block-diagonal matrices).

#### Example :

```Matrix<double, Symmetric, BlockDiagRowSym> A;

// for example you can specify A with three blocks, one with rows 0, 2
// one with rows 1, 3, 5 and one with a single row 4
Vector<IVect> BlockNum(3);
BlockNum(0).Reallocate(2); BlockNum(0)(0) = 0; BlockNum(0)(1) = 2;
BlockNum(1).Reallocate(3); BlockNum(1)(0) = 1; BlockNum(1)(1) = 3; BlockNum(1)(2) = 5;
BlockNum(2).Reallocate(1); BlockNum(2)(0) = 4;
A.SetPattern(BlockNum);

// retrieving the pointer storing all the row numbers
// it should store here (0, 2, 1, 3, 5, 4)
int* row_num = A.GetRowNumData();
```

#### Syntax

```  void AddInteraction(int i, int j, const T& val);
```

This methods tries to add a value to the entry (i, j) of the matrix. However, if this entry does not belong to the sparsity pattern of the matrix, the value will be ignored for block-diagonal matrices, and the matrix not modified. For banded and arrow matrices, an exception will be raised if the non-zero entry is outside the profile.

#### Example :

```Matrix<double, Symmetric, BlockDiagRowSym> A;
// for example you can specify A with three blocks, one with rows 0, 2
// one with rows 1, 3, 5 and one with a single row 4
Vector<IVect> BlockNum(3);
BlockNum(0).Reallocate(2); BlockNum(0)(0) = 0; BlockNum(0)(1) = 2;
BlockNum(1).Reallocate(3); BlockNum(1)(0) = 1; BlockNum(1)(1) = 3; BlockNum(1)(2) = 5;
BlockNum(2).Reallocate(1); BlockNum(2)(0) = 4;
A.SetPattern(BlockNum);
// A is set to 0, you can fill it with AddInteraction
// for example you set block 0, 2 to [1.5 0.5; 0.5 1.5]

// if you try to add a value to entry (0, 3), it will do nothing

// For banded or arrow matrices, the behaviour is slightly different
Matrix<double, General, ArrowCol> B;

// the profile of the matrix is given by the Reallocate function
int m = 10, n = 10, kl = 1, ku = 2, nrow = 3, ncol = 2;
B.Reallocate(m, n, kl, ku, nrow, ncol);

// this profile cannot be modified by the function AddInteraction
// therefore if the new entry is outside the profile, an exception is launched
// aborting the program
B.AddInteraction(0, 3, 2.3); // the program is aborted
```

#### Syntax

```  void AddInteractionRow(int i, int nb, const IVect& col, const Vector<T>& val);
```

This method tries to add values on the row i of the matrix. However, if some entries does not belong to the sparsity pattern of the matrix, these values will be ignored for block-diagonal matrices. As for AddInteraction, these values will generate an exception for banded and arrow matrices.

#### Example :

```Matrix<float, Symmetric, BlockDiagRowSym> A;
// for example you can specify A with three blocks, one with rows 0, 2
// one with rows 1, 3, 5 and one with a single row 4
Vector<IVect> BlockNum(3);
BlockNum(0).Reallocate(2); BlockNum(0)(0) = 0; BlockNum(0)(1) = 2;
BlockNum(1).Reallocate(3); BlockNum(1)(0) = 1; BlockNum(1)(1) = 3; BlockNum(1)(2) = 5;
BlockNum(2).Reallocate(1); BlockNum(2)(0) = 4;
A.SetPattern(BlockNum);
// A is set to 0, you can fill it with AddInteractionRow
IVect col(3); col(0) = 1; col(1) = 3; col(2) = 5;
Vector<double> val(3); val.Fill();
// here we add 3 values on the row 1
```

### WriteText

#### Syntax

```  void WriteText(const string& file_name) const;
void WriteText(ostream& file_stream) const;
```

This method writes the matrix on a file or on a stream. Row numbers begin with 1 so that Matlab is able to read the file with method spconvert. The format is the same as for sparse matrices (three columns, the first one with row numbers, the second with column numbers and the third one with values). For example :

1 1 1.234
1 2 -0.56
1 5 0.23
2 1 -3.54
...

#### Example :

```Matrix<float, Symmetric, BlockDiagRowSym> A;
A.WriteText("Ah.dat");
```

### Mlt

#### Syntax

```  void Mlt(const T& alpha, Matrix& A);

// for banded or arrow matrices
void Mlt(const Matrix&, const Vector& x, Vector& y);
```

This function multiplies all the values of A by a coefficient alpha. For banded and arrow matrices, this function can also be used to perform a matrix-vector product

#### Example :

```Matrix<double, General, BlockDiagRow> A;
Mlt(2.5, A);

// matrix vector product for banded or arrow matrices
Matrix<double, General, BandedCol> B;
int N = 20, kl = 2, ku = 2;
B.Reallocate(N, N, kl, ku);

Vector<double> x(N), y(N);
x.FillRand();
// computing y = B x
Mlt(B, x, y);
```

#### Syntax

```  void MltAdd(const T& alpha, const Matrix& A, const Vector& X, const T& beta, Vector& B);
```

This function performs a matrix vector product and overwrites B with beta*B + alpha*A*X.

#### Example :

```Matrix<double, General, BlockDiagRow> A;
Vector<double> B(n), X(n); B.Fill(); X.Fill();
// B <- B - A*X
```

#### Syntax

```  void Add(const T& alpha, const Matrix& A, Matrix& B);
```

This function overwrites B with B + alpha*A. If the sparsity pattern of A and B is the same, this operation will directly be performed on values. Otherwise, an intermediary conversion to sparse matrices is used to complete the linear combination. You can't combine symmetric matrix with non-symmetric matrix. For banded and arrow matrices, this operation is performed, assuming that the profile of B contains the profile of A, the profile of B is not modified during the operation. If the profile of A is too large, the program should abort.

#### Example :

```Matrix<double, General, BlockDiagRow> A, B;
// B <- B + 2*A
```

### GetInverse

#### Syntax

```  void GetInverse(Matrix& A);
```

This function overwrites B with its inverse, the operation being completed on each block.

#### Example :

```Matrix<double, General, BlockDiagRow> B;
// B <- inv(B)
GetInverse(B);
```

### Copy

#### Syntax

```  void Copy(const Matrix<T, General, ArrayRowSparse>& A, Matrix& B);
void Copy(const Matrix& A, Matrix<T, General, ArrayRowSparse>& B);
```

This function converts a sparse matrix (only ArrayRowSparse, ArrayRowSymSparse, ArrayRowComplexSparse or ArrayRowSymComplexSparse are accepted) into a block-diagonal (symmetric or not), and vice versa. For banded matrix, a sparse matrix can be converted to a banded matrix (the band width kl and ku are computed), but the reverse operation is not implemented. The conversions procedures are not available for arrow matrices.

#### Example :

```Matrix<double, General, ArrayRowSparse> A;
Matrix<double, Symmetric, BlockDiagRowSym> B;
// you can construct a block diagonal matrix by converting a sparse matrix
Copy(A, B);
// inverse it
GetInverse(B);
// and convert it back to a sparse matrix
Copy(B, A);
```

### ConvertToSparse

#### Syntax

```  void ConvertToSparse(const Matrix& A, Matrix<T, General, ArrayRowSparse>& B);
```

This function converts a block-diagonal matrix to a sparse matrix (ArrayRowSparse or ArrayRowSymSparse for symmetric matrices).

#### Example :

```Matrix<double, Symmetric, ArrayRowSymSparse> A;
Matrix<double, Symmetric, BlockDiagRowSym> B;

// the block-diagonal matrix is filled
// ...

// then you can convert it to a sparse matrix
ConvertToSparse(B, A);
```

### ConvertToBlockDiagonal

#### Syntax

```  void ConvertToBlockDiagonal(const Matrix<T, General, ArrayRowSparse>& B, const Matrix& A);

void ConvertToBlockDiagonal(const Matrix<T, General, ArrayRowSparse>& B, const Matrix& A, int m1, int m2);
```

This function converts a sparse matrix (ArrayRowSparse or ArrayRowSymSparse for symmetric matrices) to block-diagonal matrices (BlockDiagRow or BlockDiagRowSym). The extent of subscripts can be specified with the optional arguments m1 and m2.

#### Example :

```Matrix<double, Symmetric, ArrayRowSymSparse> A;
Matrix<double, Symmetric, BlockDiagRowSym> B;

// the sparse matrix is filled
// ...

// then you can convert it to a block-diagonal matrix
ConvertToSparse(A, B);

// if you want to convert only a submatrix (0, N) x (0, N) :
int N = 10;
ConvertToSparse(A, B, 0, N);
```

### GetCholesky

#### Syntax

```  void GetCholesky(Matrix&)
```

This function replaces a block-diagonal matrix by its Cholesky factor L. This function is relevant only for positive matrices.

#### Example :

```Matrix<double, Symmetric, BlockDiagRowSym> A;

// the matrix is filled
// ...

// then you can compute its Cholesky factorisation
GetCholesky(A);

// you can use this factorisation in order to solve the linear system A x = b
Vector<double> x(A.GetM()), b(A.GetM());
b.FillRand();
x = b;
SolveCholesky(SeldonNoTrans, A, x);
SolveCholesky(SeldonTrans, A, x);

// in order to compute the matrix vector product y = A x
// you cannot use MltAdd since A contains now the Cholesky factor L
// but you can use MltCholesky
Vector<double> y(x);
MltCholesky(SeldonTrans, A, y);
MltCholesky(SeldonNoTrans, A, y);
```

### SolveCholesky

#### Syntax

```  void SolveCholesky(Transpose, Matrix&, Vector&)
```

This function solves the linear system L x = b or LT x = b, the vector contains the source b on input, and the solution x on output. It is assumed that GetCholesky has been previously called in order to compute the Cholesky factor L of the initial matrix A.

#### Example :

```Matrix<double, Symmetric, BlockDiagRowSym> A;

// the matrix is filled
// ...

// then you can compute its Cholesky factorisation
GetCholesky(A);

// you can use this factorisation in order to solve the linear system A x = b
Vector<double> x(A.GetM()), b(A.GetM());
b.FillRand();
x = b;
SolveCholesky(SeldonNoTrans, A, x);
SolveCholesky(SeldonTrans, A, x);

// in order to compute the matrix vector product y = A x
// you cannot use MltAdd since A contains now the Cholesky factor L
// but you can use MltCholesky
Vector<double> y(x);
MltCholesky(SeldonTrans, A, y);
MltCholesky(SeldonNoTrans, A, y);
```

### MltCholesky

#### Syntax

```  void MltCholesky(Transpose, Matrix&, Vector&)
```

This function performs the matrix vector product y = L x or y = LT x, the vector contains x on input, and the result y on output. It is assumed that GetCholesky has been previously called in order to compute the Cholesky factor L of the initial matrix A.

#### Example :

```Matrix<double, Symmetric, BlockDiagRowSym> A;

// the matrix is filled
// ...

// then you can compute its Cholesky factorisation
GetCholesky(A);

// you can use this factorisation in order to solve the linear system A x = b
Vector<double> x(A.GetM()), b(A.GetM());
b.FillRand();
x = b;
SolveCholesky(SeldonNoTrans, A, x);
SolveCholesky(SeldonTrans, A, x);

// in order to compute the matrix vector product y = A x
// you cannot use MltAdd since A contains now the Cholesky factor L
// but you can use MltCholesky
Vector<double> y(x);
MltCholesky(SeldonTrans, A, y);
MltCholesky(SeldonNoTrans, A, y);
```

### GetKL

#### Syntax

```  int GetKL() const
```

This function returns the number of diagonals strictly below the main diagonal. For example, for a tridiagonal matrix, it will return 1.

#### Example :

```Matrix<double, General, BandedCol> A;

// The parameters kl and ku are given when the function Reallocate is called
int m = 10, n = 10, kl = 2, ku = 3;
A.Reallocate(m, n, kl, ku);

for (int d = 0; d < A.GetKL(); d++)
cout << "Sub-diagonal " << d << endl;
```

BandMatrix.cxx

### GetKU

#### Syntax

```  int GetKU() const
```

This function returns the number of diagonals strictly above the main diagonal. For example, for a tridiagonal matrix, it will return 1.

#### Example :

```Matrix<double, General, BandedCol> A;

// The parameters kl and ku are given when the function Reallocate is called
int m = 10, n = 10, kl = 2, ku = 3;
A.Reallocate(m, n, kl, ku);

for (int d = 0; d < A.GetKU(); d++)
cout << "Supra-diagonal " << d << endl;
```

BandMatrix.cxx

### GetNbLastRow

#### Syntax

```  int GetNbLastRow() const
```

This function returns the number of dense rows located at the end of the matrix.

#### Example :

```// given an arrow matrix
Matrix<double, General, ArrowCol> A;

// The number of last rows and columns are given during the allocation
int m = 10, n = 10, kl = 2, ku = 3, nrow = 2, ncol = 3;
A.Reallocate(m, n, kl, ku, nrow, ncol);

// then you can access to the number of dense rows with GetNbLastRow
for (int d = 0; d < A.GetNbLastRow(); d++)
cout << "last row " << d << endl;
```

BandMatrix.cxx

### GetNbLastCol

#### Syntax

```  int GetNbLastCol() const
```

This function returns the number of dense columns located at the end of the matrix.

#### Example :

```// given an arrow matrix
Matrix<double, General, ArrowCol> A;

// The number of last rows and columns are given during the allocation
int m = 10, n = 10, kl = 2, ku = 3, nrow = 2, ncol = 3;
A.Reallocate(m, n, kl, ku, nrow, ncol);

// then you can access to the number of dense columns with GetNbLastCol
for (int d = 0; d < A.GetNbLastCol(); d++)
cout << "last column " << d << endl;
```

BandMatrix.cxx

### ClearRow

#### Syntax

```  void ClearRow(int i)
```

This function sets to 0 all non-zero entries of the row i.

#### Example :

```// given an arrow matrix
Matrix<double, General, ArrowCol> A;

// The number of last rows and columns are given during the allocation
int m = 10, n = 10, kl = 2, ku = 3, nrow = 2, ncol = 3;
A.Reallocate(m, n, kl, ku, nrow, ncol);

// then the profile is no longer modified
// ClearRow merely sets to 0 the values of non-zero entries of a row
A.ClearRow(1); // here the row 1
```

### ClearColumn

#### Syntax

```  void ClearColumn(int i)
```

This function sets to 0 all non-zero entries of the column i.

#### Example :

```// given a block-diagonal matrix
Matrix<double, General, BlockDiagRow> A;

// then the profile is no longer modified
// ClearColumn merely sets to 0 the values of non-zero entries of a column
A.ClearColumn(1); // here the row 1
```

### Val

#### Syntax

```  void Val(int i, int j)
```

This function gives an access to the element (i, j) of a banded or arrow matrix, whereas the operator () returns the value of A(i, j) (the matrix cannot be modified).

#### Example :

```// given an arrow matrix
Matrix<double, General, ArrowCol> A;

// The number of last rows and columns are given during the allocation
int m = 10, n = 10, kl = 2, ku = 3, nrow = 2, ncol = 3;
A.Reallocate(m, n, kl, ku, nrow, ncol);

// then the values can be modified with Val
int i = 1, j = 2;
A.Val(i, j) = 2.3;
```

BandMatrix.cxx

### Get

#### Syntax

```  void Get(int i, int j)
```

This function gives an access to the element (i, j) of a banded or arrow matrix, whereas the operator () returns the value of A(i, j) (the matrix cannot be modified).

#### Example :

```// given an arrow matrix
Matrix<double, General, ArrowCol> A;

// The number of last rows and columns are given during the allocation
int m = 10, n = 10, kl = 2, ku = 3, nrow = 2, ncol = 3;
A.Reallocate(m, n, kl, ku, nrow, ncol);

// then the values can be modified with Get
int i = 1, j = 2;
A.Get(i, j) = 2.3;
```

BandMatrix.cxx

### Set

#### Syntax

```  void Set(int i, int j, double x)
```

This function modifies the element (i, j) of a banded or arrow matrix. If the element (i, j) does not belong to the sparsity pattern of the matrix, the program should be aborted.

#### Example :

```// given an arrow matrix
Matrix<double, General, ArrowCol> A;

// The number of last rows and columns are given during the allocation
int m = 10, n = 10, kl = 2, ku = 3, nrow = 2, ncol = 3;
A.Reallocate(m, n, kl, ku, nrow, ncol);

// then the values can be set with the function Set
int i = 1, j = 2;
A.Set(i, j, 2.3);
```

BandMatrix.cxx

### SetIdentity

#### Syntax

```  void SetIdentity()
```

This function fills the matrix such that it is equal to the identity matrix.

#### Example :

```// given an arrow matrix
Matrix<double, General, ArrowCol> A;

// The number of last rows and columns are given during the allocation
int m = 10, n = 10, kl = 2, ku = 3, nrow = 2, ncol = 3;
A.Reallocate(m, n, kl, ku, nrow, ncol);

// then the values can be initialized to 0 except 1 on the diagonal :
A.SetIdentity();
```

BandMatrix.cxx

### Factorize

#### Syntax

```   void Factorize();
void Factorize(IVect& pivot);
```

This function computes the LU factorisation of a banded or arrow matrix without pivoting such that no fill-in occurs. The matrix is replaced by its factors L and U. Lapack is here not called. The second syntax performs a LU factorization with pivoting (available only for banded matrices).

#### Example :

```  // given an arrow matrix
Matrix<double, General, ArrowCol> A;

// the matrix is filled
// ...

// then you can factorize the matrix
A.Factorize();

// and solve each linear system A x = b with Solve
Vector<double> x(A.GetM()), b(A.GetM());
b.FillRand();
x = b;
A.Solve(x);

```

BandMatrix.cxx

### Solve

#### Syntax

```  void Solve(Vector&)
void Solve(const IVect& pivot, Vector&)
```

This function solves the linear system A x = b, assuming that the LU factorisation has been computed by calling the function Factorize. The second syntax solves the linear system in the case where factorization with pivoting has been called (available only for banded matrices).

#### Example :

```// given an arrow matrix
Matrix<double, General, ArrowCol> A;

// the matrix is filled
// ...

// then you can factorize the matrix
A.Factorize();

// and solve each linear system A x = b with Solve
Vector<double> x(A.GetM()), b(A.GetM());
b.FillRand();
x = b;
A.Solve(x);

```

BandMatrix.cxx

### Write

#### Syntax

```  void Write(const string& file_name) const;
void Write(ostream& file_stream) const;
```

This method writes the matrix on a file or on a stream. The datas are written in binary.

#### Example :

```Matrix<float, General, BandedCol> A;
A.Write("Ah.dat");
```

BandMatrix.cxx

### GetLU

#### Syntax

``` void GetLU(Matrix&);
void GetLU(Matrix& A, Matrix& mat_lu);
void GetLU(Matrix& A, Matrix& mat_lu, bool keep_matrix);
void GetLU(Matrix& A, IVect& pivot);
```

The first syntax computes the LU factorisation of a banded/arrow matrix, the matrix being overwritten by factors L and U. The second syntax write the factors L and U in mat_lu, the matrix A being cleared. The third syntax keeps the initial matrix A if keep_matrix is set to true. The last syntax is available for banded matrices only, and uses Lapack subroutines to perform the LU factorisation. In that case the factorisation is completed with pivoting.

#### Example :

```// factorisation of A without pivoting
// A is replaced by its factors L and U
Matrix<double, General, ArrowCol> A;
GetLU(A);

// SolveLU is called to obtain the solution to A x = b
Vector<double> b(A.GetM()), x;
b.FillRand();
x = b;
SolveLU(A, x);

// if you want to keep the initial matrix :
Matrix<double, General, ArrowCol> B, mat_lu;
GetLU(B, mat_lu, true);
// B contains the initial matrix, mat_lu the factors L and U

// SolveLU again to obtain the final solution
SolveLU(mat_lu, x);

// if you desire to call Lapack subroutine, provide a pivot array
// this possibility is present for banded matrices only
Matrix<double, General, BandedCol> C;
// the Lapack subroutine performs a factorisation with pivoting
Vector<int> pivot;
GetLU(A, pivot);

// SolveLU to obtain the final solution
SolveLU(A, pivot, x);
```

BandMatrix.cxx

### SolveLU

#### Syntax

``` void SolveLU(Matrix&, Vector&);
void SolveLU(Matrix& A, IVect& pivot, Vector& x);
```

When a pivot array is provided, Lapack subroutine is called, whereas in other cases the function Solve is used. The resolution with pivoting can be completed if a factorisation with pivoting has been previously called.

#### Example :

```// factorisation of A without pivoting
// A is replaced by its factors L and U
Matrix<double, General, ArrowCol> A;
GetLU(A);

// SolveLU is called to obtain the solution to A x = b
Vector<double> b(A.GetM()), x;
b.FillRand();
x = b;
SolveLU(A, x);

// if you want to keep the initial matrix :
Matrix<double, General, ArrowCol> B, mat_lu;
GetLU(B, mat_lu, true);
// B contains the initial matrix, mat_lu the factors L and U

// SolveLU again to obtain the final solution
SolveLU(mat_lu, x);

// if you desire to call Lapack subroutine, provide a pivot array
// this possibility is present for banded matrices only
Matrix<double, General, BandedCol> C;
// the Lapack subroutine performs a factorisation with pivoting
Vector<int> pivot;
GetLU(A, pivot);

// SolveLU to obtain the final solution
SolveLU(A, pivot, x);
```

BandMatrix.cxx