Subject: [Boost-commit] svn:boost r50176 - sandbox/boost/numeric/bindings/lapack
From: thomas.klimpel_at_[hidden]
Date: 2008-12-07 08:19:55
Author: klimpel
Date: 2008-12-07 08:19:55 EST (Sun, 07 Dec 2008)
New Revision: 50176
URL: http://svn.boost.org/trac/boost/changeset/50176
Log:
added Jesse Manning's gels bindings
Added:
sandbox/boost/numeric/bindings/lapack/gels.hpp (contents, props changed)
sandbox/boost/numeric/bindings/lapack/gelsd.hpp (contents, props changed)
sandbox/boost/numeric/bindings/lapack/gelss.hpp (contents, props changed)
Added: sandbox/boost/numeric/bindings/lapack/gels.hpp
==============================================================================
--- (empty file)
+++ sandbox/boost/numeric/bindings/lapack/gels.hpp 2008-12-07 08:19:55 EST (Sun, 07 Dec 2008)
@@ -0,0 +1,209 @@
+//
+// Copyright Jesse Manning 2007
+//
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+//
+
+#ifndef BOOST_NUMERIC_BINDINGS_LAPACK_GELS_HPP
+#define BOOST_NUMERIC_BINDINGS_LAPACK_GELS_HPP
+
+#include <algorithm>
+
+#include <boost/numeric/bindings/traits/type.hpp>
+#include <boost/numeric/bindings/traits/traits.hpp>
+#include <boost/numeric/bindings/traits/type_traits.hpp>
+#include <boost/numeric/bindings/lapack/lapack.h>
+#include <boost/numeric/bindings/lapack/workspace.hpp>
+#include <boost/numeric/bindings/traits/detail/array.hpp>
+#include <boost/numeric/bindings/traits/detail/utils.hpp>
+
+
+#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
+# include <boost/static_assert.hpp>
+# include <boost/type_traits.hpp>
+#endif
+
+namespace boost { namespace numeric { namespace bindings {
+
+ namespace lapack {
+
+ ////////////////////////////////////////////////////////////////////////
+ //
+ // Linear Least Squares of an underdetermined or overdetermined matrix
+ //
+ ///////////////////////////////////////////////////////////////////////
+
+ /* gels - uses the LQ or QR factorization to solve an overdetermined
+ * or underdetermined linear system. A full rank matrix is
+ * assumed.
+ *
+ * The linear least squares system is defined by A*x=b. A is the m-by-n
+ * coefficients matrix and b is the m-by-nrhs matrix. Several
+ * right hand side vectors b and solution vectors x can be handled in
+ * a single call; they are stored as the columns of the m-by-nrhs right
+ * hand side matrix B and the n-by-nrh solution matrix x.
+ *
+ * If trans = 'N' and m >= n: find least squares solution of overdetermined system
+ * minimizes || b - A x ||2
+ *
+ * if trans = 'N' and m < n: find minimum norm solution of underdetermined system
+ * A*X = B
+ *
+ * if trans = 'T' or 'C' and m >= n: find minimum norm soltion of underdetermined system
+ * A^H*X = B
+ *
+ * if trans = 'T' or 'C' and m < n: find least squares solution of overdetermined system
+ * minimize || b - A^H x ||2
+ *
+ * Workspace is organized following the arguments in the calling sequence.
+ * optimal_workspace() : for optimizing use of blas 3 kernels
+ * minimal_workspace() : minimum size of work arrays, but does not allow for optimization
+ * of blas 3 kernels
+ * workspace( work ) : where work size must be at least min(m, n) + max(1, m, n, nrhs)
+ *
+ *
+ * Each solution vector x (length n) is returned column-wise in the rhs matrix B.
+ */
+
+ namespace detail {
+
+ inline void gels(char const trans, const int m, const int n,
+ const int nrhs, float *a, const int lda,
+ float *b, const int ldb, float *work,
+ const int lwork, int *info)
+ {
+ LAPACK_SGELS(&trans, &m, &n, &nrhs, a, &lda, b, &ldb, work, &lwork, info);
+ }
+
+ inline void gels(char const trans, const int m, const int n,
+ const int nrhs, double *a, const int lda,
+ double *b, const int ldb, double *work,
+ const int lwork, int *info)
+ {
+ LAPACK_DGELS(&trans, &m, &n, &nrhs, a, &lda, b, &ldb, work, &lwork, info);
+ }
+
+ inline void gels(char const trans, const int m, const int n,
+ const int nrhs, traits::complex_f *a, const int lda,
+ traits::complex_f *b, const int ldb, traits::complex_f *work,
+ const int lwork, int *info)
+ {
+ LAPACK_CGELS(&trans, &m, &n, &nrhs,
+ traits::complex_ptr(a), &lda,
+ traits::complex_ptr(b), &ldb,
+ traits::complex_ptr(work), &lwork, info);
+ }
+
+ inline void gels(char const trans, const int m, const int n,
+ const int nrhs, traits::complex_d *a, const int lda,
+ traits::complex_d *b, const int ldb, traits::complex_d *work,
+ const int lwork, int *info)
+ {
+ LAPACK_ZGELS(&trans, &m, &n, &nrhs,
+ traits::complex_ptr(a), &lda,
+ traits::complex_ptr(b), &ldb,
+ traits::complex_ptr(work), &lwork, info);
+ }
+
+ // generic function that calls more detailed lapack function
+ template <typename MatrA, typename VecB, typename Work>
+ int gels(const char trans, MatrA& A, VecB& b, Work& work)
+ {
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int mrhs = traits::matrix_size1(b);
+ const int nrhs = traits::matrix_size2(b);
+
+ // sanity checks
+ assert(trans == 'N' || trans == 'T' || trans == 'C');
+ assert(m >= 0 && n >= 0);
+ assert(mrhs >= 0 && nrhs >= 0);
+ assert(traits::leading_dimension(A) >= 1 && traits::leading_dimension(b) >= 1);
+ assert(mrhs == std::max(m, n));
+
+ int info;
+ detail::gels(trans,
+ traits::matrix_size1(A),
+ traits::matrix_size2(A),
+ traits::matrix_size2(b),
+ traits::matrix_storage(A),
+ traits::leading_dimension(A),
+ traits::matrix_storage(b),
+ traits::leading_dimension(b),
+ traits::vector_storage(work),
+ traits::vector_size(work),
+ &info);
+
+ return info;
+ }
+
+ // query for recommended workspace
+ template <typename MatrA, typename VecB>
+ inline
+ int gels_optimal_work(const char trans, MatrA& A, VecB& b)
+ {
+ typename MatrA::value_type work;
+ int info;
+ detail::gels(trans,
+ traits::matrix_size1(A),
+ traits::matrix_size2(A),
+ traits::matrix_size2(b),
+ traits::matrix_storage(A),
+ traits::leading_dimension(A),
+ traits::matrix_storage(b),
+ traits::leading_dimension(b),
+ &work, //traits::vector_storage(work),
+ -1, // traits::vector_size(work),
+ &info);
+
+ assert(info == 0);
+
+ int lwork = traits::detail::to_int(work);
+
+ return lwork;
+ }
+ } // namespace detail
+
+
+ template <typename MatrA, typename VecB>
+ inline
+ int gels(const char trans, MatrA& A, VecB& b, optimal_workspace)
+ {
+ // query optimal workspace size
+ int work_size = detail::gels_optimal_work(trans, A, b);
+ traits::detail::array<typename MatrA::value_type> work(work_size);
+
+ return detail::gels(trans, A, b, work);
+ }
+
+ template <typename MatrA, typename VecB>
+ inline
+ int gels(const char trans, MatrA& A, VecB& b, minimal_workspace)
+ {
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int r = traits::matrix_size2(b);
+
+ const int minmn = std::min(m, n); //m < n ? m : n;
+ const int maxmn = std::max(m, n); // m > n ? m : n;
+ const int maxdim = std::max(maxmn, r); // maxmn > r ? maxmn : r;
+
+ traits::detail::array<typename MatrA::value_type> work(minmn + std::max(1, maxdim));
+
+ return detail::gels(trans, A, b, work);
+ }
+
+ template <typename MatrA, typename VecB, typename Work>
+ inline
+ int gels(const char trans, MatrA& A, VecB& b, detail::workspace1<Work> workspace)
+ {
+ return detail::gels(trans, A, b, workspace.w_);
+ }
+
+ } // namespace lapack
+
+}}}
+
+#endif
Added: sandbox/boost/numeric/bindings/lapack/gelsd.hpp
==============================================================================
--- (empty file)
+++ sandbox/boost/numeric/bindings/lapack/gelsd.hpp 2008-12-07 08:19:55 EST (Sun, 07 Dec 2008)
@@ -0,0 +1,416 @@
+//
+// Copyright Jesse Manning 2007
+//
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+//
+
+#ifndef BOOST_NUMERIC_BINDINGS_LAPACK_GELSD_HPP
+#define BOOST_NUMERIC_BINDINGS_LAPACK_GELSD_HPP
+
+#include <algorithm>
+
+#include <boost/numeric/bindings/traits/type.hpp>
+#include <boost/numeric/bindings/traits/traits.hpp>
+#include <boost/numeric/bindings/traits/type_traits.hpp>
+#include <boost/numeric/bindings/lapack/lapack.h>
+#include <boost/numeric/bindings/lapack/workspace.hpp>
+#include <boost/numeric/bindings/traits/detail/array.hpp>
+#include <boost/numeric/bindings/traits/detail/utils.hpp>
+#include <boost/numeric/bindings/lapack/ilaenv.hpp>
+
+
+#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
+# include <boost/static_assert.hpp>
+# include <boost/type_traits.hpp>
+#endif
+
+namespace boost { namespace numeric { namespace bindings {
+
+ namespace lapack {
+
+ namespace detail {
+
+ inline void gelsd(const int m, const int n, const int nrhs,
+ float *a, const int lda, float *b, const int ldb,
+ float *s, const float rcond, int *rank, float *work,
+ const int lwork, int *iwork, int *info)
+ {
+ LAPACK_SGELSD(&m, &n, &nrhs, a, &lda, b, &ldb, s,
+ &rcond, rank, work, &lwork, iwork, info);
+ }
+
+ inline void gelsd(const int m, const int n, const int nrhs,
+ double *a, const int lda, double *b, const int ldb,
+ double *s, const double rcond, int *rank, double *work,
+ const int lwork, int *iwork, int *info)
+ {
+ LAPACK_DGELSD(&m, &n, &nrhs, a, &lda, b, &ldb, s,
+ &rcond, rank, work, &lwork, iwork, info);
+ }
+
+ inline void gelsd(const int m, const int n, const int nrhs,
+ traits::complex_f *a, const int lda, traits::complex_f *b,
+ const int ldb, float *s, const float rcond, int *rank,
+ traits::complex_f *work, const int lwork, float *rwork,
+ int *iwork, int *info)
+ {
+ LAPACK_CGELSD(&m, &n, &nrhs, traits::complex_ptr(a),
+ &lda, traits::complex_ptr(b), &ldb, s,
+ &rcond, rank, traits::complex_ptr(work),
+ &lwork, rwork, iwork, info);
+ }
+
+ inline void gelsd(const int m, const int n, const int nrhs,
+ traits::complex_d *a, const int lda, traits::complex_d *b,
+ const int ldb, double *s, const double rcond, int *rank,
+ traits::complex_d *work, const int lwork, double *rwork,
+ int *iwork, int *info)
+ {
+ LAPACK_ZGELSD(&m, &n, &nrhs, traits::complex_ptr(a),
+ &lda, traits::complex_ptr(b), &ldb, s,
+ &rcond, rank, traits::complex_ptr(work),
+ &lwork, rwork, iwork, info);
+ }
+
+ // gelsd for real type
+ template <typename MatrA, typename MatrB, typename VecS, typename Work>
+ inline int gelsd(MatrA& A, MatrB& B, VecS& s, Work& work)
+ {
+ typedef typename MatrA::value_type val_t;
+ typedef typename traits::type_traits<val_t>::real_type real_t;
+
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int nrhs = traits::matrix_size2(B);
+ const int maxmn = std::max(m, n);
+ const int minmn = std::min(m, n);
+
+ // sanity checks
+ assert(m >= 0 && n >= 0);
+ assert(nrhs >= 0);
+ assert(traits::leading_dimension(A) >= std::max(1, m));
+ assert(traits::leading_dimension(B) >= std::max(1, maxmn));
+ assert(traits::vector_size(work) >= 1);
+ assert(traits::vector_size(s) >= std::max(1, minmn));
+
+ int info;
+ const real_t rcond = -1; // use machine precision
+ int rank;
+
+ // query for maximum size of subproblems
+ const int smlsiz = ilaenv(9, "GELSD", "");
+ const int nlvl = static_cast<int>(((std::log(static_cast<float>(minmn))/std::log(2.f))/ (smlsiz+1)) + 1);
+
+ traits::detail::array<int> iwork(3*minmn*nlvl + 11*minmn);
+
+ detail::gelsd(traits::matrix_size1(A),
+ traits::matrix_size2(A),
+ traits::matrix_size2(B),
+ traits::matrix_storage(A),
+ traits::leading_dimension(A),
+ traits::matrix_storage(B),
+ traits::leading_dimension(B),
+ traits::vector_storage(s),
+ rcond,
+ &rank,
+ traits::vector_storage(work),
+ traits::vector_size(work),
+ traits::vector_storage(iwork),
+ &info);
+
+ return info;
+ }
+
+ // gelsd for complex type
+ template <typename MatrA, typename MatrB, typename VecS,
+ typename Work, typename RWork>
+ inline int gelsd(MatrA& A, MatrB& B, VecS& s, Work& work, RWork& rwork)
+ {
+ typedef typename MatrA::value_type val_t;
+ typedef typename traits::type_traits<val_t>::real_type real_t;
+
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int nrhs = traits::matrix_size2(B);
+ const int maxmn = std::max(m, n);
+ const int minmn = std::min(m, n);
+
+ // sanity checks
+ assert(m >= 0 && n >= 0);
+ assert(nrhs >= 0);
+ assert(traits::leading_dimension(A) >= std::max(1, m));
+ assert(traits::leading_dimension(B) >= std::max(1, maxmn));
+ assert(traits::vector_size(work) >= 1);
+ assert(traits::vector_size(s) >= std::max(1, minmn));
+
+ int info;
+ const real_t rcond = -1; // use machine precision
+ int rank;
+
+ // query for maximum size of subproblems
+ const int smlsiz = ilaenv(9, "GELSD", "");
+ const int nlvl = static_cast<int>(((std::log(static_cast<float>(minmn))/std::log(2.f))/ (smlsiz+1)) + 1);
+
+ traits::detail::array<int> iwork(3*minmn*nlvl + 11*minmn);
+
+ detail::gelsd(traits::matrix_size1(A),
+ traits::matrix_size2(A),
+ traits::matrix_size2(B),
+ traits::matrix_storage(A),
+ traits::leading_dimension(A),
+ traits::matrix_storage(B),
+ traits::leading_dimension(B),
+ traits::vector_storage(s),
+ rcond,
+ &rank,
+ traits::vector_storage(work),
+ traits::vector_size(work),
+ traits::vector_storage(rwork),
+ traits::vector_storage(iwork),
+ &info);
+
+ return info;
+ }
+
+ template <int N>
+ struct Gelsd { };
+
+ // specialization for gelsd real flavors (sgelsd, dgelsd)
+ template <>
+ struct Gelsd<1>
+ {
+ template <typename MatrA, typename MatrB, typename VecS>
+ inline int operator() (MatrA& A, MatrB& B, VecS& s, minimal_workspace) const
+ {
+ typedef typename MatrA::value_type val_t;
+
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int nrhs = traits::matrix_size2(B);
+
+ const int minmn = std::min(m, n); // minmn = m < n ? m : n
+ const int maxmn = std::max(m, n); // maxmn = m > n ? m : n
+ const int maxmnr = std::max(maxmn, nrhs); // maxmnr = maxmn > nrhs ? maxmn : nrhs
+
+ // query for maximum size of subproblems
+ const int smlsiz = ilaenv(9, "GELSD", "");
+ const int nlvl = static_cast<int>(((std::log(static_cast<float>(minmn))/std::log(2.f))/ (smlsiz+1)) + 1);
+
+ const int lwork = 12*minmn + 2*minmn*smlsiz + 8*minmn*nlvl +
+ minmn*nrhs + (smlsiz+1)*(smlsiz+1);
+
+ traits::detail::array<val_t> work(lwork);
+
+ return gelsd(A, B, s, work);
+ }
+
+ template <typename MatrA, typename MatrB, typename VecS>
+ inline int operator() (MatrA& A, MatrB& B, VecS& s, optimal_workspace) const
+ {
+ typedef typename MatrA::value_type val_t;
+ typedef typename traits::type_traits<val_t>::real_type real_t;
+
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int nrhs = traits::matrix_size2(B);
+
+ const int minmn = std::min(m, n); // minmn = m < n ? m : n
+ const int maxmn = std::max(m, n); // maxmn = m > n ? m : n
+ const int maxmnr = std::max(maxmn, nrhs); // maxmnr = maxmn > nrhs ? maxmn : nrhs
+
+ val_t temp_work;
+ int temp_iwork;
+
+ const real_t rcond = -1;
+ int rank;
+ int info;
+
+ // query for optimal workspace size
+ detail::gelsd(traits::matrix_size1(A),
+ traits::matrix_size2(A),
+ traits::matrix_size2(B),
+ traits::matrix_storage(A),
+ traits::leading_dimension(A),
+ traits::matrix_storage(B),
+ traits::leading_dimension(B),
+ traits::vector_storage(s),
+ rcond,
+ &rank,
+ &temp_work, //traits::vector_storage(work),
+ -1, //traits::vector_size(work),
+ &temp_iwork,
+ &info);
+
+ assert(info == 0);
+
+ const int lwork = traits::detail::to_int(temp_work);
+
+ traits::detail::array<val_t> work(lwork);
+
+ return gelsd(A, B, s, work);
+ }
+
+ template <typename MatrA, typename MatrB, typename VecS, typename Work>
+ inline int operator() (MatrA& A, MatrB& B, VecS& s, detail::workspace1<Work> workspace) const
+ {
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int nrhs = traits::matrix_size2(B);
+
+ const int minmn = std::min(m, n); // minmn = m < n ? m : n
+ const int maxmn = std::max(m, n); // maxmn = m > n ? m : n
+ const int maxmnr = std::max(maxmn, nrhs); // maxmnr = maxmn > nrhs ? maxmn : nrhs
+
+ return gelsd(A, B, s, workspace.w_);
+ }
+ };
+
+ // specialization for gelsd (cgelsd, zgelsd)
+ template <>
+ struct Gelsd<2>
+ {
+ template <typename MatrA, typename MatrB, typename VecS>
+ inline int operator() (MatrA& A, MatrB& B, VecS& s, minimal_workspace) const
+ {
+ typedef typename MatrA::value_type val_t;
+ typedef typename traits::type_traits<val_t>::real_type real_t;
+
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int nrhs = traits::matrix_size2(B);
+
+ const int minmn = std::min(m, n); // minmn = m < n ? m : n
+ const int maxmn = std::max(m, n); // maxmn = m > n ? m : n
+ const int maxmnr = std::max(maxmn, nrhs); // maxmnr = maxmn > nrhs ? maxmn : nrhs
+
+ // query for maximum size of subproblems
+ const int smlsiz = ilaenv(9, "GELSD", "");
+ const int nlvl = static_cast<int>(((std::log(static_cast<float>(minmn))/std::log(2.f))/ (smlsiz+1)) + 1);
+
+ traits::detail::array<val_t> work(2*minmn + minmn*nrhs);
+
+ const int rwork_size = 10*minmn + 2*minmn*smlsiz + 8*minmn*nlvl +
+ 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
+
+ traits::detail::array<real_t> rwork(std::max(1, rwork_size));
+
+ return gelsd(A, B, s, work, rwork);
+ }
+
+ template <typename MatrA, typename MatrB, typename VecS>
+ inline int operator() (MatrA& A, MatrB& B, VecS& s, optimal_workspace) const
+ {
+ typedef typename MatrA::value_type val_t;
+ typedef typename traits::type_traits<val_t>::real_type real_t;
+
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int nrhs = traits::matrix_size2(B);
+
+ const int minmn = std::min(m, n); // minmn = m < n ? m : n
+ const int maxmn = std::max(m, n); // maxmn = m > n ? m : n
+ const int maxmnr = std::max(maxmn, nrhs); // maxmnr = maxmn > nrhs ? maxmn : nrhs
+
+ val_t temp_work;
+ real_t temp_rwork;
+ int temp_iwork;
+
+ const real_t rcond = -1;
+ int rank;
+ int info;
+
+ // query for optimal workspace size
+ detail::gelsd(traits::matrix_size1(A),
+ traits::matrix_size2(A),
+ traits::matrix_size2(B),
+ traits::matrix_storage(A),
+ traits::leading_dimension(A),
+ traits::matrix_storage(B),
+ traits::leading_dimension(B),
+ traits::vector_storage(s),
+ rcond,
+ &rank,
+ &temp_work, //traits::vector_storage(work),
+ -1, //traits::vector_size(work),
+ &temp_rwork,
+ &temp_iwork,
+ &info);
+
+ assert(info == 0);
+
+ const int lwork = traits::detail::to_int(temp_work);
+
+ traits::detail::array<val_t> work(lwork);
+
+ // query for maximum size of subproblems
+ const int smlsiz = ilaenv(9, "GELSD", "");
+ const int nlvl = static_cast<int>(((std::log(static_cast<float>(minmn))/std::log(2.f))/ (smlsiz+1)) + 1);
+
+ const int rwork_size = 10*minmn + 2*minmn*smlsiz + 8*minmn*nlvl +
+ 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
+
+ traits::detail::array<real_t> rwork(std::max(1, rwork_size));
+
+ return gelsd(A, B, s, work, rwork);
+ }
+
+ template <typename MatrA, typename MatrB, typename VecS, typename Work, typename RWork>
+ inline int operator() (MatrA& A, MatrB& B, VecS& s, detail::workspace2<Work, RWork> workspace) const
+ {
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int nrhs = traits::matrix_size2(B);
+
+ const int minmn = std::min(m, n); // minmn = m < n ? m : n
+ const int maxmn = std::max(m, n); // maxmn = m > n ? m : n
+ const int maxmnr = std::max(maxmn, nrhs); // maxmnr = maxmn > nrhs ? maxmn : nrhs
+
+ return gelsd(A, B, s, workspace.w_, workspace.wr_);
+ }
+ };
+
+ } // detail
+
+ // gelsd
+ // Parameters:
+ // A: matrix of coefficients
+ // B: matrix of solutions (stored column-wise)
+ // s: vector to store singular values on output, length >= max(1, min(m,n))
+ // workspace: either optimal, minimal, or user supplied
+ //
+ template <typename MatrA, typename MatrB, typename VecS, typename Work>
+ inline int gelsd(MatrA& A, MatrB& B, VecS& s, Work workspace)
+ {
+ typedef typename MatrA::value_type val_t;
+
+ return detail::Gelsd<n_workspace_args<val_t>::value>() (A, B, s, workspace);
+ }
+
+ // gelsd, no singular values are returned
+ // Parameters:
+ // A: matrix of coefficients
+ // B: matrix of solutions (stored column-wise)
+ // workspace: either optimal, minimal, or user supplied
+ //
+ template <typename MatrA, typename MatrB, typename Work>
+ inline int gelsd(MatrA& A, MatrB& B, Work workspace)
+ {
+ typedef typename MatrA::value_type val_t;
+ typedef typename traits::type_traits<val_t>::real_type real_t;
+
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+
+ const int s_size = std::max(1, std::min(m,n));
+ traits::detail::array<real_t> s(s_size);
+
+ return detail::Gelsd<n_workspace_args<val_t>::value>() (A, B, s, workspace);
+ }
+
+ } // lapack
+
+}}}
+
+#endif
Added: sandbox/boost/numeric/bindings/lapack/gelss.hpp
==============================================================================
--- (empty file)
+++ sandbox/boost/numeric/bindings/lapack/gelss.hpp 2008-12-07 08:19:55 EST (Sun, 07 Dec 2008)
@@ -0,0 +1,354 @@
+//
+// Copyright Jesse Manning 2007
+//
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+//
+
+#ifndef BOOST_NUMERIC_BINDINGS_LAPACK_GELSS_HPP
+#define BOOST_NUMERIC_BINDINGS_LAPACK_GELSS_HPP
+
+#include <algorithm>
+
+#include <boost/numeric/bindings/traits/type.hpp>
+#include <boost/numeric/bindings/traits/traits.hpp>
+#include <boost/numeric/bindings/traits/type_traits.hpp>
+#include <boost/numeric/bindings/lapack/lapack.h>
+#include <boost/numeric/bindings/lapack/workspace.hpp>
+#include <boost/numeric/bindings/traits/detail/array.hpp>
+#include <boost/numeric/bindings/traits/detail/utils.hpp>
+
+
+#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
+# include <boost/static_assert.hpp>
+# include <boost/type_traits.hpp>
+#endif
+
+namespace boost { namespace numeric { namespace bindings {
+
+ namespace lapack {
+
+ namespace detail {
+
+ inline void gelss(const int m, const int n, const int nrhs,
+ float *a, const int lda, float *b, const int ldb,
+ float *s, const float rcond, int *rank, float *work,
+ const int lwork, int *info)
+ {
+ LAPACK_SGELSS(&m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, rank, work, &lwork, info);
+ }
+
+ inline void gelss(const int m, const int n, const int nrhs,
+ double *a, const int lda, double *b, const int ldb,
+ double *s, const double rcond, int *rank, double *work,
+ const int lwork, int *info)
+ {
+ LAPACK_DGELSS(&m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, rank, work, &lwork, info);
+ }
+
+ inline void gelss(const int m, const int n, const int nrhs,
+ traits::complex_f *a, const int lda, traits::complex_f *b,
+ const int ldb, float *s, const float rcond, int *rank,
+ traits::complex_f *work, const int lwork, float *rwork, int *info)
+ {
+ LAPACK_CGELSS(&m, &n, &nrhs, traits::complex_ptr(a),
+ &lda, traits::complex_ptr(b), &ldb, s,
+ &rcond, rank, traits::complex_ptr(work),
+ &lwork, rwork, info);
+ }
+
+ inline void gelss(const int m, const int n, const int nrhs,
+ traits::complex_d *a, const int lda, traits::complex_d *b,
+ const int ldb, double *s, const double rcond, int *rank,
+ traits::complex_d *work, const int lwork, double *rwork, int *info)
+ {
+ LAPACK_ZGELSS(&m, &n, &nrhs, traits::complex_ptr(a),
+ &lda, traits::complex_ptr(b), &ldb, s,
+ &rcond, rank, traits::complex_ptr(work),
+ &lwork, rwork, info);
+ }
+
+ // gelss for real type
+ template <typename MatrA, typename MatrB, typename VecS, typename Work>
+ inline int gelss(MatrA& A, MatrB& B, VecS& s, Work& work)
+ {
+ typedef typename MatrA::value_type val_t;
+ typedef typename traits::type_traits<val_t>::real_type real_t;
+
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int nrhs = traits::matrix_size2(B);
+ const int maxmn = std::max(m, n);
+ const int minmn = std::min(m, n);
+
+ // sanity checks
+ assert(m >= 0 && n >= 0);
+ assert(nrhs >= 0);
+ assert(traits::leading_dimension(A) >= std::max(1, m));
+ assert(traits::leading_dimension(B) >= std::max(1, maxmn));
+ assert(traits::vector_size(work) >= 1);
+ assert(traits::vector_size(s) >= std::max(1, minmn));
+
+ int info;
+ const real_t rcond = -1; // use machine precision
+ int rank;
+
+ detail::gelss(traits::matrix_size1(A),
+ traits::matrix_size2(A),
+ traits::matrix_size2(B),
+ traits::matrix_storage(A),
+ traits::leading_dimension(A),
+ traits::matrix_storage(B),
+ traits::leading_dimension(B),
+ traits::vector_storage(s),
+ rcond,
+ &rank,
+ traits::vector_storage(work),
+ traits::vector_size(work),
+ &info);
+
+ return info;
+ }
+
+ // gelss for complex type
+ template <typename MatrA, typename MatrB, typename VecS, typename Work, typename RWork>
+ inline int gelss(MatrA& A, MatrB& B, VecS& s, Work& work, RWork& rwork)
+ {
+ typedef typename MatrA::value_type val_t;
+ typedef typename traits::type_traits<val_t>::real_type real_t;
+
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int nrhs = traits::matrix_size2(B);
+ const int maxmn = std::max(m, n);
+ const int minmn = std::min(m, n);
+
+ // sanity checks
+ assert(m >= 0 && n >= 0);
+ assert(nrhs >= 0);
+ assert(traits::leading_dimension(A) >= std::max(1, m));
+ assert(traits::leading_dimension(B) >= std::max(1, maxmn));
+ assert(traits::vector_size(work) >= 1);
+ assert(traits::vector_size(s) >= std::max(1, minmn));
+
+ int info;
+ const real_t rcond = -1; // use machine precision
+ int rank;
+
+ detail::gelss(traits::matrix_size1(A),
+ traits::matrix_size2(A),
+ traits::matrix_size2(B),
+ traits::matrix_storage(A),
+ traits::leading_dimension(A),
+ traits::matrix_storage(B),
+ traits::leading_dimension(B),
+ traits::vector_storage(s),
+ rcond,
+ &rank,
+ traits::vector_storage(work),
+ traits::vector_size(work),
+ traits::vector_storage(rwork),
+ &info);
+
+ return info;
+ }
+
+ // default minimal workspace functor
+ template <int N>
+ struct Gelss { };
+
+ // specialization for gelss (sgelss, dgelss)
+ template <>
+ struct Gelss<1>
+ {
+ template <typename MatrA, typename MatrB, typename VecS>
+ inline int operator() (MatrA& A, MatrB& B, VecS& s, minimal_workspace) const
+ {
+ typedef typename MatrA::value_type val_t;
+
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int rhs = traits::matrix_size2(B);
+
+ const int minmn = std::min(m, n); // minmn = m < n ? m : n
+ const int maxmn = std::max(m, n); // maxmn = m > n ? m : n
+ const int maxmnr = std::max(maxmn, rhs); // maxmnr = maxmn > rhs ? maxmn : rhs
+
+ traits::detail::array<val_t> work(3*minmn + std::max(2*minmn, maxmnr));
+
+ return gelss(A, B, s, work);
+ }
+
+ template <typename MatrA, typename MatrB, typename VecS>
+ inline int operator() (MatrA& A, MatrB& B, VecS& s, optimal_workspace) const
+ {
+ typedef typename MatrA::value_type val_t;
+ typedef typename traits::type_traits<val_t>::real_type real_t;
+
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int rhs = traits::matrix_size2(B);
+
+ const int minmn = std::min(m, n); // minmn = m < n ? m : n
+ const int maxmn = std::max(m, n); // maxmn = m > n ? m : n
+ const int maxmnr = std::max(maxmn, rhs); // maxmnr = maxmn > rhs ? maxmn : rhs
+
+ val_t temp_work;
+
+ const real_t rcond = -1;
+ int rank;
+ int info;
+
+ // query for optimal workspace size
+ detail::gelss(traits::matrix_size1(A),
+ traits::matrix_size2(A),
+ traits::matrix_size2(B),
+ traits::matrix_storage(A),
+ traits::leading_dimension(A),
+ traits::matrix_storage(B),
+ traits::leading_dimension(B),
+ traits::vector_storage(s),
+ rcond,
+ &rank,
+ &temp_work, //traits::vector_storage(work),
+ -1, //traits::vector_size(work),
+ &info);
+
+ assert(info == 0);
+
+ const int lwork = traits::detail::to_int(temp_work);
+
+ traits::detail::array<val_t> work(lwork);
+
+ return gelss(A, B, s, work);
+ }
+
+ template <typename MatrA, typename MatrB, typename VecS, typename Work>
+ inline int operator() (MatrA& A, MatrB& B, VecS& s, detail::workspace1<Work> workspace) const
+ {
+ return gelss(A, B, s, workspace.w_);
+ }
+ };
+
+ // specialization for gelss (cgelss, zgelss)
+ template <>
+ struct Gelss<2>
+ {
+ template <typename MatrA, typename MatrB, typename VecS>
+ inline int operator() (MatrA& A, MatrB& B, VecS& s, minimal_workspace) const
+ {
+ typedef typename MatrA::value_type val_t;
+ typedef typename traits::type_traits<val_t>::real_type real_t;
+
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int rhs = traits::matrix_size2(B);
+
+ const int minmn = std::min(m, n); // minmn = m < n ? m : n
+ const int maxmn = std::max(m, n); // maxmn = m > n ? m : n
+ const int maxmnr = std::max(maxmn, rhs); // maxmnr = maxmn > rhs ? maxmn : rhs
+
+ traits::detail::array<val_t> work(2*minmn + maxmnr);
+ traits::detail::array<real_t> rwork(std::max(1, (5*minmn)));
+
+ return gelss(A, B, s, work, rwork);
+ }
+
+ template <typename MatrA, typename MatrB, typename VecS>
+ inline int operator() (MatrA& A, MatrB& B, VecS& s, optimal_workspace) const
+ {
+ typedef typename MatrA::value_type val_t;
+ typedef typename traits::type_traits<val_t>::real_type real_t;
+
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+ const int rhs = traits::matrix_size2(B);
+
+ const int minmn = std::min(m, n); // minmn = m < n ? m : n
+ const int maxmn = std::max(m, n); // maxmn = m > n ? m : n
+ const int maxmnr = std::max(maxmn, rhs); // maxmnr = maxmn > rhs ? maxmn : rhs
+
+ val_t temp_work;
+ real_t temp_rwork;
+
+ const real_t rcond = -1;
+ int rank;
+ int info;
+
+ // query for optimal workspace size
+ detail::gelss(traits::matrix_size1(A),
+ traits::matrix_size2(A),
+ traits::matrix_size2(B),
+ traits::matrix_storage(A),
+ traits::leading_dimension(A),
+ traits::matrix_storage(B),
+ traits::leading_dimension(B),
+ traits::vector_storage(s),
+ rcond,
+ &rank,
+ &temp_work, //traits::vector_storage(work),
+ -1, //traits::vector_size(work),
+ &temp_rwork,
+ &info);
+
+ assert(info == 0);
+
+ const int lwork = traits::detail::to_int(temp_work);
+
+ traits::detail::array<val_t> work(lwork);
+ traits::detail::array<real_t> rwork(std::max(1, (5*minmn)));
+
+ return gelss(A, B, s, work, rwork);
+ }
+
+ template <typename MatrA, typename MatrB, typename VecS, typename Work, typename RWork>
+ inline int operator() (MatrA& A, MatrB& B, VecS& s, detail::workspace2<Work, RWork> workspace) const
+ {
+ return gelss(A, B, s, workspace.w_, workspace.wr);
+ }
+ };
+
+ } // detail
+
+ // gelss
+ // Parameters:
+ // A: matrix of coefficients
+ // B: matrix of solutions (stored column-wise)
+ // s: vector to store singular values on output, length >= max(1, min(m,n))
+ // workspace: either optimal, minimal, or user supplied
+ //
+ template <typename MatrA, typename MatrB, typename VecS, typename Work>
+ inline int gelss(MatrA& A, MatrB& B, VecS& s, Work workspace)
+ {
+ typedef typename MatrA::value_type val_t;
+
+ return detail::Gelss<n_workspace_args<val_t>::value>() (A, B, s, workspace);
+ }
+
+ // gelss, no singular values are returned
+ // Parameters:
+ // A: matrix of coefficients
+ // B: matrix of solutions (stored column-wise)
+ // workspace: either optimal, minimal, or user supplied
+ //
+ template <typename MatrA, typename MatrB, typename Work>
+ inline int gelss(MatrA& A, MatrB& B, Work workspace)
+ {
+ typedef typename MatrA::value_type val_t;
+ typedef typename traits::type_traits<val_t>::real_type real_t;
+
+ const int m = traits::matrix_size1(A);
+ const int n = traits::matrix_size2(A);
+
+ const int s_size = std::max(1, std::min(m,n));
+ traits::detail::array<real_t> s(s_size);
+
+ return detail::Gelss<n_workspace_args<val_t>::value>() (A, B, s, workspace);
+ }
+
+ } // namespace lapack
+
+}}}
+
+#endif