Computations

Seldon is fully interfaced with Blas (level 1, 2 and 3), except for Blas functions involving banded matrices and rank operations. If Blas is not available to the user, a few alternative functions (same functions written in C++) may be used.

Blas

The interface is implemented in the files Seldon/computation/interfaces/Blas_*) if you have a doubt about the syntax. In the following list, C++ functions are detailed, and the equivalent BLAS subroutines mentionned in bold. Duplicate names are avoided, e.g. only daxpy is mentionned even if other subroutines saxpy, caxpy and zaxpy are interfaced as well.

  • GenRot(a, b, c, s) constructs givens plane rotation (drotg)

  • ApplyRot(X, Y, c, s) applies givens plane rotation (drot)

  • GenModifRot(a, b, c, s, param) constructs modified givens plane rotation (drotmg)

  • ApplyModifRot(X, Y, param) applies modified givens plane rotation (drotm)

  • Swap(X, Y) exchanges two vectors (dswap).

  • Copy(X, Y) copies vector X into vector Y (dcopy).

  • Mlt(alpha, X) multiplies vector X by alpha (dscal).

  • Add(alpha, X, Y) adds vector X to Y : alpha.X + Y -> Y (daxpy).

  • DotProd(X, Y) returns scalar product between two vectors (ddot, zdotu).

  • DotProdConj(X, Y) returns scalar product between two vectors, first vector being conjugated (zdotc).

  • Norm1(X) computes 1-norm of vector X (dasum).

  • Norm2(X) computes 2-norm of vector X (dnrm2).

  • GetMaxAbsIndex(X) returns index where maximum of X is reached (idamax).

  • Mlt(A, X) computes product of A by X for triangular matrices : A*X -> X (dtrmv, dtpmv).

  • MltAdd(alpha, A, X, beta, Y) computes matrix-vector product alpha*A*X + beta*Y -> Y. This function is overloaded for all types of matrices (dgemv, zhemv, zhpmv, dsymv, dspmv).

  • Rank1Update(alpha, X, Y, M) Performs rank 1 operation : alpha*X.Y' + M -> M (dger, zgeru)

  • Rank1Update(alpha, X, conj, Y, M) Performs rank 1 operation : alpha*X.Y' + M -> M by taking the conjugate of Y (zgerc)

  • Rank1Update(alpha, X, M) Performs rank 1 operation : alpha*X.X' + M -> M for symmetric/hermitian matrices (dspr, zhpr)

  • Rank1Update(alpha, X, M) Performs rank 2 operation : alpha*X.Y' + Y.X' + M -> M for symmetric/hermitian matrices (dspr2, zhpr2)

  • Solve(A, X) Solves triangular linear system A\X -> X (dtrsv, dtpsv)

  • MltAdd(alpha, A, B, beta, C) computes matrix-matrix product alpha*A*B + beta*C -> C. This function is overloaded for all types of matrices (dgemm, dsymm, zhemm).

  • Mlt(alpha, A, B) computes matrix-matrix product alpha*A*B -> B with A triangular matrix (dtrmm, dsymm, zhemm).

  • Solve(alpha, A, B) Solves triangular linear system alpha*A\B -> B with B general matrix (dtrsm, dtpsv, ztpsv )

