Commit b326caef authored by Manuel Günther's avatar Manuel Günther

Brought back all ..._ functions in C++, but not in Python

parent da5f997b
Pipeline #10452 passed with stages
in 8 minutes and 18 seconds
......@@ -15,7 +15,11 @@
double bob::math::det(const blitz::Array<double,2>& A)
{
bob::core::array::assertSameDimensionLength(A.extent(0),A.extent(1));
return bob::math::det_(A);
}
double bob::math::det_(const blitz::Array<double,2>& A)
{
// Size variable
int N = A.extent(0);
......@@ -49,7 +53,11 @@ double bob::math::det(const blitz::Array<double,2>& A)
double bob::math::slogdet(const blitz::Array<double,2>& A, int& sign)
{
bob::core::array::assertSameDimensionLength(A.extent(0),A.extent(1));
return bob::math::slogdet_(A, sign);
}
double bob::math::slogdet_(const blitz::Array<double,2>& A, int& sign)
{
// Size variable
int N = A.extent(0);
......@@ -83,3 +91,4 @@ double bob::math::slogdet(const blitz::Array<double,2>& A, int& sign)
return Udiag;
}
......@@ -57,6 +57,16 @@ void bob::math::eig(const blitz::Array<double,2>& A,
bob::core::array::assertSameShape(V,shape2);
bob::core::array::assertSameShape(D,shape1);
bob::math::eig_(A, V, D);
}
void bob::math::eig_(const blitz::Array<double,2>& A,
blitz::Array<std::complex<double>,2>& V,
blitz::Array<std::complex<double>,1>& D)
{
// Size variable
const int N = A.extent(0);
// Prepares to call LAPACK function
// Initialises LAPACK variables
const char jobvl = 'N'; // Do NOT compute left eigen-vectors
......@@ -137,6 +147,15 @@ void bob::math::eigSym(const blitz::Array<double,2>& A,
bob::core::array::assertSameShape(V,shape2);
bob::core::array::assertSameShape(D,shape1);
bob::math::eigSym_(A, V, D);
}
void bob::math::eigSym_(const blitz::Array<double,2>& A,
blitz::Array<double,2>& V, blitz::Array<double,1>& D)
{
// Size variable
const int N = A.extent(0);
// Prepares to call LAPACK function
// Initialises LAPACK variables
const char jobz = 'V'; // Get both the eigenvalues and the eigenvectors
......@@ -215,6 +234,15 @@ void bob::math::eigSym(const blitz::Array<double,2>& A, const blitz::Array<doubl
bob::core::array::assertSameShape(V,shape2);
bob::core::array::assertSameShape(D,shape1);
bob::math::eigSym_(A, B, V, D);
}
void bob::math::eigSym_(const blitz::Array<double,2>& A, const blitz::Array<double,2>& B,
blitz::Array<double,2>& V, blitz::Array<double,1>& D)
{
// Size variable
const int N = A.extent(0);
// Prepares to call LAPACK function
// Initialises LAPACK variables
const int itype = 1;
......
......@@ -35,6 +35,14 @@ void bob::math::inv(const blitz::Array<double,2>& A, blitz::Array<double,2>& B)
bob::core::array::assertSameShape(A,shapeA);
bob::core::array::assertSameShape(B,shapeA);
bob::math::inv_(A, B);
}
void bob::math::inv_(const blitz::Array<double,2>& A, blitz::Array<double,2>& B)
{
// Size variable
const int N = A.extent(0);
//////////////////////////////////////
// Prepares to call LAPACK functions
// Initializes LAPACK variables
......@@ -88,3 +96,4 @@ void bob::math::inv(const blitz::Array<double,2>& A, blitz::Array<double,2>& B)
if (!B_direct_use)
B = A_blitz_lapack;
}
......@@ -35,6 +35,12 @@ void bob::math::linsolve(const blitz::Array<double,2>& A,
bob::core::array::assertSameDimensionLength(A.extent(0), A.extent(1));
bob::core::array::assertSameDimensionLength(A.extent(1), b.extent(0));
bob::math::linsolve_(A, b, x);
}
void bob::math::linsolve_(const blitz::Array<double,2>& A,
const blitz::Array<double,1>& b, blitz::Array<double,1>& x)
{
// Defines dimensionality variables
const int N = A.extent(0);
......@@ -88,6 +94,13 @@ void bob::math::linsolve(const blitz::Array<double,2>& A,
bob::core::array::assertSameDimensionLength(A.extent(0), B.extent(0));
bob::core::array::assertSameDimensionLength(X.extent(1), B.extent(1));
bob::math::linsolve_(A, B, X);
}
void bob::math::linsolve_(const blitz::Array<double,2>& A,
const blitz::Array<double,2>& B, blitz::Array<double,2>& X)
{
// Defines dimensionality variables
const int N = A.extent(0);
const int P = X.extent(1);
......@@ -148,6 +161,13 @@ void bob::math::linsolveSympos(const blitz::Array<double,2>& A,
bob::core::array::assertSameDimensionLength(A.extent(0), A.extent(1));
bob::core::array::assertSameDimensionLength(A.extent(1), b.extent(0));
bob::math::linsolveSympos_(A, b, x);
}
void bob::math::linsolveSympos_(const blitz::Array<double,2>& A,
const blitz::Array<double,1>& b, blitz::Array<double,1>& x)
{
// Defines dimensionality variables
const int N = A.extent(0);
......@@ -202,6 +222,12 @@ void bob::math::linsolveSympos(const blitz::Array<double,2>& A,
bob::core::array::assertSameDimensionLength(A.extent(0), B.extent(0));
bob::core::array::assertSameDimensionLength(X.extent(1), B.extent(1));
bob::math::linsolveSympos_(A, B, X);
}
void bob::math::linsolveSympos_(const blitz::Array<double,2>& A,
const blitz::Array<double,2>& B, blitz::Array<double,2>& X)
{
// Defines dimensionality variables
const int N = A.extent(0);
const int P = X.extent(1);
......@@ -268,6 +294,16 @@ void bob::math::linsolveCGSympos(const blitz::Array<double,2>& A,
bob::core::array::assertSameDimensionLength(A.extent(0), N);
bob::core::array::assertSameDimensionLength(A.extent(1), N);
bob::math::linsolveCGSympos_(A, b, x, acc, max_iter);
}
void bob::math::linsolveCGSympos_(const blitz::Array<double,2>& A,
const blitz::Array<double,1>& b, blitz::Array<double,1>& x,
const double acc, const int max_iter)
{
// Dimensionality of the problem
const int N = b.extent(0);
blitz::Array<double,1> r(N), d(N), best_x(N), q(N), tmp(N);
x = 0.;
r = b;
......
......@@ -45,6 +45,17 @@ void bob::math::lu(const blitz::Array<double,2>& A, blitz::Array<double,2>& L,
bob::core::array::assertSameShape(U,shapeU);
bob::core::array::assertSameShape(P,shapeP);
bob::math::lu_(A, L, U, P);
}
void bob::math::lu_(const blitz::Array<double,2>& A, blitz::Array<double,2>& L,
blitz::Array<double,2>& U, blitz::Array<double,2>& P)
{
// Size variable
const int M = A.extent(0);
const int N = A.extent(1);
const int minMN = std::min(M,N);
// Prepares to call LAPACK function
// Initialises LAPACK variables
......@@ -105,6 +116,15 @@ void bob::math::chol(const blitz::Array<double,2>& A,
bob::core::array::assertSameDimensionLength(M,N);
bob::core::array::assertSameShape(A,L);
bob::math::chol_(A, L);
}
void bob::math::chol_(const blitz::Array<double,2>& A,
blitz::Array<double,2>& L)
{
// Size variable
const int N = A.extent(0);
// Prepares to call LAPACK function
// Initialises LAPACK variables
int info = 0;
......@@ -141,3 +161,4 @@ void bob::math::chol(const blitz::Array<double,2>& A,
blitz::secondIndex j;
L = blitz::where(i < j, 0, L);
}
......@@ -11,6 +11,13 @@
#include <bob.core/assert.h>
void bob::math::pavx(const blitz::Array<double,1>& y, blitz::Array<double,1>& ghat)
{
bob::core::array::assertSameShape(y, ghat);
assert(y.extent(0) > 0);
math::pavx_(y, ghat);
}
static size_t pavx_1(const blitz::Array<double,1>& y, blitz::Array<double,1>& ghat,
blitz::Array<size_t,1>& index, blitz::Array<size_t,1>& len)
{
......@@ -57,11 +64,8 @@ static void pavx_2(blitz::Array<double,1>& ghat, blitz::Array<size_t,1>& index,
}
}
void bob::math::pavx(const blitz::Array<double,1>& y, blitz::Array<double,1>& ghat)
void bob::math::pavx_(const blitz::Array<double,1>& y, blitz::Array<double,1>& ghat)
{
bob::core::array::assertSameShape(y, ghat);
assert(y.extent(0) > 0);
// Define working arrays: An interval of indices is represented by its left
// endpoint "index" and its length "len"
int N = y.extent(0);
......
......@@ -18,7 +18,6 @@ void bob::math::pinv(const blitz::Array<double,2>& A,
// Size variables
const int M = A.extent(0);
const int N = A.extent(1);
const int nb_singular = std::min(M,N);
// Checks zero base
bob::core::array::assertZeroBase(A);
......@@ -27,6 +26,18 @@ void bob::math::pinv(const blitz::Array<double,2>& A,
bob::core::array::assertSameDimensionLength(B.extent(0), N);
bob::core::array::assertSameDimensionLength(B.extent(1), M);
// Calls pinv_()
bob::math::pinv_(A, B, rcond);
}
void bob::math::pinv_(const blitz::Array<double,2>& A,
blitz::Array<double,2>& B, const double rcond)
{
// Size variables
const int M = A.extent(0);
const int N = A.extent(1);
const int nb_singular = std::min(M,N);
// Allocates arrays for the SVD
blitz::Array<double,2> U(M,M);
blitz::Array<double,1> sigma(nb_singular);
......@@ -35,7 +46,7 @@ void bob::math::pinv(const blitz::Array<double,2>& A,
blitz::Array<double,2> sigmai_Ut(N, M);
blitz::Array<double,1> Ut_sum(M);
// Computes the SVD
bob::math::svd(A, U, sigma, Vt);
bob::math::svd_(A, U, sigma, Vt);
// Cuts off small sigmas for the inverse similar to Numpy:
// cf. https://github.com/numpy/numpy/blob/maintenance/1.7.x/numpy/linalg/linalg.py#L1577
......@@ -57,3 +68,4 @@ void bob::math::pinv(const blitz::Array<double,2>& A,
blitz::Array<double,2> V = Vt.transpose(1,0);
bob::math::prod(V, sigmai_Ut, B);
}
......@@ -28,13 +28,22 @@ void bob::math::sqrtSymReal(const blitz::Array<double,2>& A,
bob::core::array::assertSameShape(A,shape);
bob::core::array::assertSameShape(B,shape);
bob::math::sqrtSymReal_(A, B);
}
void bob::math::sqrtSymReal_(const blitz::Array<double,2>& A,
blitz::Array<double,2>& B)
{
// Size variable
int N = A.extent(0);
// 1/ Perform the Eigenvalue decomposition of the symmetric matrix
// A = V.D.V^T, and V^-1=V^T
blitz::Array<double,2> V(N,N);
blitz::Array<double,2> Vt = V.transpose(1,0);
blitz::Array<double,1> D(N);
blitz::Array<double,2> tmp(N,N); // Cache for multiplication
bob::math::eigSym(A,V,D);
bob::math::eigSym_(A,V,D);
// 2/ Updates the diagonal matrix D, such that D=sqrt(|D|)
// |.| is used to deal with values close to zero (-epsilon)
......
......@@ -79,10 +79,10 @@ static void svd_lapack( const char jobz, const int M, const int N,
if (info != 0)
throw std::runtime_error("The LAPACK dgesdd function returned a non-zero value. You may consider using LAPACK dgsevd instead (see #171) by enabling the 'safe' option.");
}
// Defining the sign of the eigenvectors
// Approch extracted from page 8 - http://prod.sandia.gov/techlib/access-control.cgi/2007/076422.pdf
if(U[0] < 0){
// Approch extracted from page 8 - http://prod.sandia.gov/techlib/access-control.cgi/2007/076422.pdf
if(U[0] < 0){
int ucol=0; ucol= (jobz=='A')? M : std::min(M,N);
for (int i=0; i<ldu*ucol; i++){
U[i] = -1*U[i];
......@@ -114,6 +114,17 @@ void bob::math::svd(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
bob::core::array::assertSameDimensionLength(Vt.extent(0), N);
bob::core::array::assertSameDimensionLength(Vt.extent(1), N);
bob::math::svd_(A, U, sigma, Vt, safe);
}
void bob::math::svd_(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
blitz::Array<double,1>& sigma, blitz::Array<double,2>& Vt, bool safe)
{
// Size variables
const int M = A.extent(0);
const int N = A.extent(1);
const int nb_singular = std::min(M,N);
// Prepares to call LAPACK function:
// We will decompose A^T rather than A to reduce the required number of copy
// We recall that FORTRAN/LAPACK is column-major order whereas blitz arrays
......@@ -150,7 +161,7 @@ void bob::math::svd(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
// Call the LAPACK function
svd_lapack(jobz, N, M, A_lapack, lda, S_lapack, U_lapack, ldu,
VT_lapack, ldvt, safe);
VT_lapack, ldvt, safe);
// Copy singular vectors back to U, V and sigma if required
......@@ -177,6 +188,17 @@ void bob::math::svd(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
bob::core::array::assertSameDimensionLength(U.extent(1), nb_singular);
bob::core::array::assertSameDimensionLength(sigma.extent(0), nb_singular);
bob::math::svd_(A, U, sigma, safe);
}
void bob::math::svd_(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
blitz::Array<double,1>& sigma, bool safe)
{
// Size variables
const int M = A.extent(0);
const int N = A.extent(1);
const int nb_singular = std::min(M,N);
// Prepares to call LAPACK function
// Initialises LAPACK variables
......@@ -227,6 +249,16 @@ void bob::math::svd(const blitz::Array<double,2>& A, blitz::Array<double,1>& sig
// Checks and resizes if required
bob::core::array::assertSameDimensionLength(sigma.extent(0), nb_singular);
bob::math::svd_(A, sigma, safe);
}
void bob::math::svd_(const blitz::Array<double,2>& A, blitz::Array<double,1>& sigma, bool safe)
{
// Size variables
const int M = A.extent(0);
const int N = A.extent(1);
const int nb_singular = std::min(M,N);
// Prepares to call LAPACK function
// Initialises LAPACK variables
......
......@@ -20,6 +20,12 @@ namespace bob { namespace math {
* @param A The A matrix to consider (size NxN)
*/
double det(const blitz::Array<double,2>& A);
/**
* @brief Function which computes the determinant of a square matrix
* @param A The A matrix to consider (size NxN)
* @warning Does not check the input matrix
*/
double det_(const blitz::Array<double,2>& A);
/**
* @brief Function which computes the sign and (natural) logarithm of
......@@ -29,6 +35,15 @@ namespace bob { namespace math {
* (-1 if negative, 0 if zero, +1 if positive)
*/
double slogdet(const blitz::Array<double,2>& A, int& sign);
/**
* @brief Function which computes the sign and (natural) logarithm of
* the determinant of a square matrix
* @param A The A matrix to consider (size NxN)
* @param sign The (output) sign of the determinant
* (-1 if negative, 0 if zero, +1 if positive)
* @warning Does not check the input matrix
*/
double slogdet_(const blitz::Array<double,2>& A, int& sign);
}}
......
......@@ -32,6 +32,23 @@ namespace bob { namespace math {
blitz::Array<std::complex<double>,2>& V,
blitz::Array<std::complex<double>,1>& D);
/**
* @brief Function which performs an eigenvalue decomposition of a real
* matrix, using the LAPACK function <code>dgeev</code>. This version does
* <b>NOT</b> perform any checks on the input data and should be, therefore,
* faster.
*
* @param A The A matrix to decompose (size NxN)
*
* @param V The V matrix of eigenvectors (size NxN) stored in columns
*
* @param D The vector of eigenvalues (size N) (in ascending order). Note
* that this is a 1D array rather than a 2D diagonal matrix!
*/
void eig_(const blitz::Array<double,2>& A,
blitz::Array<std::complex<double>,2>& V,
blitz::Array<std::complex<double>,1>& D);
/**
* @brief Function which performs an eigenvalue decomposition of a real
* symmetric matrix, using the LAPACK function <code>dsyevd</code>.
......@@ -48,6 +65,24 @@ namespace bob { namespace math {
void eigSym(const blitz::Array<double,2>& A, blitz::Array<double,2>& V,
blitz::Array<double,1>& D);
/**
* @brief Function which performs an eigenvalue decomposition of a real
* symmetric matrix, using the LAPACK function <code>dsyevd</code>. This
* version does <b>NOT</b> perform any checks on the input data and should
* be, therefore, faster.
*
* @warning The input matrix should be symmetric.
*
* @param A The A matrix to decompose (size NxN)
*
* @param V The V matrix of eigenvectors (size NxN) stored in columns
*
* @param D The vector of eigenvalues (size N) (in ascending order). Note
* that this is a 1D array rather than a 2D diagonal matrix!
*/
void eigSym_(const blitz::Array<double,2>& A, blitz::Array<double,2>& V,
blitz::Array<double,1>& D);
/**
* @brief Computes all the eigenvalues and the eigenvectors of a real
* generalized symmetric-definite eigenproblem, of the form:
......@@ -68,6 +103,28 @@ namespace bob { namespace math {
void eigSym(const blitz::Array<double,2>& A, const blitz::Array<double,2>& B,
blitz::Array<double,2>& V, blitz::Array<double,1>& D);
/**
* @brief Computes all the eigenvalues and the eigenvectors of a real
* generalized symmetric-definite eigenproblem, of the form:
* A*x=(lambda)*B*x, using the LAPACK function <code>dsygvd</code>. This
* version does <b>NOT</b> perform any checks on the input data and should
* be, therefore, faster.
*
* @warning The input matrices A and B are assumed to be symmetric and B is
* also positive definite.
*
* @param A The A input matrix (size NxN) of the problem
*
* @param B The B input matrix (size NxN) of the problem
*
* @param V The V matrix of eigenvectors (size NxN) stored in columns
*
* @param D The vector of eigenvalues (size N) (in ascending order). Note
* that this is a 1D array rather than a 2D diagonal matrix!
*/
void eigSym_(const blitz::Array<double,2>& A, const blitz::Array<double,2>& B,
blitz::Array<double,2>& V, blitz::Array<double,1>& D);
}}
#endif /* BOB_MATH_EIG_H */
......@@ -25,14 +25,12 @@ namespace bob { namespace math {
* @param input The input blitz array
* @param g The output blitz array for the gradient
* @param dx The sample distance along the x-axis
* @warning Does not check that g has the same size as input
*/
template <typename T, typename U>
void gradient(const blitz::Array<T,1>& input, blitz::Array<U,1>& g,
void gradient_(const blitz::Array<T,1>& input, blitz::Array<U,1>& g,
const double dx=1.)
{
// Check input size
bob::core::array::assertSameShape(input, g);
const int M=input.extent(0);
// Check input
if (M<2) {
......@@ -66,6 +64,24 @@ namespace bob { namespace math {
if (dx!=1.) g *= (1./dx);
}
/**
* @brief Function which computes the gradient of a 1D signal
* The gradient is computed using central differences in the interior
* and first differences at the boundaries.
* Similar to NumPy and MATLAB gradient function
* @param input The input blitz array
* @param dx The sample distance along the x-axis
* @param g The output blitz array for the gradient
*/
template <typename T, typename U>
void gradient(const blitz::Array<T,1>& input, blitz::Array<U,1>& g,
const double dx=1.)
{
// Check input size
bob::core::array::assertSameShape(input, g);
gradient_<T,U>(input, g, dx);
}
/**
* @brief Function which computes the gradient of a 2D signal
* The gradient is computed using central differences in the interior
......@@ -76,15 +92,13 @@ namespace bob { namespace math {
* @param gx The output blitz array for the gradient along the x-axis
* @param dy The sample distance along the y-axis
* @param dx The sample distance along the x-axis
* @warning Does not check that gx and gy have the same size as input
* @warning For unsigned input datatypes, this function might not work as expected.
*/
template <typename T, typename U>
void gradient(const blitz::Array<T,2>& input, blitz::Array<U,2>& gy,
void gradient_(const blitz::Array<T,2>& input, blitz::Array<U,2>& gy,
blitz::Array<U,2>& gx, const double dy=1., const double dx=1.)
{
bob::core::array::assertSameShape(input, gy);
bob::core::array::assertSameShape(input, gx);
const int M=input.extent(0);
const int N=input.extent(1);
// Check input
......@@ -143,6 +157,27 @@ namespace bob { namespace math {
if (dy!=1.) gy *= (1./dy);
if (dx!=1.) gx *= (1./dx);
}
/**
* @brief Function which computes the gradient of a 2D signal
* The gradient is computed using central differences in the interior
* and first differences at the boundaries.
* Similar to NumPy and MATLAB gradient function
* @param input The input blitz array
* @param gy The output blitz array for the gradient along the y-axis
* @param gx The output blitz array for the gradient along the x-axis
* @param dy The sample distance along the y-axis
* @param dx The sample distance along the x-axis
* @warning For unsigned input datatypes, this function might not work as expected.
*/
template <typename T, typename U>
void gradient(const blitz::Array<T,2>& input, blitz::Array<U,2>& gy,
blitz::Array<U,2>& gx, const double dy=1., const double dx=1.)
{
// Check input size
bob::core::array::assertSameShape(input, gy);
bob::core::array::assertSameShape(input, gx);
gradient_<T,U>(input, gy, gx, dy, dx);
}
/**
* @brief Function which computes the gradient of a 3D signal
......@@ -156,17 +191,14 @@ namespace bob { namespace math {
* @param dz The sample distance along the z-axis
* @param dy The sample distance along the y-axis
* @param dx The sample distance along the x-axis
* @warning Does not check that gx and gy have the same size as input
* @warning For unsigned input datatypes, this function might not work as expected.
*/
template <typename T, typename U>
void gradient(const blitz::Array<T,3>& input, blitz::Array<U,3>& gz,
void gradient_(const blitz::Array<T,3>& input, blitz::Array<U,3>& gz,
blitz::Array<U,3>& gy, blitz::Array<U,3>& gx, const double dz=1.,
const double dy=1., const double dx=1.)
{
bob::core::array::assertSameShape(input, gz);
bob::core::array::assertSameShape(input, gy);
bob::core::array::assertSameShape(input, gx);
const int M=input.extent(0);
const int N=input.extent(1);
const int P=input.extent(2);
......@@ -248,6 +280,32 @@ namespace bob { namespace math {
if (dy!=1.) gy *= (1./dy);
if (dx!=1.) gx *= (1./dx);
}
/**
* @brief Function which computes the gradient of a 3D signal
* The gradient is computed using central differences in the interior
* and first differences at the boundaries.
* Similar to NumPy and MATLAB gradient function
* @param input The input blitz array
* @param gz The output blitz array for the gradient along the z-axis
* @param gy The output blitz array for the gradient along the y-axis
* @param gx The output blitz array for the gradient along the x-axis
* @param dz The sample distance along the z-axis
* @param dy The sample distance along the y-axis
* @param dx The sample distance along the x-axis
* @warning For unsigned input datatypes, this function might not work as expected.
*/
template <typename T, typename U>
void gradient(const blitz::Array<T,3>& input, blitz::Array<U,3>& gz,
blitz::Array<U,3>& gy, blitz::Array<U,3>& gx, const double dz=1.,
const double dy=1., const double dx=1.)
{
// Check input size
bob::core::array::assertSameShape(input, gz);
bob::core::array::assertSameShape(input, gy);
bob::core::array::assertSameShape(input, gx);
gradient_<T,U>(input, gz, gy, gx, dz, dy, dx);
}
}}
#endif /* BOB_MATH_GRADIENT_H */
......@@ -21,6 +21,7 @@ namespace bob { namespace math {
* @param B The B=inverse(A) matrix (size NxN)
*/
void inv(const blitz::Array<double,2>& A, blitz::Array<double,2>& B);
void inv_(const blitz::Array<double,2>& A, blitz::Array<double,2>& B);
}}
......
......@@ -25,6 +25,8 @@ namespace bob { namespace math {
*/
void linsolve(const blitz::Array<double,2>& A,
const blitz::Array<double,1>& b, blitz::Array<double,1>& x);
void linsolve_(const blitz::Array<double,2>& A,
const blitz::Array<double,1>& b, blitz::Array<double,1>& x);
/**
* @brief Function which solves a linear system of equation using the
......@@ -36,6 +38,8 @@ namespace bob { namespace math {
*/
void linsolve(const blitz::Array<double,2>& A,
const blitz::Array<double,2>& B, blitz::Array<double,2>& X);
void linsolve_(const blitz::Array<double,2>& A,
const blitz::Array<double,2>& B, blitz::Array<double,2>& X);
/**
* @brief Function which solves a symmetric positive definite linear
......@@ -50,6 +54,8 @@ namespace bob { namespace math {
*/
void linsolveSympos(const blitz::Array<double,2>& A,
const blitz::Array<double,1>& b, blitz::Array<double,1>& x);
void linsolveSympos_(const blitz::Array<double,2>& A,
const blitz::Array<double,1>& b, blitz::Array<double,1>& x);
/**
* @brief Function which solves a symmetric positive definite linear
......@@ -64,6 +70,8 @@ namespace bob { namespace math {
*/
void linsolveSympos(const blitz::Array<double,2>& A,
const blitz::Array<double,2>& B, blitz::Array<double,2>& X);
void linsolveSympos_(const blitz::Array<double,2>& A,
const blitz::Array<double,2>& B, blitz::Array<double,2>& X);
/**
......@@ -81,6 +89,9 @@ namespace bob { namespace math {
void linsolveCGSympos(const blitz::Array<double,2>& A,
const blitz::Array<double,1>& b, blitz::Array<double,1>& x,
const double acc, const int max_iter);
void linsolveCGSympos_(const blitz::Array<double,2>& A,
const blitz::Array<double,1>& b, blitz::Array<double,1>& x,
const double acc, const int max_iter);
}}
......
......@@ -22,6 +22,8 @@ namespace bob { namespace math {
*/
void lu(const blitz::Array<double,2>& A, blitz::Array<double,2>& L,
blitz::Array<double,2>& U, blitz::Array<double,2>& P);
void lu_(const blitz::Array<double,2>& A, blitz::Array<double,2>& L,
blitz::Array<double,2>& U, blitz::Array<double,2>& P);
/**
* @brief Performs the Cholesky decomposition of a real symmetric
......@@ -33,6 +35,7 @@ namespace bob { namespace math {
* @param L The L lower-triangular matrix of the decomposition
*/
void chol(const blitz::Array<double,2>& A, blitz::Array<double,2>& L);
void chol_(const blitz::Array<double,2>& A, blitz::Array<double,2>& L);
/**
* @}
*/
......
......@@ -33,6 +33,7 @@ namespace bob { namespace math {
* Duembgen is that the data points of the vector y can not be weighted.
*/
void pavx(const blitz::Array<double,1>& y, blitz::Array<double,1>& ghat);
void pavx_(const blitz::Array<double,1>& y, blitz::Array<double,1>& ghat);