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

First version of updated linsolve; seems to have an issue

parent 50d5c18e
......@@ -22,8 +22,8 @@ extern "C" void dgesv_( const int *N, const int *NRHS, double *A,
extern "C" void dposv_( const char* uplo, const int *N, const int *NRHS,
double *A, const int *lda, double *B, const int *ldb, int *info);
void bob::math::linsolve(const blitz::Array<double,2>& A, blitz::Array<double,1>& x,
const blitz::Array<double,1>& b)
void bob::math::linsolve(const blitz::Array<double,2>& A,
const blitz::Array<double,1>& b, blitz::Array<double,1>& x)
{
// Check x and b
bob::core::array::assertZeroBase(x);
......@@ -35,12 +35,6 @@ void bob::math::linsolve(const blitz::Array<double,2>& A, blitz::Array<double,1>
bob::core::array::assertSameDimensionLength(A.extent(0), A.extent(1));
bob::core::array::assertSameDimensionLength(A.extent(1), b.extent(0));
bob::math::linsolve_(A, x, b);
}
void bob::math::linsolve_(const blitz::Array<double,2>& A, blitz::Array<double,1>& x,
const blitz::Array<double,1>& b)
{
// Defines dimensionality variables
const int N = A.extent(0);
......@@ -82,8 +76,8 @@ void bob::math::linsolve_(const blitz::Array<double,2>& A, blitz::Array<double,1
}
void bob::math::linsolve(const blitz::Array<double,2>& A, blitz::Array<double,2>& X,
const blitz::Array<double,2>& B)
void bob::math::linsolve(const blitz::Array<double,2>& A,
const blitz::Array<double,2>& B, blitz::Array<double,2>& X)
{
// Checks dimensionality and zero base
bob::core::array::assertZeroBase(A);
......@@ -94,12 +88,6 @@ void bob::math::linsolve(const blitz::Array<double,2>& A, blitz::Array<double,2>
bob::core::array::assertSameDimensionLength(A.extent(0), B.extent(0));
bob::core::array::assertSameDimensionLength(X.extent(1), B.extent(1));
bob::math::linsolve_(A, X, B);
}
void bob::math::linsolve_(const blitz::Array<double,2>& A, blitz::Array<double,2>& X,
const blitz::Array<double,2>& B)
{
// Defines dimensionality variables
const int N = A.extent(0);
const int P = X.extent(1);
......@@ -148,7 +136,7 @@ void bob::math::linsolve_(const blitz::Array<double,2>& A, blitz::Array<double,2
void bob::math::linsolveSympos(const blitz::Array<double,2>& A,
blitz::Array<double,1>& x, const blitz::Array<double,1>& b)
const blitz::Array<double,1>& b, blitz::Array<double,1>& x)
{
// Check x and b
bob::core::array::assertZeroBase(x);
......@@ -160,12 +148,6 @@ 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, x, b);
}
void bob::math::linsolveSympos_(const blitz::Array<double,2>& A,
blitz::Array<double,1>& x, const blitz::Array<double,1>& b)
{
// Defines dimensionality variables
const int N = A.extent(0);
......@@ -208,8 +190,8 @@ void bob::math::linsolveSympos_(const blitz::Array<double,2>& A,
x = x_blitz_lapack;
}
void bob::math::linsolveSympos(const blitz::Array<double,2>& A, blitz::Array<double,2>& X,
const blitz::Array<double,2>& B)
void bob::math::linsolveSympos(const blitz::Array<double,2>& A,
const blitz::Array<double,2>& B, blitz::Array<double,2>& X)
{
// Checks dimensionality and zero base
bob::core::array::assertZeroBase(A);
......@@ -220,12 +202,6 @@ void bob::math::linsolveSympos(const blitz::Array<double,2>& A, blitz::Array<dou
bob::core::array::assertSameDimensionLength(A.extent(0), B.extent(0));
bob::core::array::assertSameDimensionLength(X.extent(1), B.extent(1));
bob::math::linsolveSympos_(A, X, B);
}
void bob::math::linsolveSympos_(const blitz::Array<double,2>& A, blitz::Array<double,2>& X,
const blitz::Array<double,2>& B)
{
// Defines dimensionality variables
const int N = A.extent(0);
const int P = X.extent(1);
......@@ -275,8 +251,9 @@ void bob::math::linsolveSympos_(const blitz::Array<double,2>& A, blitz::Array<do
void bob::math::linsolveCGSympos(const blitz::Array<double,2>& A, blitz::Array<double,1>& x,
const blitz::Array<double,1>& b, const double acc, const int 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);
......@@ -291,15 +268,6 @@ void bob::math::linsolveCGSympos(const blitz::Array<double,2>& A, blitz::Array<d
bob::core::array::assertSameDimensionLength(A.extent(0), N);
bob::core::array::assertSameDimensionLength(A.extent(1), N);
bob::math::linsolveCGSympos_(A, x, b, acc, max_iter);
}
void bob::math::linsolveCGSympos_(const blitz::Array<double,2>& A, blitz::Array<double,1>& x,
const blitz::Array<double,1>& b, 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;
......@@ -347,4 +315,3 @@ void bob::math::linsolveCGSympos_(const blitz::Array<double,2>& A, blitz::Array<
// TODO return best_res and number of iterations?
//double res = best_res;
}
......@@ -9,7 +9,7 @@
#define BOB_MATH_CONFIG_H
/* Macros that define versions and important names */
#define BOB_MATH_API_VERSION 0x0200
#define BOB_MATH_API_VERSION 0x0201
#ifdef BOB_IMPORT_VERSION
......
......@@ -23,10 +23,8 @@ namespace bob { namespace math {
* @param x The x vector of the system A*x=b which will be updated
* at the end of the function.
*/
void linsolve(const blitz::Array<double,2>& A, blitz::Array<double,1>& x,
const blitz::Array<double,1>& b);
void linsolve_(const blitz::Array<double,2>& A, blitz::Array<double,1>& x,
const blitz::Array<double,1>& b);
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,10 +34,8 @@ namespace bob { namespace math {
* @param X The X matrix of the system A*X=B which will be updated
* at the end of the function (size NxP).
*/
void linsolve(const blitz::Array<double,2>& A, blitz::Array<double,2>& X,
const blitz::Array<double,2>& B);
void linsolve_(const blitz::Array<double,2>& A, blitz::Array<double,2>& X,
const blitz::Array<double,2>& B);
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
......@@ -53,9 +49,7 @@ namespace bob { namespace math {
* at the end of the function (size N)
*/
void linsolveSympos(const blitz::Array<double,2>& A,
blitz::Array<double,1>& x, const blitz::Array<double,1>& b);
void linsolveSympos_(const blitz::Array<double,2>& A,
blitz::Array<double,1>& x, const blitz::Array<double,1>& b);
const blitz::Array<double,1>& b, blitz::Array<double,1>& x);
/**
* @brief Function which solves a symmetric positive definite linear
......@@ -69,9 +63,7 @@ namespace bob { namespace math {
* at the end of the function (size NxP)
*/
void linsolveSympos(const blitz::Array<double,2>& A,
blitz::Array<double,2>& X, const blitz::Array<double,2>& B);
void linsolveSympos_(const blitz::Array<double,2>& A,
blitz::Array<double,2>& X, const blitz::Array<double,2>& B);
const blitz::Array<double,2>& B, blitz::Array<double,2>& X);
/**
......@@ -86,12 +78,10 @@ namespace bob { namespace math {
* norm(Ax-b)/norm(b) < acc
* @param max_iter The maximum number of iterations
*/
void linsolveCGSympos(const blitz::Array<double,2>& A, blitz::Array<double,1>& x,
const blitz::Array<double,1>& b, const double acc, const int max_iter);
void linsolveCGSympos_(const blitz::Array<double,2>& A, blitz::Array<double,1>& x,
const blitz::Array<double,1>& b, 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);
}}
#endif /* BOB_MATH_LINSOLVE_H */
This diff is collapsed.
/**
* @author Andre Anjos <andre.anjos@idiap.ch>
* @date Wed 4 Dec 15:26:54 2013
* @date Wed 4 Dec 15:26:54 2013
*
* @brief Declaration of components relevant for main.cpp
*/
......@@ -8,8 +8,9 @@
#include <Python.h>
PyObject* py_linsolve(PyObject*, PyObject* args, PyObject* kwargs);
PyObject* py_linsolve_nocheck(PyObject*, PyObject* args, PyObject* kwargs);
PyObject* py_linsolve_sympos(PyObject*, PyObject* args, PyObject* kwargs);
PyObject* py_linsolve_sympos_nocheck(PyObject*, PyObject* args, PyObject* kwargs);
PyObject* py_linsolve_cg_sympos(PyObject*, PyObject* args, PyObject* kwargs);
PyObject* py_linsolve_cg_sympos_nocheck(PyObject*, PyObject* args, PyObject* kwargs);
extern bob::extension::FunctionDoc s_linsolve;
extern bob::extension::FunctionDoc s_linsolve_sympos;
extern bob::extension::FunctionDoc s_linsolve_cg_sympos;
......@@ -82,107 +82,9 @@ static bob::extension::FunctionDoc s_kullback_leibler = bob::extension::Function
.add_return("dist", "float", "The Kullback-Leibler divergence value for the given histograms.")
;
static bob::extension::FunctionDoc s_linsolve = bob::extension::FunctionDoc(
"linsolve",
"Solves the linear system :math:`Ax=b` and returns the result in :math:`x`.",
"This method uses LAPACK's ``dgesv`` generic solver. "
"You can use this method in two different formats. "
"The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
"The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
)
.add_prototype("A, b", "x")
.add_prototype("A, b, x")
.add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
.add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
.add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
.add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
;
static bob::extension::FunctionDoc s_linsolve_nocheck = bob::extension::FunctionDoc(
"linsolve_",
"Solves the linear system :math:`Ax=b` and returns the result in :math:`x`.",
".. warning:: This variant does not perform any checks on the input matrices and is faster then :py:func:`linsolve`. "
"Use it when you are sure your input matrices sizes match.\n\n"
"This method uses LAPACK's ``dgesv`` generic solver. "
"You can use this method in two different formats. "
"The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
"The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
)
.add_prototype("A, b", "x")
.add_prototype("A, b, x")
.add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
.add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
.add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
.add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
;
static bob::extension::FunctionDoc s_linsolve_sympos = bob::extension::FunctionDoc(
"linsolve_sympos",
"Solves the linear system :math:`Ax=b` and returns the result in :math:`x` for symmetric :math:`A` matrix.",
"This method uses LAPACK's ``dposv`` solver, assuming :math:`A` is a symmetric positive definite matrix. "
"You can use this method in two different formats. "
"The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
"The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
)
.add_prototype("A, b", "x")
.add_prototype("A, b, x")
.add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
.add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
.add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
.add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
;
static bob::extension::FunctionDoc s_linsolve_sympos_nocheck = bob::extension::FunctionDoc(
"linsolve_sympos_",
"Solves the linear system :math:`Ax=b` and returns the result in :math:`x` for symmetric :math:`A` matrix.",
".. warning:: This variant does not perform any checks on the input matrices and is faster then :py:func:`linsolve_sympos`. "
"Use it when you are sure your input matrices sizes match.\n\n"
"This method uses LAPACK's ``dposv`` solver, assuming :math:`A` is a symmetric positive definite matrix. "
"You can use this method in two different formats. "
"The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
"The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
)
.add_prototype("A, b", "x")
.add_prototype("A, b, x")
.add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
.add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
.add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
.add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
;
static bob::extension::FunctionDoc s_linsolve_cg_sympos = bob::extension::FunctionDoc(
"linsolve_cg_sympos",
"Solves the linear system :math:`Ax=b` using conjugate gradients and returns the result in :math:`x` for symmetric :math:`A` matrix.",
"This method uses the conjugate gradient solver, assuming :math:`A` is a symmetric positive definite matrix. "
"You can use this method in two different formats. "
"The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
"The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
)
.add_prototype("A, b", "x")
.add_prototype("A, b, x")
.add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
.add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
.add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
.add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
;
static bob::extension::FunctionDoc s_linsolve_cg_sympos_nocheck = bob::extension::FunctionDoc(
"linsolve_cg_sympos_",
"Solves the linear system :math:`Ax=b` using conjugate gradients and returns the result in :math:`x` for symmetric :math:`A` matrix.",
".. warning:: This variant does not perform any checks on the input matrices and is faster then :py:func:`linsolve_cg_sympos`. "
"Use it when you are sure your input matrices sizes match.\n\n"
"This method uses the conjugate gradient solver, assuming :math:`A` is a symmetric positive definite matrix. "
"You can use this method in two different formats. "
"The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
"The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
)
.add_prototype("A, b", "x")
.add_prototype("A, b, x")
.add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
.add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
.add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
.add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
;
static bob::extension::FunctionDoc s_pavx = bob::extension::FunctionDoc(
"pavx",
......@@ -417,36 +319,18 @@ static PyMethodDef module_methods[] = {
METH_VARARGS|METH_KEYWORDS,
s_linsolve.doc()
},
{
s_linsolve_nocheck.name(),
(PyCFunction)py_linsolve_nocheck,
METH_VARARGS|METH_KEYWORDS,
s_linsolve_nocheck.doc()
},
{
s_linsolve_sympos.name(),
(PyCFunction)py_linsolve_sympos,
METH_VARARGS|METH_KEYWORDS,
s_linsolve_sympos.doc()
},
{
s_linsolve_sympos_nocheck.name(),
(PyCFunction)py_linsolve_sympos_nocheck,
METH_VARARGS|METH_KEYWORDS,
s_linsolve_sympos_nocheck.doc()
},
{
s_linsolve_cg_sympos.name(),
(PyCFunction)py_linsolve_cg_sympos,
METH_VARARGS|METH_KEYWORDS,
s_linsolve_cg_sympos.doc()
},
{
s_linsolve_cg_sympos_nocheck.name(),
(PyCFunction)py_linsolve_cg_sympos_nocheck,
METH_VARARGS|METH_KEYWORDS,
s_linsolve_cg_sympos_nocheck.doc()
},
{
s_pavx.name(),
(PyCFunction)py_pavx,
......
......@@ -30,7 +30,7 @@ def test_linsolve():
x1 = numpy.ndarray((3,), 'float64')
# Computes the solution
linsolve(A,x1,b)
linsolve(A,b,x1)
x2 = linsolve(A,b)
# Compare to reference
......@@ -54,7 +54,7 @@ def test_linsolveSympos():
x1 = numpy.ndarray((3,), 'float64')
# Computes the solution
linsolve_sympos(A,x1,b)
linsolve_sympos(A,b,x1)
x2 = linsolve_sympos(A,b)
# Compare to reference
......@@ -81,7 +81,7 @@ def test_linsolveCGSympos():
# Computes the solution
eps = 1e-6
max_iter = 1000
linsolve_cg_sympos(A,x1,b,eps,max_iter)
linsolve_cg_sympos(A,b,x1,eps,max_iter)
x2 = linsolve_cg_sympos(A,b,eps,max_iter)
# Compare to reference
......
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