Remember that, before the inclusion of Seldon (#include "Seldon.hxx"), you must define SELDON_WITH_BLAS. For backward compatibility, SELDON_WITH_CBLAS is also accepted, but not recommended. If you use a Makefile (such as the file Makefile.LINUX), the flag is activated if USE_BLAS is set to YES. Below is an example: 

#define SELDON_WITH_BLAS
#include "Seldon.hxx"

using namespace Seldon;

// complex vectors
int n = 10;
ZVect X(n), Y(n); 
X. Fill(); Y.Fill();

// scalar product X'.Y (conjugate of X)
complex<double> alpha = DotProdConj(X, Y);
// norm of the vector
double norm_vec = Norm2(X);
// alpha*X + Y -> Y
Add(alpha, X, Y);

// matrices must have the good dimension before multiplication
int n = 4, k = 5;
Matrix<double> A(m,n), B(n,k), C(m,k);
MltAdd(1.0, A, B, 2.0, C);

// you can multiply with transpose
A.Reallocate(n,m); A.Fill();
MltAdd(1.0, SeldonTrans, A, SeldonNoTrans, B, 2.0, C);

A comprehensive test of Blas interface is done in file test/program/blas_test.cpp.

C++ functions

C++ functions (that do not call Blas) are available in Seldon/computation/basic_functions/ *. The syntax is the same as for the functions in the interface to Blas. The functions Add, DotProd , DotProdConj , Mlt and MltAdd are written in C++ for any type of matrix and vector, and those functions can be used without using the Blas interface.

Lapack

The interface is implemented in the files Seldon/computation/interfaces/Lapack_*) if you have a doubt about the syntax. The following C++ names have been chosen (in bold, name of Lapack subroutines)

  • GetLU(A, pivot) factorization of matrix A with pivot. (dgetrf, dsytrf, dsptrf, zhetrf, zhptrf)

  • SolveLU(A, pivot, X) solves linear system, assuming that GetLU has been called before. ( dgetrs , dsytrs, dsptrs, zhetrs, zhptrs, dtrtrs, dtptrs )

  • ReciprocalConditionNumber(A, pivot, norm, anorm) returns the inverse of condition number of matrix A, assuming that GetLU has been called before (dgecon, dsycon, dspcon, zhecon , zhpcon , dtrcon , dtpcon )

  • RefineSolutionLU(A, Alu, pivot, x, b, ferr, berr) improves the solution computed by SolveLU and provides error bounds and backward error estimates. (dgerfs, dsyrfs, dsprfs, zherfs, zhprfs)

  • GetInverse(A) computes the inverse of A. This method first calls GetLU in order to factorize the matrix before forming the inverse. (dgetri, dsytri, dsptri, zhetri, zhptri, dtrtri , dtptri )

  • GetScalingFactors(A, row_scale, col_scale, row_cond_number, col_cond_number, amax) computes row and column scalings intended to equilibrate a matrix and reduce its condition number. (dgeequ)

  • GetAndSolveLU(A, pivot, X) factorization of matrix A, and resolution of linear system A X = b.

  • GetQR(A, tau) QR factorization of matrix A. (dgeqrf)

  • GetLQ(A, tau) LQ factorization of matrix A. (dgelqf)

  • GetQ_FromQR(A, tau) explicits matrix Q issued from QR factorization of matrix A. It is assumed that GetQR has been called before. (dorgqr)

  • MltQ_FromQR(A, tau, x) multiplies by matrix Q issued from QR factorization of matrix A. It is assumed that GetQR has been called before. (dormqr)

  • SolveQR(A, tau, x) solves least-squares problem by using QR factorization of matrix A. It is assumed that GetQR has been called before.

  • SolveLQ(A, tau, x) solves least-squares problem by using LQ factorization of matrix A. It is assumed that GetLQ has been called before.

  • GetEigenvalues(A, lambda_real, lambda_imag) computes eigenvalues (real part and imaginary part) of unsymmetric real matrix. (dgeev )

  • GetEigenvalues(A, lambda) computes eigenvalues of matrix. (zgeev, dsyev, zheev, dspev, zhpev)

  • GetEigenvaluesEigenvectors(A, lambda_real, lambda_imag, W) computes eigenvalues (real part and imaginary part) and eigenvectors of unsymmetric real matrix.

  • GetEigenvaluesEigenvectors(A, lambda, W) computes eigenvalues and eigenvectors of matrix.

  • GetEigenvalues(A, B, lambda) solves generalized eigenvalue problem. (zggev, dsygv, zhegv, dspgv, zhpgv)

  • GetEigenvaluesEigenvectors(A, B, lambda, W) solves generalized eigenvalue problem.

  • GetSVD(A, lambda, U, V) computes singular value decomposition of matrix A. (dgesvd)

A comprehensive test of Lapack interface is done in file test/program/lapack_test.cpp.