BLAS Functions

Those functions are available both for dense and sparse vectors. In the case of dense vectors, Blas subroutines are called if SELDON_WITH_BLAS (or, for backward compatibility, SELDON_WITH_CBLAS) is defined.

Mlt multiplies the elements of the vector/matrix by a scalar, or performs a matrix-vector product
MltAdd performs a matrix-vector product or a matrix-matrix product
Add adds two vectors or two matrices
Copy copies a vector into another one, or a matrix into another matrix
Swap exchanges two vectors
SwapPointer exchanges two vectors (only adresses are exchanged)
GenRot computes givens rotation
GenModifRot computes givens roation
ApplyRot applies givens rotation
ApplyModifRot applies givens rotation
DotProd scalar product between two vectors
DotProdConj scalar product between two vectors, first vector being conjugated
Conjugate conjugates all elements of a vector
GetMaxAbsIndex returns index where highest absolute value is reached
Norm1 returns 1-norm of a vector
Norm2 returns 2-norm of a vector
Rank1Update Adds a contribution X.Y' to a matrix
Rank2Update Adds a contribution X.Y' + Y.X' to a symmetric matrix
Solve solves a triangular system

Mlt

Syntax :

  void Mlt(T0, Vector&);
  void Mlt(T0, Matrix&);
  void Mlt(T0, Array3D&);  

void Mlt(const Matrix&, const Vector&, Vector&); void Mlt(SeldonTrans, const Matrix&, const Matrix&, Matrix&); void Mlt(const Matrix&, const Matrix&, Matrix&);

This function multiplies all elements of vector/matrix/array by a constant. It can also perform a matrix-vector product. This function is implemented both for dense and sparse vectors/matrices.

Example :

Vector<double> X(3);
X.Fill();
// we multiply all elements by 2
Mlt(2.0, X);

// similar method for matrix
Matrix<double> A(3, 3);
A.Fill();
// we multiply all elements by -3.4
Mlt(-3.4, A);

// similar stuff for a 3-D array
Array3D<double> C(3, 4, 6);
C.Fill();
Mlt(-2.5, C);

// and sparse matrix as well
Matrix<double, General, ArrayRowSparse> Asp(3,3), Bsp(3, 3), Csp(3, 3);
Asp(0, 1) = 1.4;
Asp(1, 1) = -0.3; 
Asp(2, 0) = 2.1;
Mlt(2.1, Asp);

// with Mlt, you can compute Y = Asp*X or Y = A*X
Vector<double> Y(3);
Mlt(A, X, Y);
Mlt(Asp, X, Y);

// matrix-matrix product
// available for dense and sparse matrices
Mlt(Asp, Bsp, Csp);

// you can multiply with the transpose matrix
// Y = A^T x
Mlt(SeldonTrans, A, X, Y);

// for triangular matrices, the result is overwritten in X
// X = T X
Matrix<double, General, RowUpTriangPacked> T(3, 3);
Mlt(T, X);

// multiplication with transpose
Mlt(SeldonTrans, SeldonNonUnit, T, X);

// multiplication with T, assuming that diagonal of T is equal to 1
X = transpose(T) X
Mlt(SeldonNoTrans, SeldonUnit, T, X);

// multiplication triangular/ dense matrices
// M = T M
Matrix<double> M(3, 5);
Mlt(T, M);

Related topics :

Vector operators
Matrix operators
MltAdd

Location :

Functions_Vector.cxx
Functions_MatVect.cxx
Functions_Matrix.cxx
Functions_MatrixArray.cxx
Blas_1.cxx
Blas_2.cxx
Blas_3.cxx

MltAdd

Syntax :

  void MltAdd(T0, const Matrix&, const Vector&, T0, Vector&);
  void MltAdd(T0, SeldonTrans, const Matrix&, const Vector&, T0, Vector&);
  void MltAdd(T0, const Matrix&, const Matrix&, T0, Matrix&);
  void MltAdd(T0, SeldonTrans, const Matrix&, SeldonTrans, const Matrix&,
              T0, Matrix&);

