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 |
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 |
void Mlt(T0, Vector&); void Mlt(T0, Matrix&); void Mlt(Matrix&, const Vector&, Vector&);
This function multiplies all elements of vector/matrix by a constant. It can also perform a matrix-vector product. This function is implemented both for dense and sparse vectors/matrices.
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);
// and sparse matrix as well
Matrix<double, General, ArrayRowSparse> Asp(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);
Vector operators
Matrix operators
MltAdd
Functions_Vector.cxx
Functions_MatVect.cxx
Functions_MatrixArray.cxx
Blas_1.cxx
Blas_2.cxx
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. The matrix-matrix product is not available for sparse matrices.
// 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
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);
Vector operators
Matrix operators
Mlt
Functions_MatVect.cxx
Functions_Matrix.cxx
Functions_MatrixArray.cxx
Blas_2.cxx
Blas_3.cxx
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.
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);
Functions_Vector.cxx
Functions_Matrix.cxx
Functions_MatrixArray.cxx
Blas_1.cxx
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.
// 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);
Functions_Vector.cxx
Blas_1.cxx
Matrix_Conversions.cxx
void Swap(Vector& Vector&);
The two vectors are exchanged. If BLAS is enabled, the two dense vectors have to be of the same size.
// vectors need to have the same size
Vector<float> V(3), W(3);
V.Fill();
W.Fill(1.0);
Swap(V, W);
Functions_Vector.cxx
Blas_1.cxx
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).
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_);
Functions_Vector.cxx
Blas_1.cxx
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.
// 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);
Functions_Vector.cxx
Blas_1.cxx
void Conjugate(Vector&);
Each component of the vector is conjugated. In the case of real numbers, the vector is not modified.
// complex vector
Vector<complex<double> > X(10);
X.FillRand();
// we take the conjugate
Conjugate(X);
Functions_Vector.cxx
Blas_1.cxx
int GetMaxAbsIndex(const Vector&);
This function returns the index for which the vector reached its highest absolute value.
// complex vector
Vector<complex<double> > X(10);
X.FillRand();
int imax = GetMaxAbsIndex(X);
// maximum ?
double maximum = abs(X(imax));
Functions_Vector.cxx
Blas_1.cxx
T Norm1(const Vector&);
This function returns the 1-norm of the vector.
// complex vector
Vector<complex<double> > X(10);
X.Fill();
// 1-norm
double norme = Norm1(X);
Functions_Vector.cxx
Blas_1.cxx
T Norm2(const Vector&);
This function returns the 2-norm of the vector.
// complex vector
Vector<complex<double> > X(10);
X.Fill();
// 2-norm
double norme = Norm2(X);
Functions_Vector.cxx
Blas_1.cxx
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.
// 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);
void Rank2Update(T, const Vector&, const Vector&, Matrix&);
This function adds a contribution of rank 2 to a symmetric or hermitian matrix.
// 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);
void Solve(const Matrix&, Vector&); void Solve(SeldonTrans, SeldonNonUnit, const Matrix&, Vector&);
This function solves a triangular linear system.
// 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
Solve(SeldonNoTrans, SeldonUnit, T, B);