[ublas] trsm ATLAS numerical bindings

Subject: [ublas] trsm ATLAS numerical bindings
From: Jörn Ungermann (j.ungermann_at_[hidden])
Date: 2009-04-15 07:16:02

Hi all!

We use the the trsm function as part of the calculation of K^T*(S^-1)*K,
whereby we have a Cholesky decomposition of S.
I found to my astonishment, that this is not included in the ATLAS
numerical bindings of uBLAS (and that also the existing uBLAS tsm is not
actually simple to use starting from a valid cBLAS call...)

Is there something wrong with the following code, which adds ATLAS
bindings for trsm and how do I get this into the standard release of the
bindings?
(I have to admit that I have no deep understanding of all these
trait-thingies...)

Kind regards,
Jörn

    namespace detail {

      inline
      void trsm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE
Side,
                       const enum CBLAS_UPLO Uplo, const enum
CBLAS_TRANSPOSE TransA,
                       const enum CBLAS_DIAG Diag, const int M, const
int N,
                       const float alpha, const float *A, const int lda,
                       float *B, const int ldb)
      {
        cblas_strsm (Order, Side, Uplo, TransA, Diag, M, N,
                     alpha, A, lda, B, ldb);
      }

      inline
      void trsm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE
Side,
                const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE
TransA,
                const enum CBLAS_DIAG Diag, const int M, const int N,
                const double alpha, const double *A, const int lda,
                double *B, const int ldb)
      {
        cblas_dtrsm (Order, Side, Uplo, TransA, Diag, M, N,
                     alpha, A, lda, B, ldb);
      }

      inline
      void trsm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE
Side,
                const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE
TransA,
                const enum CBLAS_DIAG Diag, const int M, const int N,
                traits::complex_f const& alpha, traits::complex_f const*
A, const int lda,
                traits::complex_f* B, const int ldb)
      {
        cblas_ctrsm (Order, Side, Uplo, TransA, Diag, M, N,
                     static_cast<void const*> (&alpha),
                     static_cast<void const*> (A), lda,
                     static_cast<void*> (B), ldb);
      }

      inline
      void trsm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE
Side,
                const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE
TransA,
                const enum CBLAS_DIAG Diag, const int M, const int N,
                traits::complex_d const& alpha, traits::complex_d const*
A, const int lda,
                traits::complex_d* B, const int ldb)
      {
        cblas_ztrsm (Order, Side, Uplo, TransA, Diag, M, N,
                     static_cast<void const*> (&alpha),
                     static_cast<void const*> (A), lda,
                     static_cast<void*> (B), ldb);
      }
    }

    template<typename T, typename MatrA, typename MatrB>
    inline
    void trsm(CBLAS_SIDE const Side, CBLAS_UPLO const Uplo,
              CBLAS_TRANSPOSE const TransA, CBLAS_DIAG const Diag,
              double const& alpha, DenseMatrix const& a, DenseMatrix& b)
    {
      // Mostly copied and adapted from ATLAS numeric bindings for dsymm
function.
      int const ma = int(traits::matrix_size1 (a));
      int const na = int(traits::matrix_size2 (a));
      int const mb = int(traits::matrix_size1 (b));
      int const nb = int(traits::matrix_size2 (b));

      assert(ma == na);

      CBLAS_ORDER const stor_ord
        = enum_cast<CBLAS_ORDER const>
        (storage_order<
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
           typename traits::matrix_traits<MatrA>::ordering_type
#else
           typename MatrA::orientation_category
#endif
         >::value);

      assert((Side == CblasLeft && mb == ma) || (Side == CblasRight &&
nb == ma));
      detail::trsm (stor_ord, Side, Uplo, TransA, Diag, mb, nb, alpha,
#ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
                   traits::matrix_storage(a),
#else
                   traits::matrix_storage_const(a),
#endif
                   int(traits::leading_dimension(a)),
#ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
                   traits::matrix_storage(b),
#else
                   traits::matrix_storage_const(b),
#endif
                   int(traits::leading_dimension(b)));
    }

-- 
Jörn Ungermann
Dipl.-Mathematiker
ICG-1, Forschungszentrum Jülich
Tel.: [phone]