diff --git a/bob/math/cpp/gsvd.cpp b/bob/math/cpp/gsvd.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..47e0a5fc0dfa84ec19d6f8fcd3b47b75c8b00c64
--- /dev/null
+++ b/bob/math/cpp/gsvd.cpp
@@ -0,0 +1,328 @@
+/**
+ * @date Tue Jan 10 21:14 2016 +0100
+ * @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#include <stdexcept>
+#include <boost/shared_array.hpp>
+
+#include <bob.math/gsvd.h>
+
+#include <bob.core/assert.h>
+#include <bob.core/check.h>
+#include <bob.core/array_copy.h>
+#include <bob.math/linear.h>
+
+
+// Declaration of the external LAPACK function (Divide and conquer SVD)
+extern "C" void dggsvd3_(const char *jobu,
+ 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,
+ const int *ldb,
+ double *alpha,
+ double *beta,
+ double *U,
+ const int *ldu,
+ double *V,
+ const int *ldv,
+ double *Q,
+ const int *ldq,
+ double *work,
+ const int *lwork,
+ int *iwork,
+ int *info);
+
+
+void bob::math::gsvd( blitz::Array<double,2>& A,
+ blitz::Array<double,2>& B,
+ blitz::Array<double,2>& U,
+ blitz::Array<double,2>& V,
+ blitz::Array<double,2>& zeroR,
+ blitz::Array<double,2>& Q,
+ blitz::Array<double,2>& X,
+ blitz::Array<double,2>& C,
+ blitz::Array<double,2>& S)
+{
+ const char jobu = 'U';
+ const char jobv = 'V';
+ const char jobq = 'Q';
+
+ // Size variables
+ 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
+
+ // 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
+ // are row-major order by default.
+
+ //A_lapack = A^T
+ blitz::Array<double,2> A_blitz_lapack(bob::core::array::ccopy(const_cast<blitz::Array<double,2>&>(A).transpose(1,0)));
+ double* A_lapack = A_blitz_lapack.data();
+
+
+ //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
+ 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();
+
+ const int lwork_query = -1;
+ double work_query;
+
+ boost::shared_array<int> iwork(new int[N]);
+
+ int info = 0;
+ // A/ Queries the optimal size of the working array
+
+ dggsvd3_(&jobu,
+ &jobv,
+ &jobq,
+ &M,
+ &N,
+ &P,
+ &K,
+ &L,
+ A_lapack,
+ &lda,
+ B_lapack,
+ &ldb,
+ C_lapack,
+ S_lapack,
+ U_lapack,
+ &ldu,
+ V_lapack,
+ &ldv,
+ Q_lapack,
+ &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,
+ &jobq,
+ &M,
+ &N,
+ &P,
+ &K,
+ &L,
+ A_lapack,
+ &lda,
+ B_lapack,
+ &ldb,
+ C_lapack,
+ S_lapack,
+ U_lapack,
+ &ldu,
+ V_lapack,
+ &ldv,
+ Q_lapack,
+ &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 )
+ K+L-M ( 0 0 0 R33 )
+
+
+ where R11, R12, R13, R22, R23 = A(1:M, N-K-L+1:N)
+ where R33 = B(M-K+1:L,N+M-K-L+1:N)
+
+ */
+ int r = K + L; //Dimension of R
+ C.resize(M, r);
+ S.resize(P, r);
+ S = 0;
+ C = 0;
+ zeroR.resize(r, N);
+ zeroR = 0;
+ 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)
+ 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
+ // COPY AND PASTE
+
+ //2.1 Preparing C
+ // K L
+ // 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);
+ 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(0,L-1)), C_diag);
+ C(blitz::Range(K, M-1), blitz::Range(K, K+L-1)) = C_diag;
+
+ //2.2 Preparing S
+ // K L
+ //D2 = L ( 0 S )
+ // P-L ( 0 0 )
+
+ // A - diag(S) part
+ // Swap
+ bob::math::swap_(S_1d, iwork.get(), K, std::min(M,r));
+ blitz::Array<double,2> S_diag (L,L); S_diag = 0;
+ bob::math::diag(S_1d(blitz::Range(0,L-1)), S_diag);
+ S(blitz::Range(0, L-1), blitz::Range(K, K+L-1)) = S_diag;
+
+ }
+ else{
+
+ //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)
+ 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)
+ zeroR(blitz::Range(M, r - 1), blitz::Range(N-L+M-K,N-1)) = B_blitz_lapack.transpose(1,0)(blitz::Range(M-K, L-1) , blitz::Range(N+M-K-L, 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
+ // COPY AND PASTE
+
+ //2.1 Preparing C, where C=diag( ALPHA(K+1), ... , ALPHA(M) ),
+ // K M-K K+L-M
+ // D1 = K ( I 0 0 )
+ // M-K ( 0 C 0 )
+
+ // A - Identity part
+ if (K>0){
+ blitz::Array<double,2> I (K,K); I = 0;
+ bob::math::eye_(I);
+ C(blitz::Range(0, K-1), blitz::Range(0, K-1)) = I;
+ }
+
+ // B - diag(C) part
+ // Swaping
+ blitz::Array<double,1> C_1d_cropped(M-K); C_1d_cropped = 0;
+ C_1d_cropped = C_1d(blitz::Range(K,K+M-1));
+ bob::math::swap_(C_1d_cropped, iwork.get(), K, std::min(M,r));
+ blitz::Array<double,2> C_diag (M,M); C_diag = 0;
+ bob::math::diag(C_1d_cropped, C_diag);
+ C(blitz::Range(K,M-1), blitz::Range(K,M-1)) = C_diag;
+
+
+ //2.2 Preparing S
+ // K M-K K+L-M
+ // D2 = M-K ( 0 S 0 )
+ // K+L-M ( 0 0 I )
+ // P-L ( 0 0 0 )
+
+ // 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);
+ S(blitz::Range(M-K,L-1), blitz::Range(M,K+L-1)) = I;
+ }
+
+ // B - diag(S) part
+ // 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);
+ S(blitz::Range(0,M-K-1), blitz::Range(K,M-1)) = S_diag;
+ }
+
+ // Transposing U and V
+ blitz::Array<double,2> Ut(
+ bob::core::array::ccopy(const_cast<blitz::Array<double,2>&>(U).transpose(1,0)));
+ U = Ut;
+
+ blitz::Array<double,2> Vt(
+ bob::core::array::ccopy(const_cast<blitz::Array<double,2>&>(V).transpose(1,0)));
+ V = Vt;
+
+ // Swaping U
+ bob::math::swap_(U, iwork.get(), K, std::min(M,r));
+ // Swaping V
+ bob::math::swap_(V, iwork.get(), K, std::min(M,r));
+
+ //Computing 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));
+
+}
diff --git a/bob/math/gsvd.cpp b/bob/math/gsvd.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..aa20a8156b06b67c1e39b8fa278d5230640ab81d
--- /dev/null
+++ b/bob/math/gsvd.cpp
@@ -0,0 +1,77 @@
+/**
+ * @date Tue Jan 10 21:14 2016 +0100
+ * @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
+ *
+ * @brief Binds the Generalized SVD
+ */
+
+#include "gsvd.h"
+#include <bob.blitz/cppapi.h>
+#include <bob.blitz/cleanup.h>
+#include <bob.math/gsvd.h>
+
+PyObject* py_gsvd (PyObject*, PyObject* args, PyObject* kwds) {
+
+ /* Parses input arguments in a single shot */
+ static const char* const_kwlist[] = { "A", "B", 0 /* Sentinel */ };
+ static char** kwlist = const_cast<char**>(const_kwlist);
+
+ PyBlitzArrayObject* A = 0;
+ PyBlitzArrayObject* B = 0;
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&O&", kwlist,
+ &PyBlitzArray_Converter, &A,
+ &PyBlitzArray_Converter, &B
+ ))
+ return 0;
+
+ auto A_ = make_safe(A);
+ auto B_ = make_safe(B);
+
+ if (A->ndim != 2 || A->type_num != NPY_FLOAT64) {
+ PyErr_Format(PyExc_TypeError, "`A` matrix only supports 2D 64-bit float array");
+ return 0;
+ }
+
+ if (B->ndim != 2 || B->type_num != NPY_FLOAT64) {
+ PyErr_Format(PyExc_TypeError, "`B` matrix only supports 2D 64-bit float array");
+ return 0;
+ }
+
+ auto A_bz = PyBlitzArrayCxx_AsBlitz<double,2>(A);
+ auto B_bz = PyBlitzArrayCxx_AsBlitz<double,2>(B);
+
+ const int M = A_bz->extent(0);
+ const int N = A_bz->extent(1);
+ const int P = B_bz->extent(0);
+
+ // Creating the output matrices
+ blitz::Array<double,2> U(M, M); U=0;
+ blitz::Array<double,2> V(P, P); V=0;
+ blitz::Array<double,2> Q(N, N); Q=0;
+ blitz::Array<double,2> zeroR(N, N); zeroR=0;
+ blitz::Array<double,2> X(N, N); X=0;
+ blitz::Array<double,2> C(N, N); C=0;
+ blitz::Array<double,2> S(N, N); S=0;
+
+
+ try {
+ bob::math::gsvd(*A_bz,*B_bz,U,V,zeroR,Q,X,C,S);
+ return Py_BuildValue("NNNNN",
+ PyBlitzArrayCxx_AsConstNumpy(U),
+ PyBlitzArrayCxx_AsConstNumpy(V),
+ PyBlitzArrayCxx_AsConstNumpy(X),
+ PyBlitzArrayCxx_AsConstNumpy(C),
+ PyBlitzArrayCxx_AsConstNumpy(S));
+ }
+ catch (std::exception& e) {
+ PyErr_SetString(PyExc_RuntimeError, e.what());
+ }
+ catch (...) {
+ PyErr_SetString(PyExc_RuntimeError, "norminv failed: unknown exception caught");
+ }
+
+
+ return 0;
+
+}
diff --git a/bob/math/gsvd.h b/bob/math/gsvd.h
new file mode 100644
index 0000000000000000000000000000000000000000..cc38f9fc38360e041536a83cca82c2e6ad3f33e8
--- /dev/null
+++ b/bob/math/gsvd.h
@@ -0,0 +1,10 @@
+/**
+ * @author Andre Anjos <andre.anjos@idiap.ch>
+ * @date Thu 5 Dec 12:01:57 2013
+ *
+ * @brief Declaration of components relevant for main.cpp
+ */
+
+#include <Python.h>
+
+PyObject* py_gsvd(PyObject*, PyObject* args, PyObject* kwds);
diff --git a/bob/math/include/bob.math/gsvd.h b/bob/math/include/bob.math/gsvd.h
new file mode 100644
index 0000000000000000000000000000000000000000..4da252418ec7ffde991138cc8639dafb6e0be3fb
--- /dev/null
+++ b/bob/math/include/bob.math/gsvd.h
@@ -0,0 +1,107 @@
+/**
+ * @date Tue Jan 10 21:14 2016 +0100
+ * @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
+ *
+ * @brief This file defines a function to determine the SVD decomposition
+ * a 2D blitz array using LAPACK.
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#ifndef BOB_MATH_GSVD_H
+#define BOB_MATH_GSVD_H
+
+#include <blitz/array.h>
+
+namespace bob { namespace math {
+/**
+ * @ingroup MATH
+ * @{
+ */
+
+/**
+ * @brief Function which performs the Generalized Singular Value Decomposition
+ * using the 'simple' driver routine dggsvd3 of LAPACK.
+ * Just remembering that given A and B as input we have:
+ *
+ * A = U C X and B = V S X
+ * where X = [0,R]Q^T
+ *
+ *
+ *
+ * @warning The output blitz::array sigma should have the correct
+ * size, with zero base index. Checks are performed.
+ * @param A The A matrix to decompose (size MxP)
+ * @param B The B matrix to decompose (size NxP)
+ * @warning A and B must have the same number of columns, but may have different numbers of rows. If A is m-by-p and B is n-by-p, then U is m-by-m, V is n-by-n, X is p-by-q, C is m-by-q
+ * and S is n-by-q, where q = min(m+n,p).
+ * @param U The U matrix (MxM)
+ * @param V The V matrix (NxN)
+ * @param zeroR The X matrix (rxN)
+ * @param Q The Q matrix (Q)
+ * @param C The C matrix (MxQ)
+ * @param S The S matrix (NxQ)
+ */
+void gsvd(blitz::Array<double,2>& A,
+ blitz::Array<double,2>& B,
+
+ blitz::Array<double,2>& U,
+ blitz::Array<double,2>& V,
+ blitz::Array<double,2>& zeroR,
+ blitz::Array<double,2>& Q,
+ blitz::Array<double,2>& X,
+ blitz::Array<double,2>& C,
+ blitz::Array<double,2>& S
+ );
+ /**
+ * @brief Swaping using the LAPACK variable iWork
+ * http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga4a187519e5c71da3b3f67c85e9baf0f2.html#ga4a187519e5c71da3b3f67c85e9baf0f2
+ *
+ * @param A 1D Matrix
+ * @param indexes Sorting information
+ * @param n number of operations
+ */
+ template<typename T>
+ void swap_(blitz::Array<T,1>& A, int* indexes, int begin, int end) {
+ T aux = 0;
+ int fortran_index = 0;
+ for (int i=0; i<A.extent(0); i++){
+ fortran_index = indexes[i]-1;
+ aux = A(i);
+ A(i) = A(fortran_index);
+ A(fortran_index) = aux;
+ }
+ }
+
+
+ /**
+ * @brief Swaping using the LAPACK variable iWork
+ * http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga4a187519e5c71da3b3f67c85e9baf0f2.html#ga4a187519e5c71da3b3f67c85e9baf0f2
+ *
+ * @param A D Matrix
+ * @param indexes Sorting information
+ * @param n number of operations
+ */
+ template<typename T>
+ void swap_(blitz::Array<T,2>& A, int* indexes, int begin, int end) {
+ blitz::Array<T,1> aux(A.extent(1)); aux = 0;
+ int fortran_index = 0;
+ for (int i=begin; i<end; i++){
+ fortran_index = indexes[i]-1;
+ aux = A(blitz::Range::all(), i);
+ A(blitz::Range::all(), i) = A(blitz::Range::all(), fortran_index);
+ A(blitz::Range::all(), fortran_index) = aux;
+ }
+ }
+
+
+
+
+
+
+/**
+ * @}
+ */
+}}
+
+#endif /* BOB_MATH_GSVD_H */
diff --git a/bob/math/main.cpp b/bob/math/main.cpp
index 4f70e351c17edf0a8e537a44af1b5479eb92aecc..e78fe28de0135d7fdf1e9a78f94f803175404202 100644
--- a/bob/math/main.cpp
+++ b/bob/math/main.cpp
@@ -20,6 +20,7 @@
#include "norminv.h"
#include "scatter.h"
#include "lp_interior_point.h"
+#include "gsvd.h"
static bob::extension::FunctionDoc s_histogram_intersection = bob::extension::FunctionDoc(
"histogram_intersection",
@@ -355,6 +356,26 @@ static bob::extension::FunctionDoc s_scatters_nocheck = bob::extension::Function
;
+static bob::extension::FunctionDoc s_gsvd = bob::extension::FunctionDoc(
+ "gsvd",
+ "Computes the Generalized SVD",
+ "Computes the Generalized SVD. The output of this function is similar with the one found in Matlab\n"
+ "[U,V,X,C,S] = gsvd(A,B) returns unitary matrices :math:`U` and :math:`V`, the square matrix :math:`X` (which is :math:`[0 R] Q^{T}`), and nonnegative diagonal matrices :math:`C` and :math:`S` such that:\n\n"
+ ".. math:: C^{T}C + S^{T}S = I \n"
+ ".. math:: A = (XC^{T}U^{T})^{T}\n"
+ ".. math:: B = (XS^{T}V^{T})^{T}\n"
+ )
+ .add_prototype("A, B", "")
+ .add_parameter("A", "[array_like (float, 2D)]", "Must be :math:`m \\times n`")
+ .add_parameter("B", "[array_like (float, 2D)]", "Must be :math:`p \\times n`")
+ .add_return("U", "[array_like (float, 2D)]", "Contains a :math:`m \\times m` orthogonal matrix.")
+ .add_return("V", "[array_like (float, 2D)]", "Contains a :math:`n \\times n` orthogonal matrix.")
+ .add_return("X", "[array_like (float, 2D)]", "Contains a :math:`p \\times q` matrix, where :math:`p=\\min(m+n,p)` and :math:`X=[0, R] Q^{T}` (Check LAPACK documentation).")
+ .add_return("C", "[array_like (float, 2D)]", "Contains a :math:`p \\times q` matrix.")
+ .add_return("S", "[array_like (float, 2D)]", "Contains a :math:`n \\times q` matrix.")
+;
+
+
static PyMethodDef module_methods[] = {
{
s_histogram_intersection.name(),
@@ -470,6 +491,13 @@ static PyMethodDef module_methods[] = {
METH_VARARGS|METH_KEYWORDS,
s_scatters_nocheck.doc()
},
+ {
+ s_gsvd.name(),
+ (PyCFunction)py_gsvd,
+ METH_VARARGS|METH_KEYWORDS,
+ s_gsvd.doc()
+ },
+
{0} /* Sentinel */
};
diff --git a/bob/math/test_gsvd.py b/bob/math/test_gsvd.py
new file mode 100644
index 0000000000000000000000000000000000000000..1a93b3ca61ec83afe890f6a537f39fe5fe9b7dac
--- /dev/null
+++ b/bob/math/test_gsvd.py
@@ -0,0 +1,67 @@
+#!/usr/bin/env python
+# Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
+# Sun Jan 15 19:12:43 CET 2017
+#
+
+"""
+Tests GSVD
+
+ Basically these tests test the GSVD relation.
+ Given 2 matrices A and B GSVD(A,B) = [U,V,X,C,S] where,
+
+ A= (X * C.T * U^T)^T and
+ B= (X * S.T * V^T)^T and
+
+ C**2 + S**2 = 1
+
+"""
+
+import bob.math
+import numpy
+import nose.tools
+numpy.random.seed(10)
+
+
+def gsvd_relations(A,B):
+ [U,V,X,C,S] = bob.math.gsvd(A, B)
+
+ # Cheking the relation C**2 + S**2 = 1
+ I = numpy.eye(A.shape[1])
+ I_check = numpy.dot(C.T, C) + numpy.dot(S.T, S)
+ nose.tools.eq_( (abs(I-I_check) < 1e-10).all(), True )
+
+ # Cheking the relation A= (X * C.T * U^T)^T
+ A_check = numpy.dot(numpy.dot(X,C.T),U.T).T
+ nose.tools.eq_( (abs(A-A_check) < 1e-10).all(), True )
+
+ # Cheking the relation B= (X * S.T * V^T)^T
+ B_check = numpy.dot(numpy.dot(X,S.T),V.T).T
+ nose.tools.eq_( (abs(B-B_check) < 1e-10).all(), True )
+
+ del U,V,X,C,S
+
+
+def test_first_case():
+ """
+ Testing the first scenario:
+ M-K-L >= 0 (check http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga4a187519e5c71da3b3f67c85e9baf0f2.html#ga4a187519e5c71da3b3f67c85e9baf0f2)
+ """
+
+ A = numpy.random.rand(10,10)
+ B = numpy.random.rand(790,10)
+
+ gsvd_relations(A, B)
+
+
+def test_second_case():
+ """
+ Testing the second scenario:
+ M-K-L < 0 (check http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga4a187519e5c71da3b3f67c85e9baf0f2.html#ga4a187519e5c71da3b3f67c85e9baf0f2)
+ """
+
+ A = numpy.random.rand(4,5)
+ B = numpy.random.rand(11,5)
+
+ gsvd_relations(A, B)
+
+
diff --git a/doc/py_api.rst b/doc/py_api.rst
index aa88ba3b3870ab08c62daa27c09b409451c0dbf9..252233d3e0a009b637dd39cf2de4cc17448f1c99 100644
--- a/doc/py_api.rst
+++ b/doc/py_api.rst
@@ -15,7 +15,19 @@ Summary
bob.math.LPInteriorPoint
bob.math.LPInteriorPointShortstep
bob.math.LPInteriorPointLongstep
-
+ bob.math.chi_square
+ bob.math.gsvd
+ bob.math.histogram_intersection
+ bob.math.kullback_leibler
+ bob.math.linsolve
+ bob.math.linsolve_cg_sympos
+ bob.math.linsolve_sympos
+ bob.math.norminv
+ bob.math.pavx
+ bob.math.pavxWidth
+ bob.math.pavxWidthHeight
+ bob.math.scatter
+ bob.math.scatters
Details
.......
diff --git a/setup.py b/setup.py
index e93ec9260100e1d2ecba7d869f29de8004e20326..8935cc7c6e3b2a1c9a33b089c5438f46c98379ae 100644
--- a/setup.py
+++ b/setup.py
@@ -179,6 +179,7 @@ setup(
"bob/math/cpp/pavx.cpp",
"bob/math/cpp/pinv.cpp",
"bob/math/cpp/svd.cpp",
+ "bob/math/cpp/gsvd.cpp",
"bob/math/cpp/sqrtm.cpp",
],
version = version,
@@ -195,6 +196,7 @@ setup(
"bob/math/linsolve.cpp",
"bob/math/pavx.cpp",
"bob/math/norminv.cpp",
+ "bob/math/gsvd.cpp",
"bob/math/scatter.cpp",
"bob/math/lp_interior_point.cpp",
"bob/math/main.cpp",