Commit 8a084ad8 authored by Manuel Günther's avatar Manuel Günther

Removed the last occurrences of ..._ and ..._nocheck functions (now hopefully really)

parent 21c0c84a
Pipeline #10435 passed with stages
in 7 minutes and 54 seconds
......@@ -18,25 +18,25 @@
// Declaration of the external LAPACK function (Divide and conquer SVD)
extern "C" void dggsvd3_(const char *jobu,
const char *jobv,
const char *jobv,
const char *jobq,
const int *M,
const int *N,
const int *P,
const int *K,
const int *L,
double *A,
const int *lda,
double *B,
double *A,
const int *lda,
double *B,
const int *ldb,
double *alpha,
double *beta,
double *U,
const int *ldu,
const int *ldu,
double *V,
const int *ldv,
const int *ldv,
double *Q,
const int *ldq,
const int *ldq,
double *work,
const int *lwork,
int *iwork,
......@@ -61,16 +61,16 @@ void bob::math::gsvd( blitz::Array<double,2>& A,
const int M = A.extent(0);
const int N = A.extent(1);
const int P = B.extent(0);
const int lda = std::max(1,M);
const int ldb = std::max(1,P);
const int ldu = std::max(1,M);
const int ldv = std::max(1,P);
const int ldq = std::max(1,N);
int K = 0; //out
int L = 0; //out
int L = 0; //out
// Prepares to call LAPACK function:
// We will decompose A^T rather than A and B^T rather than B to reduce the required number of copy
// We recall that FORTRAN/LAPACK is column-major order whereas blitz arrays
......@@ -84,21 +84,21 @@ void bob::math::gsvd( blitz::Array<double,2>& A,
//B_lapack = B^T
blitz::Array<double,2> B_blitz_lapack(bob::core::array::ccopy(const_cast<blitz::Array<double,2>&>(B).transpose(1,0)));
double* B_lapack = B_blitz_lapack.data();
// U. We will trainpose this one in the end
double *U_lapack = U.data();
// V. We will trainpose this one in the end
double *V_lapack = V.data();
// For Q, we will return the transpose
// For Q, we will return the transpose
double *Q_lapack = Q.data();
// In LAPACK C and S is 1-d. Our code makes it diagonal
blitz::Array<double,1> C_1d(N); C_1d = 0;
double *C_lapack = C_1d.data();
blitz::Array<double,1> S_1d(N); S_1d = 0;
double *S_lapack = S_1d.data();
......@@ -109,81 +109,81 @@ void bob::math::gsvd( blitz::Array<double,2>& A,
int info = 0;
// A/ Queries the optimal size of the working array
dggsvd3_(&jobu,
&jobv,
&jobv,
&jobq,
&M,
&N,
&P,
&K,
&L,
A_lapack,
&lda,
B_lapack,
A_lapack,
&lda,
B_lapack,
&ldb,
C_lapack,
S_lapack,
U_lapack,
&ldu,
&ldu,
V_lapack,
&ldv,
&ldv,
Q_lapack,
&ldq,
&ldq,
&work_query,
&lwork_query,
iwork.get(),
&info);
if (info != 0)
throw std::runtime_error("The LAPACK dggsvd3 function returned a non-zero value during the checking");
// B/ Computes
const int lwork = static_cast<int>(work_query);
boost::shared_array<double> work(new double[lwork]);
dggsvd3_(&jobu,
&jobv,
&jobv,
&jobq,
&M,
&N,
&P,
&K,
&L,
A_lapack,
&lda,
B_lapack,
A_lapack,
&lda,
B_lapack,
&ldb,
C_lapack,
S_lapack,
U_lapack,
&ldu,
&ldu,
V_lapack,
&ldv,
&ldv,
Q_lapack,
&ldq,
&ldq,
work.get(),
&lwork,
iwork.get(),
&info);
if (info != 0)
throw std::runtime_error("The LAPACK dggsvd3 function returned a non-zero value during the computation.");
/*
/*
According to the website http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga4a187519e5c71da3b3f67c85e9baf0f2.html#ga4a187519e5c71da3b3f67c85e9baf0f2
if (M-K-L >= 0)
N-K-L K L
( 0 R ) = K ( 0 R11 R12 )
L ( 0 0 R22 )
Where R11, R12, R22, = A(1:K+L,N-K-L+1:N)
else
N-K-L K M-K K+L-M
( 0 R ) = K ( 0 R11 R12 R13 )
M-K ( 0 0 R22 R23 )
......@@ -204,11 +204,11 @@ void bob::math::gsvd( blitz::Array<double,2>& A,
X.resize(r,N);
if (M-K-L >= 0){
//1. First we need the [0 R], which has the shape (N-K-L + K + L, K+L)
//1. First we need the [0 R], which has the shape (N-K-L + K + L, K+L)
zeroR(blitz::Range(0, r-1), blitz::Range(N-r, N-1)) = A_blitz_lapack.transpose(1,0)(blitz::Range(0,r-1), blitz::Range(N-r,N-1));
// 2. Now we have to deal with C and S according to http://www.netlib.org/lapack/lug/node36.html
// In the end C is m-by-r and is p-by-r, both are real, nonnegative and diagonal and C'C + S'S = I ,
// They have the following structure
// They have the following structure
// COPY AND PASTE
//2.1 Preparing C
......@@ -216,16 +216,16 @@ void bob::math::gsvd( blitz::Array<double,2>& A,
// C = K ( I 0 )
// L ( 0 C )
// M-K-L ( 0 0 )
// A - Identity part
if (K>0){
blitz::Array<double,2> I (K,K); I= 0;
bob::math::eye_(I);
bob::math::eye(I);
C(blitz::Range(0, K-1), blitz::Range(0, K-1)) = I;
}
// B - diag(C) part. Here the C is LxL
// Swaping
bob::math::swap_(C_1d, iwork.get(), K, std::min(M,r));
blitz::Array<double,2> C_diag (L,L); C_diag = 0;
bob::math::diag(C_1d(blitz::Range(K,K+L-1)), C_diag);
......@@ -248,7 +248,7 @@ void bob::math::gsvd( blitz::Array<double,2>& A,
//1. First we need the [0 R], which has the shape (N-K-L K M-K K+L-M , k + M-K K+L-M )
//A. First part of R is in A(1:M, N-K-L+1:N)
//A. First part of R is in A(1:M, N-K-L+1:N)
zeroR(blitz::Range(0,M-1), blitz::Range(N-K-L, N-1)) = A_blitz_lapack.transpose(1,0)(blitz::Range(0,M-1) , blitz::Range(N-K-L,N-1));
//B. Second part of R is in B(M-K+1:L,N+M-K-L+1:N)
......@@ -256,7 +256,7 @@ void bob::math::gsvd( blitz::Array<double,2>& A,
//2. Now we have to deal with C and S according to http://www.netlib.org/lapack/lug/node36.html
// In the end C is m-by-r and is p-by-r, both are real, nonnegative and diagonal and C'C + S'S = I ,
// They have the following structure
// They have the following structure
// COPY AND PASTE
//2.1 Preparing C, where C=diag( ALPHA(K+1), ... , ALPHA(M) ),
......@@ -267,7 +267,7 @@ void bob::math::gsvd( blitz::Array<double,2>& A,
// A - Identity part
if (K>0){
blitz::Array<double,2> I (K,K); I = 0;
bob::math::eye_(I);
bob::math::eye(I);
C(blitz::Range(0, K-1), blitz::Range(0, K-1)) = I;
}
......@@ -290,7 +290,7 @@ void bob::math::gsvd( blitz::Array<double,2>& A,
// A - Identity part
if (K+L-M>0){
blitz::Array<double,2> I (K+L-M,K+L-M); I= 0;
bob::math::eye_(I);
bob::math::eye(I);
S(blitz::Range(M-K,L-1), blitz::Range(M,K+L-1)) = I;
}
......@@ -298,7 +298,7 @@ void bob::math::gsvd( blitz::Array<double,2>& A,
// Swaping
blitz::Array<double,1> S_1d_cropped(M-K); S_1d_cropped = 0;
S_1d_cropped = S_1d(blitz::Range(K,K+M-1));
bob::math::swap_(S_1d_cropped, iwork.get(), K, std::min(M,r));
blitz::Array<double,2> S_diag (M,M); S_diag = 0;
bob::math::diag(S_1d_cropped, S_diag);
......@@ -320,10 +320,10 @@ void bob::math::gsvd( blitz::Array<double,2>& A,
bob::math::swap_(V, iwork.get(), K, std::min(M,r));
//Computing X
bob::math::prod_(zeroR, Q, X);
bob::math::prod(zeroR, Q, X);
blitz::Array<double,2> Xt(
bob::core::array::ccopy(const_cast<blitz::Array<double,2>&>(X).transpose(1,0)));
X = Xt;
bob::math::swap_(X, iwork.get(), K, std::min(M,r));
}
}
......@@ -283,7 +283,7 @@ void bob::math::linsolveCGSympos(const blitz::Array<double,2>& A,
while (n_iter < max_iter && delta > acc*acc*delta0)
{
// q = A*d
bob::math::prod_(A, d, q);
bob::math::prod(A, d, q);
// alpha = delta/(d'*q);
double alpha = delta / bob::math::dot(d,q);
......
......@@ -46,5 +46,5 @@ void bob::math::sqrtSymReal(const blitz::Array<double,2>& A,
blitz::firstIndex i;
blitz::secondIndex j;
tmp = V(i,j) * D(j); // tmp = V.sqrt(D)
bob::math::prod_(tmp, Vt, B); // B = V.sqrt(D).V^T
bob::math::prod(tmp, Vt, B); // B = V.sqrt(D).V^T
}
......@@ -18,32 +18,11 @@
namespace bob { namespace math {
/**
* @brief Performs the matrix multiplication C=A*B
*
* @warning No checks are performed on the array sizes and is recommended
* only in scenarios where you have previously checked conformity and is
* focused only on speed.
*
* @param A The A matrix (left element of the multiplication) (size MxN)
* @param B The B matrix (right element of the multiplication) (size NxP)
* @param C The resulting matrix (size MxP)
*/
template<typename T1, typename T2, typename T3>
void prod_(const blitz::Array<T1,2>& A, const blitz::Array<T2,2>& B,
blitz::Array<T3,2>& C) {
blitz::firstIndex i;
blitz::secondIndex j;
blitz::thirdIndex k;
C = blitz::sum(A(i,k) * B(k,j), k);
}
/**
* @brief Performs the matrix multiplication C=A*B
*
* The input and output data have their sizes checked and this method will
* raise an appropriate exception if that is not cased. If you know that the
* input and output matrices conform, use the prod_() variant.
* raise an appropriate exception if that is not cased.
*
* @param A The A matrix (left element of the multiplication) (size MxN)
* @param B The B matrix (right element of the multiplication) (size NxP)
......@@ -62,34 +41,17 @@ namespace bob { namespace math {
bob::core::array::assertSameDimensionLength(A.extent(0), C.extent(0));
bob::core::array::assertSameDimensionLength(B.extent(1), C.extent(1));
prod_(A, B, C);
}
/**
* @brief Performs the matrix-vector multiplication c=A*b
*
* @warning No checks are performed on the array sizes and is recommended
* only in scenarios where you have previously checked conformity and is
* focused only on speed.
*
* @param A The A matrix (left element of the multiplication) (size MxN)
* @param b The b vector (right element of the multiplication) (size N)
* @param c The resulting vector (size M)
*/
template<typename T1, typename T2, typename T3>
void prod_(const blitz::Array<T1,2>& A, const blitz::Array<T2,1>& b,
blitz::Array<T3,1>& c) {
blitz::firstIndex i;
blitz::secondIndex j;
c = blitz::sum(A(i,j) * b(j), j);
blitz::thirdIndex k;
C = blitz::sum(A(i,k) * B(k,j), k);
}
/**
* @brief Performs the matrix-vector multiplication c=A*b
*
* The input and output data have their sizes checked and this method will
* raise an appropriate exception if that is not cased. If you know that the
* input and output matrices conform, use the prod_() variant.
* raise an appropriate exception if that is not cased.
*
* @param A The A matrix (left element of the multiplication) (size MxN)
* @param b The b vector (right element of the multiplication) (size N)
......@@ -107,34 +69,17 @@ namespace bob { namespace math {
bob::core::array::assertZeroBase(c);
bob::core::array::assertSameDimensionLength(c.extent(0), A.extent(0));
prod_(A, b, c);
}
/**
* @brief Performs the vector-matrix multiplication c=a*B
*
* @warning No checks are performed on the array sizes and is recommended
* only in scenarios where you have previously checked conformity and is
* focused only on speed.
*
* @param a The a vector (left element of the multiplication) (size M)
* @param B The B matrix (right element of the multiplication) (size MxN)
* @param c The resulting vector (size N)
*/
template<typename T1, typename T2, typename T3>
void prod_(const blitz::Array<T1,1>& a, const blitz::Array<T2,2>& B,
blitz::Array<T3,1>& c) {
blitz::firstIndex i;
blitz::secondIndex j;
c = blitz::sum(a(j) * B(j,i), j);
c = blitz::sum(A(i,j) * b(j), j);
}
/**
* @brief Performs the vector-matrix multiplication c=a*B
*
* The input and output data have their sizes checked and this method will
* raise an appropriate exception if that is not cased. If you know that the
* input and output matrices conform, use the prod_() variant.
* raise an appropriate exception if that is not cased.
*
* @param a The a vector (left element of the multiplication) (size M)
* @param B The B matrix (right element of the multiplication) (size MxN)
......@@ -152,34 +97,16 @@ namespace bob { namespace math {
bob::core::array::assertZeroBase(c);
bob::core::array::assertSameDimensionLength(c.extent(0), B.extent(1));
prod_(a, B, c);
}
/**
* @brief Performs the outer product between two vectors generating a matrix.
*
* @warning No checks are performed on the array sizes and is recommended
* only in scenarios where you have previously checked conformity and is
* focused only on speed.
*
* @param a The a vector (left element of the multiplication) (size M)
* @param b The b matrix (right element of the multiplication) (size M)
* @param C The resulting matrix (size MxM)
*/
template<typename T1, typename T2, typename T3>
void prod_(const blitz::Array<T1,1>& a, const blitz::Array<T2,1>& b,
blitz::Array<T3,2>& C) {
blitz::firstIndex i;
blitz::secondIndex j;
C = a(i) * b(j);
c = blitz::sum(a(j) * B(j,i), j);
}
/**
* @brief Performs the outer product between two vectors generating a matrix.
*
* The input and output data have their sizes checked and this method will
* raise an appropriate exception if that is not cased. If you know that the
* input and output matrices conform, use the prod_() variant.
* raise an appropriate exception if that is not cased.
*
* @param a The a vector (left element of the multiplication) (size M)
* @param b The b matrix (right element of the multiplication) (size M)
......@@ -197,23 +124,9 @@ namespace bob { namespace math {
bob::core::array::assertSameDimensionLength(C.extent(0), a.extent(0));
bob::core::array::assertSameDimensionLength(C.extent(1), b.extent(0));
prod_(a, b, C);
}
/**
* @brief Function which computes the dot product <a,b> between two 1D blitz
* array.
*
* @warning No checks are performed on the array sizes and is recommended
* only in scenarios where you have previously checked conformity and is
* focused only on speed.
*
* @param a The a vector (size N)
* @param b The b vector (size N)
*/
template<typename T1, typename T2>
T1 dot_(const blitz::Array<T1,1>& a, const blitz::Array<T2,1>& b) {
return blitz::sum(a * b);
blitz::firstIndex i;
blitz::secondIndex j;
C = a(i) * b(j);
}
/**
......@@ -234,24 +147,9 @@ namespace bob { namespace math {
bob::core::array::assertZeroBase(b);
bob::core::array::assertSameDimensionLength(a.extent(0),b.extent(0));
return dot_(a, b);
return blitz::sum(a * b);
}
/**
* @brief Computes the trace of a square matrix (the sum of all elements in
* the main diagonal).
*
* @warning No checks are performed on the array extent sizes and is
* recommended only in scenarios where you have previously checked conformity
* and is focused only on speed.
*
* @param A The input square matrix (size NxN)
*/
template<typename T> T trace_(const blitz::Array<T,2>& A) {
blitz::firstIndex i;
return blitz::sum(A(i,i));
}
/**
* @brief Computes the trace of a square matrix (the sum of all elements in
* the main diagonal).
......@@ -267,28 +165,10 @@ namespace bob { namespace math {
bob::core::array::assertZeroBase(A);
bob::core::array::assertSameDimensionLength(A.extent(0),A.extent(1));
return trace_(A);
}
/**
* @brief Computes the euclidean norm of a vector
*/
template<typename T> inline double norm(const blitz::Array<T,1>& array) {
return std::sqrt(blitz::sum(blitz::pow2(array)));
blitz::firstIndex i;
return blitz::sum(A(i,i));
}
/**
* @brief Normalizes a vector 'i' and outputs the normalized vector in 'o'.
*
* @warning This version of the normalize() method does not check for length
* consistencies and is given as an API for cases in which you have done
* already the check and is focused on speed.
*/
template<typename T1, typename T2> void normalize_
(const blitz::Array<T1,1>& i, blitz::Array<T2,1>& o) {
o = i / norm(i);
}
/**
* @brief Normalizes a vector 'i' and outputs the normalized vector in
* itself.
......@@ -311,25 +191,9 @@ namespace bob { namespace math {
blitz::Array<T2,1>& o) {
// Check input
bob::core::array::assertSameDimensionLength(i.extent(0),o.extent(0));
normalize_(i, o);
o = i / norm(i);
}
/**
* @brief Generates an eye 2D matrix. If the matrix is squared, then it
* returns the identity matrix.
*
* @warning No checks are performed on the array and is recommended
* only in scenarios where you have previously checked conformity and is
* focused only on speed.
*
* @param A The 2D destination matrix (size MxN)
*/
template<typename T>
void eye_(blitz::Array<T,2>& A) {
A = 0.;
for(int i=0; i<std::min(A.extent(0), A.extent(1)); ++i)
A(i,i) = 1.;
}
/**
* @brief Generates an eye 2D matrix. If the matrix is squared, then it
......@@ -340,24 +204,9 @@ namespace bob { namespace math {
template<typename T>
void eye(blitz::Array<T,2>& A) {
bob::core::array::assertZeroBase(A);
eye_(A);
}
/**
* @brief Generates a 2D square diagonal matrix from a 1D vector.
*
* @warning No checks are performed on the array sizes and is recommended
* only in scenarios where you have previously checked conformity and is
* focused only on speed.
*
* @param d The 1D vector which contains the diagonal (size N)
* @param A The 2D destination matrix (size NxN)
*/
template<typename T>
void diag_(const blitz::Array<T,1>& d, blitz::Array<T,2>& A) {
A = 0.;
for(int i=0; i<A.extent(0); ++i)
A(i,i) = d(i);
for(int i=0; i<std::min(A.extent(0), A.extent(1)); ++i)
A(i,i) = 1.;
}
/**
......@@ -372,23 +221,9 @@ namespace bob { namespace math {
bob::core::array::assertZeroBase(A);
bob::core::array::assertSameDimensionLength(d.extent(0),A.extent(0));
bob::core::array::assertSameDimensionLength(A.extent(0),A.extent(1));
diag_(d, A);
}
/**
* @brief Extracts the 1D diagonal from a 2D matrix.
*
* @warning No checks are performed on the array sizes and is recommended
* only in scenarios where you have previously checked conformity and is
* focused only on speed.
*
* @param A The 2D matrix matrix (size NxM)
* @param d The 1D vector which will contain the diagonal (size min(M,N))
*/
template<typename T>
void diag_(const blitz::Array<T,2>& A, blitz::Array<T,1>& d) {
blitz::firstIndex i;
d = A(i,i);
A = 0.;
for(int i=0; i<A.extent(0); ++i)
A(i,i) = d(i);
}
/**
......@@ -403,7 +238,8 @@ namespace bob { namespace math {
bob::core::array::assertZeroBase(d);
const int dim_d = std::min(A.extent(0),A.extent(1));
bob::core::array::assertSameDimensionLength(d.extent(0),dim_d);
diag_(A, d);
blitz::firstIndex i;
d = A(i,i);
}
}}
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment