Skip to content
This repository has been archived by the owner on Oct 19, 2020. It is now read-only.

Commit

Permalink
Add obvious improvements from current bindings svn.
Browse files Browse the repository at this point in the history
  • Loading branch information
@inducer
inducer committed Nov 12, 2008
1 parent adf39e6 commit 3fd47cf
Show file tree
Hide file tree
Showing 54 changed files with 2,258 additions and 254 deletions.
4 changes: 4 additions & 0 deletions boost/numeric/bindings/lapack/ormqr.hpp
View file
Expand Up @@ -72,6 +72,7 @@ namespace boost { namespace numeric { namespace bindings {
int const ldc, float* work, int const lwork,
int& info)
{
assert ( trans=='N' || trans=='T' );
LAPACK_SORMQR (&side, &trans, &m, &n, &k,
a, &lda,
tau,
Expand All @@ -87,6 +88,7 @@ namespace boost { namespace numeric { namespace bindings {
int const ldc, double* work, int const lwork,
int& info)
{
assert ( trans=='N' || trans=='T' );
LAPACK_DORMQR (&side, &trans, &m, &n, &k,
a, &lda,
tau,
Expand All @@ -102,6 +104,7 @@ namespace boost { namespace numeric { namespace bindings {
int const ldc, traits::complex_f* work, int const lwork,
int& info)
{
assert ( trans=='N' || trans=='C' );
LAPACK_CUNMQR (&side, &trans, &m, &n, &k,
traits::complex_ptr(a), &lda,
traits::complex_ptr(tau),
Expand All @@ -117,6 +120,7 @@ namespace boost { namespace numeric { namespace bindings {
int const ldc, traits::complex_d* work, int const lwork,
int& info)
{
assert ( trans=='N' || trans=='C' );
LAPACK_ZUNMQR (&side, &trans, &m, &n, &k,
traits::complex_ptr(a), &lda,
traits::complex_ptr(tau),
Expand Down
209 changes: 209 additions & 0 deletions boost/numeric/bindings/lapack/ptsv.hpp
View file
@@ -0,0 +1,209 @@
/*
*
* Copyright (c) Karl Meerbergen 2008
*
* 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_PTSV_HPP
#define BOOST_NUMERIC_BINDINGS_LAPACK_PTSV_HPP

#include <boost/numeric/bindings/traits/type_traits.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>
#endif

#include <cassert>

namespace boost { namespace numeric { namespace bindings {

namespace lapack {

/////////////////////////////////////////////////////////////////////
//
// system of linear equations A * X = B
// with A tridiagonal symmetric or Hermitian positive definite matrix
//
/////////////////////////////////////////////////////////////////////

/*
* ptsv() computes the solution to a system of linear equations
* A*X = B, where A is an N-by-N Hermitian positive definite tridiagonal
* matrix, and X and B are N-by-NRHS matrices.
*
* A is factored as A = L*D*L**T, and the factored form of A is then
* used to solve the system of equations.
*/

namespace detail {

inline
void ptsv ( int const n, int const nrhs,
float* d, float* e, float* b, int const ldb, int* info)
{
LAPACK_SPTSV (&n, &nrhs, d, e, b, &ldb, info);
}

inline
void ptsv ( int const n, int const nrhs,
double* d, double* e, double* b, int const ldb, int* info)
{
LAPACK_DPTSV (&n, &nrhs, d, e, b, &ldb, info);
}

inline
void ptsv ( int const n, int const nrhs,
float* d, traits::complex_f* e, traits::complex_f* b, int const ldb,
int* info)
{
LAPACK_CPTSV (&n, &nrhs, d, traits::complex_ptr(e), traits::complex_ptr(b), &ldb, info);
}

inline
void ptsv ( int const n, int const nrhs,
double* d, traits::complex_d* e, traits::complex_d* b, int const ldb,
int* info)
{
LAPACK_ZPTSV (&n, &nrhs, d, traits::complex_ptr(e), traits::complex_ptr(b), &ldb, info);
}

}

template <typename D, typename E, typename B>
inline int ptsv( D& d, E& e, B& b ) {
int const n = traits::vector_size(d) ;
assert( n==traits::vector_size(e)+1 ) ;
assert( n==traits::matrix_num_rows(b) ) ;

int info ;
detail::ptsv( n, traits::matrix_num_columns (b)
, traits::vector_storage(d)
, traits::vector_storage(e)
, traits::matrix_storage(b)
, traits::leading_dimension(b)
, &info
) ;
return info ;
} // ptsv()


/*
* pttrf() computes the L * D * L^H factorization of a Hermitian
* positive definite tridiagonal matrix A. The factorization may also
* be regarded as having the form A = U^H * D *U.
*/

namespace detail {

inline
void pttrf ( int const n, float* d, float* e, int* info) {
LAPACK_SPTTRF ( &n, d, e, info) ;
}

inline
void pttrf ( int const n, double* d, double* e, int* info) {
LAPACK_DPTTRF ( &n, d, e, info);
}

inline
void pttrf ( int const n, float* d, traits::complex_f* e, int* info)
{
LAPACK_CPTTRF ( &n, d, traits::complex_ptr(e), info);
}

inline
void pttrf ( int const n, double* d, traits::complex_d* e, int* info)
{
LAPACK_ZPTTRF ( &n, d, traits::complex_ptr(e), info);
}

}

template <typename D, typename E>
inline
int pttrf (D& d, E& e) {
int const n = traits::vector_size (d);
assert (n == traits::vector_size (e) + 1);
int info;
detail::pttrf ( n, traits::vector_storage(d), traits::vector_storage(e), &info);
return info;
}


/*
* pttrs() solves a tridiagonal system of the form
* A * X = B
* using the factorization A = U^H * D * U or A = L * D * L^H computed by pttrf().
* D is a diagonal matrix specified in the vector D, U (or L) is a unit
* bidiagonal matrix whose superdiagonal (subdiagonal) is specified in
* the vector E, and X and B are N by NRHS matrices.
*/

namespace detail {

inline
void pttrs (char const uplo, int const n, int const nrhs,
float const* d, float const* e, float* b, int const ldb, int* info)
{
LAPACK_SPTTRS (&n, &nrhs, d, e, b, &ldb, info);
}

inline
void pttrs (char const uplo, int const n, int const nrhs,
double const* d, double const* e, double* b, int const ldb, int* info)
{
LAPACK_DPTTRS (&n, &nrhs, d, e, b, &ldb, info);
}

inline
void pttrs (char const uplo, int const n, int const nrhs,
float const* d,
traits::complex_f const* e,
traits::complex_f* b, int const ldb, int* info)
{
LAPACK_CPTTRS (&uplo, &n, &nrhs, d,
traits::complex_ptr (e),
traits::complex_ptr (b), &ldb, info);
}

inline
void pttrs (char const uplo, int const n, int const nrhs,
double const* d,
traits::complex_d const* e,
traits::complex_d* b, int const ldb, int* info)
{
LAPACK_ZPTTRS (&uplo, &n, &nrhs, d,
traits::complex_ptr (e),
traits::complex_ptr (b), &ldb, info);
}

}

template <typename D, typename E, typename MatrB>
inline
int pttrs (char uplo, D const& d, E const& e, MatrB& b) {
int const n = traits::vector_size (d);
assert (n == traits::vector_size (e) + 1);
assert (n == traits::matrix_num_rows (b));

int info;
detail::pttrs (uplo, n, traits::matrix_num_columns (b),
traits::vector_storage (d),
traits::vector_storage (e),
traits::matrix_storage (b),
traits::leading_dimension (b),
&info);
return info;
} // pttrs()

}

}}}

#endif
76 changes: 76 additions & 0 deletions boost/numeric/bindings/lapack/syev.hpp
View file
Expand Up @@ -148,8 +148,84 @@ namespace boost { namespace numeric { namespace bindings {
return syev(jobz, uplo, a, w, optimal_workspace());
} // syev()

//
// With UPLO integrated in matrix type
//
template <typename A, typename W>
inline
int syev (char jobz, A& a, W& w, optimal_workspace ) {
typedef typename A::value_type value_type ;

int const n = traits::matrix_size1 (a);
char uplo = traits::matrix_uplo_tag( a ) ;
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
typedef typename traits::matrix_traits<A>::matrix_structure matrix_structure ;
BOOST_STATIC_ASSERT( (boost::mpl::or_< boost::is_same< matrix_structure, traits::symmetric_t >
, boost::is_same< matrix_structure, traits::hermitian_t >
>::value)
) ;
#endif

traits::detail::array<value_type> work( std::max<int>(1,34*n) );
return detail::syev(jobz, uplo, a, w, work);
} // syev()


// Function that allocates work arrays
template <typename A, typename W>
inline
int syev (char jobz, A& a, W& w, minimal_workspace ) {
typedef typename A::value_type value_type ;

int const n = traits::matrix_size1 (a);
char uplo = traits::matrix_uplo_tag( a ) ;
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
typedef typename traits::matrix_traits<A>::matrix_structure matrix_structure ;
BOOST_STATIC_ASSERT( (boost::mpl::or_< boost::is_same< matrix_structure, traits::symmetric_t >
, boost::is_same< matrix_structure, traits::hermitian_t >
>::value)
) ;
#endif
traits::detail::array<value_type> work( std::max<int>(1,3*n-1) );
return detail::syev(jobz, uplo, a, w, work);
} // syev()


// Function that allocates work arrays
template <typename A, typename W, typename Work>
inline
int syev (char jobz, A& a, W& w, detail::workspace1<Work> workspace ) {
typedef typename A::value_type value_type ;
char uplo = traits::matrix_uplo_tag( a ) ;
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
typedef typename traits::matrix_traits<A>::matrix_structure matrix_structure ;
BOOST_STATIC_ASSERT( (boost::mpl::or_< boost::is_same< matrix_structure, traits::symmetric_t >
, boost::is_same< matrix_structure, traits::hermitian_t >
>::value)
) ;
#endif
return detail::syev(jobz, uplo, a, w, workspace.w_);
} // syev()

// Function without workarray as argument
template <typename A, typename W>
inline
int syev (char jobz, A& a, W& w) {
char uplo = traits::matrix_uplo_tag( a ) ;
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
typedef typename traits::matrix_traits<A>::matrix_structure matrix_structure ;
BOOST_STATIC_ASSERT( (boost::mpl::or_< boost::is_same< matrix_structure, traits::symmetric_t >
, boost::is_same< matrix_structure, traits::hermitian_t >
>::value)
) ;
#endif
return syev(jobz, uplo, a, w, optimal_workspace());
} // syev()

}



}}}

#endif
2 changes: 1 addition & 1 deletion boost/numeric/bindings/mumps/mumps_driver.hpp
View file
Expand Up @@ -9,6 +9,6 @@
#ifndef BOOST_NUMERIC_BINDINGS_MUMPS_MUMPS_DRIVER_HPP
#define BOOST_NUMERIC_BINDINGS_MUMPS_MUMPS_DRIVER_HPP

#include <boost/numeric/bindings/mumps/mumps_driver_4_6_4.hpp>
#include <boost/numeric/bindings/mumps/mumps_driver_4_8_0.hpp>

#endif
6 changes: 3 additions & 3 deletions boost/numeric/bindings/mumps/mumps_driver_4_6_4.hpp
View file
Expand Up @@ -210,8 +210,8 @@ namespace boost { namespace numeric { namespace bindings { namespace mumps {
template <typename M>
void matrix_integer_data( mumps<M>& data, M& m ) {
BOOST_STATIC_ASSERT( (1 == boost::numeric::bindings::traits::sparse_matrix_traits<M>::index_base) ) ;
data.n = boost::numeric::bindings::traits::spmatrix_size1( m ) ;
assert( boost::numeric::bindings::traits::spmatrix_size2( m ) == data.n ) ;
data.n = boost::numeric::bindings::traits::spmatrix_num_rows( m ) ;
assert( boost::numeric::bindings::traits::spmatrix_num_columns( m ) == data.n ) ;

data.nz = boost::numeric::bindings::traits::spmatrix_num_nonzeros( m ) ;
detail::indices( typename boost::numeric::bindings::traits::sparse_matrix_traits<M>::ordering_type(), data.irn, data.jcn, m ) ;
Expand All @@ -238,7 +238,7 @@ namespace boost { namespace numeric { namespace bindings { namespace mumps {
template <typename M, typename X>
void rhs_sol_value_data( mumps<M>& data, X& x ) {
data.rhs = detail::cast_2_mumps( boost::numeric::bindings::traits::matrix_storage( x ) ) ;
data.nrhs = boost::numeric::bindings::traits::matrix_size2( x ) ;
data.nrhs = boost::numeric::bindings::traits::matrix_num_columns( x ) ;
data.lrhs = boost::numeric::bindings::traits::leading_dimension( x ) ;
} // matrix_rhs_sol_value_data()

Expand Down

0 comments on commit 3fd47cf

Please sign in to comment.