This function performs a matrix-vector product or a matrix-matrix product. You can specify a product with the transpose of the matrix, or the conjugate of the transpose. The size of matrices need to be compatible in order to complete the matrix-matrix product.

Example :

// matrix-vector product
Matrix<double> A(3,3);
Vector<double> X(3), B(3);
A.Fill();
X.Fill();
B.Fill(1.5);
double beta = 2, alpha = 3;
// computation of B = beta*B + alpha*A*X
MltAdd(alpha, A, X, beta, B);

// same function for sparse matrices
Matrix<double, General, RowSparse> Asp(3, 3, nnz, values, ptr, ind);
// computation of B = beta*B + alpha*Asp*X
MltAdd(alpha, Asp, X, beta, B);

// you can compute B = beta*B + alpha*transpose(Asp)*X
MltAdd(alpha, SeldonTrans, Asp, X, beta, B);

// similar method for matrix-matrix product
// available for sparse and dense matrices
Matrix<double> M(3,4), N(4,3);
// computation of A = beta*A + alpha*M*N
MltAdd(alpha, M, N, beta, A);

// and you can use SeldonTrans, SeldonNoTrans, SeldonConjTrans

// computation of A = beta*A + alpha*transpose(M)*N
MltAdd(alpha, SeldonTrans, M, SeldonNoTrans, N, beta, A);

// SeldonConjTrans is meaningfull for complex matrices
Matrix<complex<double> > Ac(3, 3), Bc(3, 3), Cc(3, 3);

// computation of Ac = beta * Ac
//                     + alpha * conjugate(transpose(Bc)) * transpose(Cc)
MltAdd(alpha, SeldonConjTrans, Bc, SeldonTrans, Cc, beta, Ac);

Related topics :

Vector operators
Matrix operators
Mlt

Location :

Functions_MatVect.cxx
Functions_Matrix.cxx
Functions_MatrixArray.cxx
Blas_2.cxx
Blas_3.cxx

Add

Syntax :

  void Add(T0, const Vector&, Vector&);
  void Add(T0, const Matrix&, Matrix&);

This function adds two vectors or two matrices. This function is available both for sparse or dense vectors/matrices. The size of the vectors or the matrices to add should be the same.

Example :

Vector<double> X(3), Y(3);
X.Fill(); Y.Fill(1);
double alpha = 2;
// computation of Y = Y + alpha*X
Add(alpha, X, Y);

// similar method for matrix
Matrix<double> A(3,3), B(3,3);
A.Fill();
B.SetIdentity();
Add(alpha, A, B);

// and sparse matrix as well
Matrix<double, General, ArrayRowSparse> Asp(3,3), Bsp(3,3);
Asp(0, 1) = 1.4;
Asp(1, 1) = -0.3; 
Asp(2, 0) = 2.1;
Bsp.SetIdentity();
Add(alpha, Asp, Bsp);

Location :

Functions_Vector.cxx
Functions_Matrix.cxx
Functions_MatrixArray.cxx
Blas_1.cxx

Copy

Syntax :

  void Copy(const Vector&, Vector&);
  void Copy(const Matrix& Matrix&);

A vector is copied into another one. If BLAS is enabled, the two dense vectors have to be of the same size. The operator = is more flexible. For sparse matrices, the function Copy is overloaded to enable the conversion between the different formats.

Example :

// vectors need to have the same size
Vector<float> V(3), W(3);
V.Fill();
Copy(V, W);
// it is equivalent to write W = V;
W = V;

// for sparse matrices, you can use Copy as a convert tool
Matrix<double, General, ArrayRowSparse> A(3, 3);
A.AddInteraction(0, 2, 2.5);
A.AddInteraction(1, 1, -1.0);

