[Boost-commit] svn:boost r50176 - sandbox/boost/numeric/bindings/lapack

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