From: karl.meerbergen_at_[hidden]
Date: 2007-08-28 03:06:59
Author: karlmeerbergen
Date: 2007-08-28 03:06:57 EDT (Tue, 28 Aug 2007)
New Revision: 39030
URL: http://svn.boost.org/trac/boost/changeset/39030
Log:
added bindings to LAPACK's *steqr and *sytrd
Added:
sandbox/boost/numeric/bindings/lapack/steqr.hpp (contents, props changed)
sandbox/boost/numeric/bindings/lapack/sytrd.hpp (contents, props changed)
sandbox/libs/numeric/bindings/lapack/test/ublas_steqr.cpp (contents, props changed)
sandbox/libs/numeric/bindings/lapack/test/ublas_sytrd.cpp (contents, props changed)
Text files modified:
sandbox/boost/numeric/bindings/lapack/lapack.h | 23 +++++++++++++++++++----
sandbox/boost/numeric/bindings/lapack/lapack_names.h | 9 +++++++++
2 files changed, 28 insertions(+), 4 deletions(-)
Modified: sandbox/boost/numeric/bindings/lapack/lapack.h
==============================================================================
--- sandbox/boost/numeric/bindings/lapack/lapack.h (original)
+++ sandbox/boost/numeric/bindings/lapack/lapack.h 2007-08-28 03:06:57 EDT (Tue, 28 Aug 2007)
@@ -365,26 +365,41 @@
int* ifst, const int * ilst, int* info );
+ /* Hermitian tridiagonal matrices */
+
+ void LAPACK_SSTEQR( char const* compz, int const* n, float* d, float* E, float* z, int const* ldz, float* work, int* info ) ;
+ void LAPACK_DSTEQR( char const* compz, int const* n, double* d, double* E, double* z, int const* ldz, double* work, int* info ) ;
+
/* Hermitian banded matrices */
void LAPACK_SSBEV( char const* jobz, char const* uplo, int const* n,
int const* kd, float* ab, int const* ldab, float* w,
- float* z, int const* ldz, float* work, int const* info );
+ float* z, int const* ldz, float* work, int* info );
void LAPACK_DSBEV( char const* jobz, char const* uplo, int const* n,
int const* kd, double* ab, int const* ldab, double* w,
- double* z, int const* ldz, double* work, int const* info );
+ double* z, int const* ldz, double* work, int* info );
void LAPACK_CHBEV( char const* jobz, char const* uplo, int const* n,
int const* kd, fcomplex_t* ab, int const* ldab, float* w,
fcomplex_t* z, int const* ldz, fcomplex_t* work,
- float* rwork, int const* info );
+ float* rwork, int* info );
void LAPACK_ZHBEV( char const* jobz, char const* uplo, int const* n,
int const* kd, dcomplex_t* ab, int const* ldab, double* w,
dcomplex_t* z, int const* ldz, dcomplex_t* work,
- double* rwork, int const* info );
+ double* rwork, int* info );
+
+
+ /*********************************************************************/
+ /* Auxiliary routines for eigenvalue problems */
+ /*********************************************************************/
+
+ void LAPACK_SSYTRD( char const* uplo, int const* n, float* a, int const* lda, float* d,
+ float* e, float* tau, float* work, int const* lwork, int* INFO ) ;
+ void LAPACK_DSYTRD( char const* uplo, int const* n, double* a, int const* lda, double* d,
+ double* e, double* tau, double* work, int const* lwork, int* INFO ) ;
/*********************************************************************/
/* SVD */
Modified: sandbox/boost/numeric/bindings/lapack/lapack_names.h
==============================================================================
--- sandbox/boost/numeric/bindings/lapack/lapack_names.h (original)
+++ sandbox/boost/numeric/bindings/lapack/lapack_names.h 2007-08-28 03:06:57 EDT (Tue, 28 Aug 2007)
@@ -181,6 +181,12 @@
#define LAPACK_ZHBEV FORTRAN_ID( zhbev )
+/********************************************/
+/* eigenproblems for tridiagonal matrices */
+
+#define LAPACK_SSTEQR FORTRAN_ID( ssteqr )
+#define LAPACK_DSTEQR FORTRAN_ID( dsteqr )
+
/********************************************/
/* QR factorization */
@@ -197,6 +203,9 @@
#define LAPACK_CUNMQR FORTRAN_ID( cunmqr )
#define LAPACK_ZUNMQR FORTRAN_ID( zunmqr )
+#define LAPACK_SSYTRD FORTRAN_ID( ssytrd )
+#define LAPACK_DSYTRD FORTRAN_ID( dsytrd )
+
/********************************************/
/* SVD */
Added: sandbox/boost/numeric/bindings/lapack/steqr.hpp
==============================================================================
--- (empty file)
+++ sandbox/boost/numeric/bindings/lapack/steqr.hpp 2007-08-28 03:06:57 EDT (Tue, 28 Aug 2007)
@@ -0,0 +1,108 @@
+//
+// Copyright Karl Meerbergen 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_STEQR_HPP
+#define BOOST_NUMERIC_BINDINGS_LAPACK_STEQR_HPP
+
+#include <boost/numeric/bindings/lapack/workspace.hpp>
+#include <boost/numeric/bindings/traits/vector_traits.hpp>
+#include <boost/numeric/bindings/traits/matrix_traits.hpp>
+#include <boost/numeric/bindings/traits/type_traits.hpp>
+#include <boost/numeric/bindings/traits/detail/array.hpp>
+#include <boost/numeric/bindings/traits/traits.hpp>
+#include <boost/numeric/bindings/lapack/lapack.h>
+
+#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
+# include <boost/static_assert.hpp>
+# include <boost/type_traits/is_same.hpp>
+#endif
+
+
+ /********************************************************************/
+ /* eigenvalue problems */
+ /********************************************************************/
+
+ /* tridiagonal symmetric */
+
+
+namespace boost { namespace numeric { namespace bindings { namespace lapack {
+
+ namespace detail {
+
+ inline
+ void steqr ( char compz, int n, float* d, float* e, float* z, int ldz, float* work, int& info )
+ {
+ LAPACK_SSTEQR( &compz, &n, d, e, z, &ldz, work, &info ) ;
+ }
+
+ inline
+ void steqr ( char compz, int n, double* d, double* e, double* z, int ldz, double* work, int& info )
+ {
+ LAPACK_DSTEQR( &compz, &n, d, e, z, &ldz, work, &info ) ;
+ }
+
+ } // namespace detail
+
+
+ template <typename D, typename E, typename Z, typename W>
+ inline
+ int steqr( char compz, D& d, E& e, Z& z, W& work ) {
+
+ int const n = traits::vector_size (d);
+ assert( traits::vector_size (e) == n-1 );
+ assert( traits::matrix_size1 (z) == n );
+ assert( traits::matrix_size2 (z) == n );
+ assert( compz=='N' || compz=='V' || compz=='I' );
+
+ int lwork = traits::vector_size( work ) ;
+
+ int info;
+ detail::steqr( compz, n,
+ traits::vector_storage( d ),
+ traits::vector_storage( e ),
+ traits::matrix_storage( z ),
+ traits::leading_dimension( z ),
+ traits::vector_storage( work ),
+ info ) ;
+ return info;
+ } // steqr()
+
+
+ template <typename D, typename E, typename Z>
+ inline
+ int steqr( char compz, D& d, E& e, Z& z, optimal_workspace ) {
+ int lwork = 0 ;
+ if (compz != 'N') lwork = 2 * traits::vector_size( d ) - 2 ;
+
+ traits::detail::array<typename traits::vector_traits<D>::value_type> work( lwork );
+
+ return steqr( compz, d, e, z, work ) ;
+ }
+
+
+ template <typename D, typename E, typename Z>
+ inline
+ int steqr( char compz, D& d, E& e, Z& z, minimal_workspace ) {
+ int lwork = 1 ;
+ if (compz != 'N') lwork = 2 * traits::vector_size( e ) ;
+
+ traits::detail::array<typename traits::vector_traits<D>::value_type> work( lwork );
+
+ return steqr( compz, d, e, z, work ) ;
+ }
+
+ template <typename D, typename E, typename Z>
+ inline
+ int steqr( char compz, D& d, E& e, Z& z ) {
+ return steqr( compz, d, e, z, minimal_workspace() ) ;
+ }
+
+
+}}}}
+
+#endif // BOOST_NUMERIC_BINDINGS_LAPACK_GBSV_HPP
Added: sandbox/boost/numeric/bindings/lapack/sytrd.hpp
==============================================================================
--- (empty file)
+++ sandbox/boost/numeric/bindings/lapack/sytrd.hpp 2007-08-28 03:06:57 EDT (Tue, 28 Aug 2007)
@@ -0,0 +1,111 @@
+//
+// Copyright Karl Meerbergen 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_SYTRD_HPP
+#define BOOST_NUMERIC_BINDINGS_LAPACK_SYTRD_HPP
+
+#include <boost/numeric/bindings/lapack/workspace.hpp>
+#include <boost/numeric/bindings/traits/vector_traits.hpp>
+#include <boost/numeric/bindings/traits/matrix_traits.hpp>
+#include <boost/numeric/bindings/traits/type_traits.hpp>
+#include <boost/numeric/bindings/traits/detail/array.hpp>
+#include <boost/numeric/bindings/traits/traits.hpp>
+#include <boost/numeric/bindings/lapack/lapack.h>
+
+#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
+# include <boost/static_assert.hpp>
+# include <boost/type_traits/is_same.hpp>
+#endif
+
+
+ /********************************************************************/
+ /* eigenvalue problems */
+ /********************************************************************/
+
+ /* tridiagonal symmetric */
+
+
+namespace boost { namespace numeric { namespace bindings { namespace lapack {
+
+ namespace detail {
+
+ inline
+ void sytrd ( char uplo, int n, float* a, int lda, float* d, float* e, float* tau, float* work, int lwork, int& info )
+ {
+ LAPACK_SSYTRD( &uplo, &n, a, &lda, d, e, tau, work, &lwork, &info ) ;
+ }
+
+ inline
+ void sytrd ( char uplo, int n, double* a, int lda, double* d, double* e, double* tau, double* work, int lwork, int& info )
+ {
+ LAPACK_DSYTRD( &uplo, &n, a, &lda, d, e, tau, work, &lwork, &info ) ;
+ }
+
+ } // namespace detail
+
+
+ template <typename A, typename D, typename E, typename Tau, typename W>
+ inline
+ int sytrd( char uplo, A& a, D& d, E& e, Tau& tau, W& work ) {
+
+ int const n = traits::matrix_size1 (a);
+ assert( traits::matrix_size2 (a) == n );
+ assert( traits::vector_size (d) == n );
+ assert( traits::vector_size (e) == n-1 );
+ assert( traits::vector_size (tau) == n-1 );
+ assert( uplo=='U' || uplo=='L' );
+
+ int lwork = traits::vector_size( work ) ;
+ assert( lwork >= 1 );
+
+ int info;
+ detail::sytrd( uplo, n,
+ traits::matrix_storage( a ),
+ traits::leading_dimension( a ),
+ traits::vector_storage( d ),
+ traits::vector_storage( e ),
+ traits::vector_storage( tau ),
+ traits::vector_storage( work ), lwork,
+ info ) ;
+ return info;
+ } // sytrd()
+
+
+ template <typename A, typename D, typename E, typename Tau>
+ inline
+ int sytrd( char uplo, A& a, D& d, E& e, Tau& tau, optimal_workspace=optimal_workspace() ) {
+ int info;
+ detail::sytrd( uplo, traits::matrix_size1( a ),
+ traits::matrix_storage( a ),
+ traits::leading_dimension( a ),
+ traits::vector_storage( d ),
+ traits::vector_storage( e ),
+ traits::vector_storage( tau ),
+ traits::vector_storage( tau ), -1,
+ info ) ;
+ if (info) return info ;
+ int lwork = * traits::vector_storage( tau ) ;
+
+ traits::detail::array<typename traits::vector_traits<D>::value_type> work( lwork );
+
+ return sytrd( uplo, a, d, e, tau, work ) ;
+ }
+
+
+ template <typename A, typename D, typename E, typename Tau>
+ inline
+ int sytrd( char uplo, A& a, D& d, E& e, Tau& tau, minimal_workspace ) {
+ int lwork = 1 ;
+ traits::detail::array<typename traits::vector_traits<D>::value_type> work( lwork );
+
+ return sytrd( uplo, a, d, e, tau, work ) ;
+ }
+
+}}}}
+
+#endif // BOOST_NUMERIC_BINDINGS_LAPACK_GBSV_HPP
Added: sandbox/libs/numeric/bindings/lapack/test/ublas_steqr.cpp
==============================================================================
--- (empty file)
+++ sandbox/libs/numeric/bindings/lapack/test/ublas_steqr.cpp 2007-08-28 03:06:57 EDT (Tue, 28 Aug 2007)
@@ -0,0 +1,72 @@
+// Permission to copy, use, modify, sell and
+// distribute this software is granted provided this copyright notice appears
+// in all copies. This software is provided "as is" without express or implied
+// warranty, and with no claim as to its suitability for any purpose.
+// Copyright Toon Knapen, Karl Meerbergen
+//
+
+#include "../../blas/test/random.hpp"
+
+#include <boost/numeric/bindings/lapack/steqr.hpp>
+#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
+#include <boost/numeric/bindings/traits/ublas_vector.hpp>
+#include <boost/numeric/ublas/io.hpp>
+#include <algorithm>
+#include <limits>
+#include <iostream>
+
+
+namespace ublas = boost::numeric::ublas;
+namespace lapack = boost::numeric::bindings::lapack;
+
+
+template <typename T>
+int do_value_type() {
+ const int n = 10 ;
+
+ typedef typename boost::numeric::bindings::traits::type_traits<T>::real_type real_type ;
+ typedef std::complex< real_type > complex_type ;
+
+ typedef ublas::matrix<T, ublas::column_major> matrix_type ;
+ typedef ublas::vector<T> vector_type ;
+
+ // Set matrix
+ matrix_type z( n, n );
+ vector_type d( n ), e ( n - 1 ) ;
+
+ std::fill( d.begin(), d.end(), 2.0 ) ;
+ std::fill( e.begin(), e.end(), -1.0 ) ;
+
+ // Compute eigendecomposition.
+ lapack::steqr( 'I', d, e, z ) ;
+
+ for ( int i=0; i<d.size(); ++i) {
+ T sum( 0.0 ) ;
+ for (int j=0; j<d.size(); ++j) {
+ sum += z(i,j)*z(i,j) * d(j) ;
+ }
+ if (std::abs( sum - 2.0 ) > 10 * std::numeric_limits<T>::epsilon() ) return 1 ;
+
+ if (i>0) {
+ sum = 0.0 ;
+ for (int j=0; j<d.size(); ++j) {
+ sum += z(i-1,j)*z(i,j) * d(j) ;
+ }
+ if (std::abs( sum + 1.0 ) > 10 * std::numeric_limits<T>::epsilon() ) return 1 ;
+ }
+ }
+
+ return 0 ;
+} // do_value_type()
+
+
+
+int main() {
+ // Run tests for different value_types
+ if (do_value_type<float>()) return 255;
+ if (do_value_type<double>()) return 255;
+
+ std::cout << "Regression test succeeded\n" ;
+ return 0;
+}
+
Added: sandbox/libs/numeric/bindings/lapack/test/ublas_sytrd.cpp
==============================================================================
--- (empty file)
+++ sandbox/libs/numeric/bindings/lapack/test/ublas_sytrd.cpp 2007-08-28 03:06:57 EDT (Tue, 28 Aug 2007)
@@ -0,0 +1,72 @@
+// Permission to copy, use, modify, sell and
+// distribute this software is granted provided this copyright notice appears
+// in all copies. This software is provided "as is" without express or implied
+// warranty, and with no claim as to its suitability for any purpose.
+// Copyright Toon Knapen, Karl Meerbergen
+//
+
+#include "../../blas/test/random.hpp"
+
+#include <boost/numeric/bindings/lapack/sytrd.hpp>
+#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
+#include <boost/numeric/bindings/traits/ublas_vector.hpp>
+#include <boost/numeric/ublas/io.hpp>
+#include <algorithm>
+#include <limits>
+#include <iostream>
+
+
+namespace ublas = boost::numeric::ublas;
+namespace lapack = boost::numeric::bindings::lapack;
+
+
+template <typename T>
+int do_value_type() {
+ const int n = 10 ;
+
+ typedef typename boost::numeric::bindings::traits::type_traits<T>::real_type real_type ;
+ typedef std::complex< real_type > complex_type ;
+
+ typedef ublas::matrix<T, ublas::column_major> matrix_type ;
+ typedef ublas::vector<T> vector_type ;
+
+ // Set matrix
+ matrix_type a( n, n );
+ vector_type d( n ), e ( n - 1 ), tau( n - 1 ) ;
+
+ for (int i=0; i<n; ++i ) {
+ for (int j=0; j<n; ++j ) {
+ a(j,i) = 0.0 ;
+ }
+ }
+
+ a(0,0) = 2.0 ;
+ for (int i=1; i<n; ++i ) {
+ a(i,i) = 2.0 ;
+ a(i-1,i) = -1.0 ;
+ }
+
+ // Compute eigendecomposition.
+ lapack::sytrd( 'U', a, d, e, tau ) ;
+
+ for ( int i=0; i<d.size(); ++i) {
+ if (std::abs( d(i) - 2.0 ) > 10 * std::numeric_limits<T>::epsilon() ) return 1 ;
+ }
+ for ( int i=0; i<e.size(); ++i) {
+ if (std::abs( e(i) + 1.0 ) > 10 * std::numeric_limits<T>::epsilon() ) return 1 ;
+ }
+
+ return 0 ;
+} // do_value_type()
+
+
+
+int main() {
+ // Run tests for different value_types
+ if (do_value_type<float>()) return 255;
+ if (do_value_type<double>()) return 255;
+
+ std::cout << "Regression test succeeded\n" ;
+ return 0;
+}
+