Matrix<double, General, RowSparse> B;
Copy(A, B);

Related topics :

Vector operators

Location :

Functions_Vector.cxx
Blas_1.cxx
Matrix_Conversions.cxx

Swap

Syntax :

  void Swap(Vector& Vector&);

The two vectors are exchanged. If BLAS is enabled, the two dense vectors have to be of the same size.

Example :

// vectors need to have the same size
Vector<float> V(3), W(3);
V.Fill(); 
W.Fill(1.0);
Swap(V, W);

Related topics :

Vector operators

SwapPointer

Syntax :

  void SwapPointer(Vector& Vector&);

The two vectors are exchanged (only pointers are exchanged).

Example :

// vectors do not need to have the same size
Vector<float> V(3), W(5);
V.Fill(); 
W.Fill(1.0);
Swap(V, W);

Related topics :

Vector operators

Location :

Functions_Vector.cxx
Blas_1.cxx

GenRot, ApplyRot, GenModifRot, ApplyModifRot

Syntax :

  void GenRot(T&, T&, T&, T&);
  void ApplyRot(Vector&, Vector&, T, T); 
  void GenModifRot(T&, T&, T&, T&, T*);
  void ApplyModifRot(Vector&, Vector&, T*); 

This function constructs the Givens rotation G = [cos_, sin_; -sin_, cos_] that zeros the second entry of the two vector (a,b).

Example :

double a, b, cos_, sin_;
a = 2.0;
b = 1.0;
// we compute cos_ and sin_ so that
// G = [cos_, sin_; -sin_, cos_] zeros second entry of vector [a;b]
// i.e. G*[a;b] = [r;0]
GenRot(a, b, cos_, sin_);

// then we can apply rotation G to vectors X, Y
// G*[x_i, y_i] for all i
Vector<double> X(10), Y(10);
X.FillRand(); 
Y.FillRand();
ApplyRot(X, Y, cos_, sin_);

Location :

Functions_Vector.cxx
Blas_1.cxx

DotProd, DotProdConj

Syntax :

  T DotProd(const Vector&, const Vector&);
  T DotProdConj(const Vector&, const Vector&);

This function returns the scalar product between two vectors. For DotProdConj, the first vector is conjugated. In the case of dense vectors, they have to be of the same size.

Example :

// we construct two vectors of same size
Vector<double> X(10), Y(10);
X.FillRand(); 
Y.FillRand();

// computation of X*Y
double scal = DotProd(X, Y);

// for complex numbers
Vector<complex<double> > Xc(10), Yc(10);
Xc.Fill(); Yc.Fill();
// computation of conj(X)*Y
complex<double> scalc = DotProdConj(X, Y);

Location :

Functions_Vector.cxx
Blas_1.cxx

Conjugate

Syntax :

  void Conjugate(Vector&);
  void Conjugate(Matrix&);

Each component of the vector (or matrix) is conjugated. In the case of real numbers, the vector (or matrix) is not modified.

Example :

// complex vector
Vector<complex<double> > X(10);
X.FillRand(); 

// we take the conjugate
Conjugate(X);

// complex matrix
Matrix<complex<double> > A(10, 8);
A.FillRand(); 

// we take the conjugate
Conjugate(A);

Location :

Functions_Vector.cxx
Blas_1.cxx

GetMaxAbsIndex

Syntax :

   int GetMaxAbsIndex(const Vector&);

This function returns the index for which the vector reached its highest absolute value.

Example :

// complex vector
Vector<complex<double> > X(10);
X.FillRand(); 

int imax = GetMaxAbsIndex(X);
// maximum ?
double maximum = abs(X(imax));

Location :

Functions_Vector.cxx
Blas_1.cxx

Norm1

Syntax :

   T Norm1(const Vector&);

This function returns the 1-norm of the vector.

Example :

// complex vector
Vector<complex<double> > X(10);
X.Fill(); 

// 1-norm
double norme = Norm1(X);

