Block-diagonal and banded matrices

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;

Location :

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;

Location :

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

Location :

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;

Location :

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;

Location :

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;

Location :

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));
  }

Location :

MatrixBlockDiagonal.cxx
BandMatrix.cxx

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[0] = 1.5;
// ...


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

data = B.GetData();

Location :

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;

Location :

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];

Location :

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];

Location :

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]];

Location :

MatrixBlockDiagonal.cxx
OptMatrixBlockDiagonal.cxx
BandMatrix.cxx

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();

Location :

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

Location :

MatrixBlockDiagonal.cxx
OptMatrixBlockDiagonal.cxx
BandMatrix.cxx

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();

Location :

MatrixBlockDiagonal.cxx
OptMatrixBlockDiagonal.cxx
BandMatrix.cxx

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();

Location :

MatrixBlockDiagonal.cxx
OptMatrixBlockDiagonal.cxx
BandMatrix.cxx

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
A.AddInteraction(0, 2, 2.34);
// ...

Location :

MatrixBlockDiagonal.cxx
OptMatrixBlockDiagonal.cxx
BandMatrix.cxx

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
A.AddInteraction(0, 2, 2.34);
// ...

// 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

Location :

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.

Location :

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;

Location :

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();

Location :

MatrixBlockDiagonal.cxx
OptMatrixBlockDiagonal.cxx
BandMatrix.cxx

AddInteraction

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]
A.AddInteraction(0, 0, 1.5);
A.AddInteraction(0, 2, 0.5);
A.AddInteraction(2, 0, 0.5);
A.AddInteraction(2, 2, 1.5);

// if you try to add a value to entry (0, 3), it will do nothing
A.AddInteraction(0, 3, 1.5);

// 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

Location :

AddInteractionRow

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
A.AddInteractionRow(1, 3, col, val);

Location :

MatrixBlockDiagonal.cxx
OptMatrixBlockDiagonal.cxx
BandMatrix.cxx

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");

Location :

MatrixBlockDiagonal.cxx
OptMatrixBlockDiagonal.cxx
BandMatrix.cxx

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

Location :

MatrixBlockDiagonal.cxx
OptMatrixBlockDiagonal.cxx
BandMatrix.cxx

MltAdd

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
MltAdd(2.0, A, X, 1.0, B);

Location :

MatrixBlockDiagonal.cxx
OptMatrixBlockDiagonal.cxx
BandMatrix.cxx

Add

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
Add(2.0, A, B);

Location :

MatrixBlockDiagonal.cxx
OptMatrixBlockDiagonal.cxx
BandMatrix.cxx

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

Location :

MatrixBlockDiagonal.cxx
OptMatrixBlockDiagonal.cxx
BandMatrix.cxx

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

Location :

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

Location :

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

Location :

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

Location :

SolveCholesky

Syntax

  void SolveCholesky(Transpose, Matrix&, Vector&)

This function solves the linear system L x = b or L^T 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);

Location :

MltCholesky

Syntax

  void MltCholesky(Transpose, Matrix&, Vector&)

This function performs the matrix vector product y = L x or y = L^T 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);

Location :

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

// then you can access to kl with GetKL
for (int d = 0; d < A.GetKL(); d++)
  cout << "Sub-diagonal " << d << endl;

Location :

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

// then you can access to ku with GetKU
for (int d = 0; d < A.GetKU(); d++)
  cout << "Supra-diagonal " << d << endl;

Location :

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;

Location :

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;

Location :

MatrixBlockDiagonal.cxx 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

Location :

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

Location :

MatrixBlockDiagonal.cxx BandMatrix.cxx

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;

Location :

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;

Location :

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

Location :

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();

Location :

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

Location :

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

Location :

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");

Location :

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

Location :

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

Location :