Location :

Functions_Vector.cxx
Blas_1.cxx

Norm2

Syntax :

   T Norm2(const Vector&);

This function returns the 2-norm of the vector.

Example :

// complex vector
Vector<complex<double> > X(10);
X.Fill(); 

// 2-norm
double norme = Norm2(X);

Location :

Functions_Vector.cxx
Blas_1.cxx

Rank1Update

Syntax :

  void Rank1Update(T, const Vector&, const Vector&, Matrix&);
  void Rank1Update(T, const Vector&, SeldonConj, const Vector&, Matrix&);
  void Rank1Update(T, const Vector&, Matrix&);

This function adds a contribution of rank 1 to matrix.

Example :

// two complex vectors
Vector<complex<double> > X(10), Y(10);
X.FillRand(); 
Y.FillRand(); 

// complex matrix
Matrix<complex<double> > A(10, 10);
A.SetIdentity();
complex<double> alpha(1, 1);

// we compute A = A + alpha*X*transpose(Y)
Rank1Update(alpha, X, Y, A);
// we can also compute A = A + alpha*X*conjugate(transpose(Y))
Rank1Update(alpha, X, SeldonConj, Y, A);

// for hermitian matrices, use of X only
Matrix<complex<double>, General, RowHermPacked> B(10, 10);
// we compute B = B + alpha*X*conjugate(transpose(X))
Rank1Update(alpha, X, B);

// for real numbers, use of symmetric matrices
Vector<double> Xr(10);
Xr.FillRand(); 
Matrix<double, Symmetric, RowSymPacked> Br(10, 10);
Br.Zero();
double beta = 2;
// computation of Br = Br + beta*Xr*tranpose(Xr)
Rank1Update(beta, Xr, Br);

Location :

Blas_2.cxx

Rank2Update

Syntax :

  void Rank2Update(T, const Vector&, const Vector&, Matrix&);

This function adds a contribution of rank 2 to a symmetric or Hermitian matrix.

Example :

// two complex vectors
Vector<complex<double> > X(10), Y(10);
X.FillRand(); 
Y.FillRand(); 

// hermitian matrix
Matrix<complex<double>, General, RowHermPacked> A(10, 10);
A.SetIdentity();
complex<double> alpha(1, 1);

// we compute A = A + alpha * (X * conjugate(transpose(Y))
//                + Y * conjugate(transpose(X)))
Rank2Update(alpha, X, Y, A);

// for real numbers, use of symmetric matrices
Vector<double> Xr(10), Yr(10);
Xr.FillRand(); 
Yr.FillRand(); 
Matrix<double, Symmetric, RowSymPacked> Br(10, 10);
Br.Zero();
double beta = 2;
// computation of Br = Br + beta*(Xr*transpose(Yr) + Yr*transpose(Xr))
Rank2Update(beta, Xr, Yr, Br);

Location :

Blas_2.cxx

Solve

Syntax :

  void Solve(const Matrix&, Vector&);
  void Solve(SeldonTrans, SeldonNonUnit, const Matrix&, Vector&);

This function solves a triangular linear system.

Example :

// right hand side
Vector<double> B(3);
B.Fill(); 

// triangular matrix
Matrix<double, General, RowUpTriangPacked> T(3, 3);
T.Fill();

// now you can solve T*X = B, the result overwrites B
Solve(T, B);

// if you want to solve transpose(T)*X = B
Solve(SeldonTrans, SeldonNonUnit, T, B);

// SeldonUnit can be used if you are assuming that T has 
// a unitary diagonal (with 1)

// we force unitary diagonal
for (int i = 0; i < T.GetM(); i++)
  T(i, i) = 1.0;

// then we can call Solve with SeldonUnit to solve T*X = B
// B can be also a dense matrix
Solve(SeldonNoTrans, SeldonUnit, T, B);

Location :

Blas_2.cxx
Blas_3.cxx