Ublas :

Date view	Thread view	Subject view	Author view

From: Karl Meerbergen (Karl.Meerbergen_at_[hidden])
Date: 2006-05-12 02:25:45

Next message: Thomas Lemaire: "Re: [ublas] minor bug in lapack bindings gesdd"
Previous message: Thomas Lemaire: "Re: [ublas] minor bug in lapack bindings gesdd"
In reply to: Thomas Lemaire: "Re: [ublas] minor bug in lapack bindings gesdd"
Next in thread: Thomas Lemaire: "Re: [ublas] minor bug in lapack bindings gesdd"
Reply: Thomas Lemaire: "Re: [ublas] minor bug in lapack bindings gesdd"

I had a look at it. I think I fixed it correctly. I am unable to commit,
because a connection problem with the server.
Attached is the corrected file. You can use it in the maintime.

Karl

Thomas Lemaire wrote:

>is anybody able to fix it ?
>
>On Monday 08 May 2006 19:49, Thomas Lemaire wrote:
>
>
>>dear list,
>>
>>when using gesdd I get the following error:
>>boost-sandbox/boost/numeric/bindings/lapack/gesdd.hpp:181: error: 'cerr' is
>>not a member of 'std'
>>
>>Indeed, gesdd uses std:cerr to report an error, grepping in the bindings
>>reveals that this is the only place where std:cerr is used. Maybe this
>>error should be changed to some assert() call.
>>In the meantime, #include <iostream>
>>
>>cheers,
>>
>>
>
>
>

/*
*
* Copyright (c) Toon Knapen & Kresimir Fresl 2003
*
* Permission to copy, modify, use and distribute this software
* for any non-commercial or commercial purpose is granted provided
* that this license appear on all copies of the software source code.
*
* Authors assume no responsibility whatsoever for its use and makes
* no guarantees about its quality, correctness or reliability.
*
* KF acknowledges the support of the Faculty of Civil Engineering,
* University of Zagreb, Croatia.
*
*/

#ifndef BOOST_NUMERIC_BINDINGS_LAPACK_GESDD_HPP
#define BOOST_NUMERIC_BINDINGS_LAPACK_GESDD_HPP

#include <boost/numeric/bindings/traits/type_traits.hpp>
#include <boost/numeric/bindings/traits/traits.hpp>
#include <boost/numeric/bindings/lapack/lapack.h>
#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/is_same.hpp>
#endif
#include <cassert>

namespace boost { namespace numeric { namespace bindings {

namespace lapack {

    ///////////////////////////////////////////////////////////////////
    //
    // singular value decomposition
    //
    ///////////////////////////////////////////////////////////////////

    /*
     * (divide and conquer driver)
     * gesdd() computes the singular value decomposition (SVD) of
     * M-by-N matrix A, optionally computing the left and/or right
     * singular vectors, by using divide-and-conquer method.
     * The SVD is written
     *
     * A = U * S * V^T or A = U * S * V^H
     *
     * where S is an M-by-N matrix which is zero except for its min(m,n)
     * diagonal elements, U is an M-by-M orthogonal/unitary matrix, and V
     * is an N-by-N orthogonal/unitary matrix. The diagonal elements of S
     * are the singular values of A; they are real and non-negative, and
     * are returnede in descending order. The first min(m,n) columns of
     * U and V are the left and right singular vectors of A. (Note that
     * the routine returns V^T or V^H, not V.
     */

namespace detail {

      inline
      void gesdd (char const jobz, int const m, int const n,
                  float* a, int const lda,
                  float* s, float* u, int const ldu,
                  float* vt, int const ldvt,
                  float* work, int const lwork, float* /* dummy */,
                  int* iwork, int* info)
      {
        LAPACK_SGESDD (&jobz, &m, &n, a, &lda, s,
                       u, &ldu, vt, &ldvt, work, &lwork, iwork, info);
      }

      inline
      void gesdd (char const jobz, int const m, int const n,
                  double* a, int const lda,
                  double* s, double* u, int const ldu,
                  double* vt, int const ldvt,
                  double* work, int const lwork, double* /* dummy */,
                  int* iwork, int* info)
      {
        LAPACK_DGESDD (&jobz, &m, &n, a, &lda, s,
                       u, &ldu, vt, &ldvt, work, &lwork, iwork, info);
      }

      inline
      void gesdd (char const jobz, int const m, int const n,
                  traits::complex_f* a, int const lda,
                  float* s, traits::complex_f* u, int const ldu,
                  traits::complex_f* vt, int const ldvt,
                  traits::complex_f* work, int const lwork,
                  float* rwork, int* iwork, int* info)
      {
        LAPACK_CGESDD (&jobz, &m, &n,
                       traits::complex_ptr (a), &lda, s,
                       traits::complex_ptr (u), &ldu,
                       traits::complex_ptr (vt), &ldvt,
                       traits::complex_ptr (work), &lwork,
                       rwork, iwork, info);
      }

      inline
      void gesdd (char const jobz, int const m, int const n,
                  traits::complex_d* a, int const lda,
                  double* s, traits::complex_d* u, int const ldu,
                  traits::complex_d* vt, int const ldvt,
                  traits::complex_d* work, int const lwork,
                  double* rwork, int* iwork, int* info)
      {
        LAPACK_ZGESDD (&jobz, &m, &n,
                       traits::complex_ptr (a), &lda, s,
                       traits::complex_ptr (u), &ldu,
                       traits::complex_ptr (vt), &ldvt,
                       traits::complex_ptr (work), &lwork,
                       rwork, iwork, info);
      }

      inline
      int gesdd_min_work (float, char jobz, int m, int n) {
        int minmn = m < n ? m : n;
        int maxmn = m < n ? n : m;
        int m3 = 3 * minmn;
        int m4 = 4 * minmn;
        int minw;
        if (jobz == 'N') {
          // leading comments:
          // LWORK >= 3*min(M,N) + max(max(M,N), 6*min(M,N))
          // code:
          // LWORK >= 3*min(M,N) + max(max(M,N), 7*min(M,N))
          int m7 = 7 * minmn;
          minw = maxmn < m7 ? m7 : maxmn;
          minw += m3;
        }
        if (jobz == 'O') {
          // LWORK >= 3*min(M,N)*min(M,N)
          // + max(max(M,N), 5*min(M,N)*min(M,N)+4*min(M,N))
          int m5 = 5 * minmn * minmn + m4;
          minw = maxmn < m5 ? m5 : maxmn;
          minw += m3 * minmn;
        }
        if (jobz == 'S' || jobz == 'A') {
          // LWORK >= 3*min(M,N)*min(M,N)
          // + max(max(M,N), 4*min(M,N)*min(M,N)+4*min(M,N)).
          int m44 = m4 * minmn + m4;
          minw = maxmn < m44 ? m44 : maxmn;
          minw += m3 * minmn;
        }
        return minw;
      }
      inline
      int gesdd_min_work (double, char jobz, int m, int n) {
        int minmn = m < n ? m : n;
        int maxmn = m < n ? n : m;
        int m3 = 3 * minmn;
        int m4 = 4 * minmn;
        int minw;
        assert( jobz=='N' || jobz=='O' || jobz=='S' || jobz=='A' ) ;
        if (jobz == 'N') {
          // leading comments:
          // LWORK >= 3*min(M,N) + max(max(M,N), 6*min(M,N))
          // code:
          // LWORK >= 3*min(M,N) + max(max(M,N), 7*min(M,N))
          int m7 = 7 * minmn;
          minw = maxmn < m7 ? m7 : maxmn;
          minw += m3;
        }
        else if (jobz == 'O') {
          // LWORK >= 3*min(M,N)*min(M,N)
          // + max(max(M,N), 5*min(M,N)*min(M,N)+4*min(M,N))
          int m5 = 5 * minmn * minmn + m4;
          minw = maxmn < m5 ? m5 : maxmn;
          minw += m3 * minmn;
        }
        else if (jobz == 'S' || jobz == 'A') {
          // LWORK >= 3*min(M,N)*min(M,N)
          // + max(max(M,N), 4*min(M,N)*min(M,N)+4*min(M,N)).
          int m44 = m4 * minmn + m4;
          minw = maxmn < m44 ? m44 : maxmn;
          minw += m3 * minmn;
        }
        return minw;
      }
      inline
      int gesdd_min_work (traits::complex_f, char jobz, int m, int n) {
        int minmn = m < n ? m : n;
        int maxmn = m < n ? n : m;
        int m2 = 2 * minmn;
        int minw = m2 + maxmn;
        if (jobz == 'N')
          // LWORK >= 2*min(M,N)+max(M,N)
          ;
        if (jobz == 'O')
          // LWORK >= 2*min(M,N)*min(M,N) + 2*min(M,N) + max(M,N)
          minw += m2 * minmn;
        if (jobz == 'S' || jobz == 'A')
          // LWORK >= min(M,N)*min(M,N) + 2*min(M,N) + max(M,N)
          minw += minmn * minmn;
        return minw;
      }
      inline
      int gesdd_min_work (traits::complex_d, char jobz, int m, int n) {
        int minmn = m < n ? m : n;
        int maxmn = m < n ? n : m;
        int m2 = 2 * minmn;
        int minw = m2 + maxmn;
        if (jobz == 'N')
          // LWORK >= 2*min(M,N)+max(M,N)
          ;
        if (jobz == 'O')
          // LWORK >= 2*min(M,N)*min(M,N) + 2*min(M,N) + max(M,N)
          minw += m2 * minmn;
        if (jobz == 'S' || jobz == 'A')
          // LWORK >= min(M,N)*min(M,N) + 2*min(M,N) + max(M,N)
          minw += minmn * minmn;
        return minw;
      }

      inline
      int gesdd_rwork (float, char, int, int) { return 1; }
      inline
      int gesdd_rwork (double, char, int, int) { return 1; }
      inline
      int gesdd_rwork (traits::complex_f, char jobz, int m, int n) {
        int minmn = m < n ? m : n;
        int minw;
        if (jobz == 'N')
          // LWORK >= 7*min(M,N)
          minw = 7 * minmn;
        else
          // LRWORK >= 5*min(M,N)*min(M,N) + 5*min(M,N)
          minw = 5 * (minmn * minmn + minmn);
        return minw;
      }
      inline
      int gesdd_rwork (traits::complex_d, char jobz, int m, int n) {
        int minmn = m < n ? m : n;
        int minw;
        if (jobz == 'N')
          // LWORK >= 7*min(M,N)
          minw = 7 * minmn;
        else
          // LRWORK >= 5*min(M,N)*min(M,N) + 5*min(M,N)
          minw = 5 * (minmn * minmn + minmn);
        return minw;
      }

      inline
      int gesdd_iwork (int m, int n) {
        int minmn = m < n ? m : n;
        return 8 * minmn;
      }

} // detail

    template <typename MatrA>
    inline
    int gesdd_work (char const q, char const jobz, MatrA const& a)
    {

#ifdef BOOST_NUMERIC_BINDINGS_LAPACK_2
      assert (q == 'M');
#else
      assert (q == 'M' || q == 'O');
#endif
      assert (jobz == 'N' || jobz == 'O' || jobz == 'A' || jobz == 'S');

      int const m = traits::matrix_size1 (a);
      int const n = traits::matrix_size2 (a);
      int lw = -13;

#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
typedef typename traits::matrix_traits<MatrA>::value_type val_t;
#else
typedef typename MatrA::value_type val_t;
#endif

if (q == 'M')
lw = detail::gesdd_min_work (val_t(), jobz, m, n);

#ifndef BOOST_NUMERIC_BINDINGS_LAPACK_2
      MatrA& a2 = const_cast<MatrA&> (a);
      if (q == 'O') {
        // traits::detail::array<val_t> w (1);
        val_t w;
        int info;
        detail::gesdd (jobz, m, n,
                       traits::matrix_storage (a2),
                       traits::leading_dimension (a2),
                       0, // traits::vector_storage (s),
                       0, // traits::matrix_storage (u),
                       m, // traits::leading_dimension (u),
                       0, // traits::matrix_storage (vt),
                       n, // traits::leading_dimension (vt),
                       &w, // traits::vector_storage (w),
                       -1, // traits::vector_size (w),
                       0, // traits::vector_storage (rw),
                       0, // traits::vector_storage (iw),
                       &info);
        assert (info == 0);

        lw = traits::detail::to_int (w);
        // // lw = traits::detail::to_int (w[0]);
        /*
         * is there a bug in LAPACK? or in Mandrake's .rpm?
         * if m == 3, n == 4 and jobz == 'N' (real A),
         * gesdd() returns optimal size == 1 while minimum size == 27
         */
        // int lwo = traits::detail::to_int (w);
        // int lwmin = detail::gesdd_min_work (val_t(), jobz, m, n);
        // lw = lwo < lwmin ? lwmin : lwo;
      }
#endif

      return lw;
    }

    template <typename MatrA>
    inline
    int gesdd_rwork (char jobz, MatrA const& a) {

assert (jobz == 'N' || jobz == 'O' || jobz == 'A' || jobz == 'S');

#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
typedef typename traits::matrix_traits<MatrA>::value_type val_t;
#else
typedef typename MatrA::value_type val_t;
#endif

      return detail::gesdd_rwork (val_t(), jobz,
                                  traits::matrix_size1 (a),
                                  traits::matrix_size2 (a));
    }

    template <typename MatrA>
    inline
    int gesdd_iwork (MatrA const& a) {

      return detail::gesdd_iwork (traits::matrix_size1 (a),
                                  traits::matrix_size2 (a));
    }

    template <typename MatrA, typename VecS,
              typename MatrU, typename MatrV, typename VecW, typename VecIW>
    inline
    int gesdd (char const jobz, MatrA& a,
               VecS& s, MatrU& u, MatrV& vt, VecW& w, VecIW& iw)
    {

#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
      BOOST_STATIC_ASSERT((boost::is_same<
        typename traits::matrix_traits<MatrA>::matrix_structure,
        traits::general_t
>::value));
      BOOST_STATIC_ASSERT((boost::is_same<
        typename traits::matrix_traits<MatrU>::matrix_structure,
        traits::general_t
>::value));
      BOOST_STATIC_ASSERT((boost::is_same<
        typename traits::matrix_traits<MatrV>::matrix_structure,
        traits::general_t
>::value));

      BOOST_STATIC_ASSERT(
        (boost::is_same<
          typename traits::matrix_traits<MatrA>::value_type, float
>::value
        ||
        boost::is_same<
          typename traits::matrix_traits<MatrA>::value_type, double
>::value));
#endif

      int const m = traits::matrix_size1 (a);
      int const n = traits::matrix_size2 (a);
      int const minmn = m < n ? m : n;

      assert (minmn == traits::vector_size (s));
      assert ((jobz == 'N')
              || ((jobz == 'O' || jobz == 'A') && m >= n)
              || ((jobz == 'O' || jobz == 'A')
                  && m < n
                  && m == traits::matrix_size2 (u))
              || (jobz == 'S' && minmn == traits::matrix_size2 (u)));
      assert ((jobz == 'N' && traits::leading_dimension (u) >= 1)
              || (jobz == 'O'
                  && m >= n
                  && traits::leading_dimension (u) >= 1)
              || (jobz == 'O'
                  && m < n
                  && traits::leading_dimension (u) >= m)
              || (jobz == 'A' && traits::leading_dimension (u) >= m)
              || (jobz == 'S' && traits::leading_dimension (u) >= m));
      assert (n == traits::matrix_size2 (vt));
      assert ((jobz == 'N' && traits::leading_dimension (vt) >= 1)
              || (jobz == 'O'
                  && m < n
                  && traits::leading_dimension (vt) >= 1)
              || (jobz == 'O'
                  && m >= n
                  && traits::leading_dimension (vt) >= n)
              || (jobz == 'A' && traits::leading_dimension (vt) >= n)
              || (jobz == 'S' && traits::leading_dimension (vt) >= minmn));
#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
      typedef typename traits::matrix_traits<MatrA>::value_type val_t;
#else
      typedef typename MatrA::value_type val_t;
#endif
      assert (traits::vector_size (w)
>= detail::gesdd_min_work (val_t(), jobz, m, n));
      assert (traits::vector_size (iw) >= detail::gesdd_iwork (m, n));

      int info;
      detail::gesdd (jobz, m, n,
                     traits::matrix_storage (a),
                     traits::leading_dimension (a),
                     traits::vector_storage (s),
                     traits::matrix_storage (u),
                     traits::leading_dimension (u),
                     traits::matrix_storage (vt),
                     traits::leading_dimension (vt),
                     traits::vector_storage (w),
                     traits::vector_size (w),
                     0, // dummy argument
                     traits::vector_storage (iw),
                     &info);
      return info;
    }

    template <typename MatrA, typename VecS, typename MatrU,
              typename MatrV, typename VecW, typename VecRW, typename VecIW>
    inline
    int gesdd (char const jobz, MatrA& a,
               VecS& s, MatrU& u, MatrV& vt, VecW& w, VecRW& rw, VecIW& iw)
    {

      int const m = traits::matrix_size1 (a);
      int const n = traits::matrix_size2 (a);
      int const minmn = m < n ? m : n;

      int info;
      detail::gesdd (jobz, m, n,
                     traits::matrix_storage (a),
                     traits::leading_dimension (a),
                     traits::vector_storage (s),
                     traits::matrix_storage (u),
                     traits::leading_dimension (u),
                     traits::matrix_storage (vt),
                     traits::leading_dimension (vt),
                     traits::vector_storage (w),
                     traits::vector_size (w),
                     traits::vector_storage (rw),
                     traits::vector_storage (iw),
                     &info);
      return info;
    }

    template <typename MatrA, typename VecS, typename MatrU, typename MatrV>
    inline
    int gesdd (char const opt, char const jobz,
               MatrA& a, VecS& s, MatrU& u, MatrV& vt)
    {

      int const m = traits::matrix_size1 (a);
      int const n = traits::matrix_size2 (a);
#ifndef NDEBUG
      int const minmn = m < n ? m : n;
#endif // NDEBUG

#ifdef BOOST_NUMERIC_BINDINGS_LAPACK_2
assert (opt == 'M');
#else
assert (opt == 'M' || opt == 'O');
#endif

#ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
      typedef typename traits::matrix_traits<MatrA>::value_type val_t;
#else
      typedef typename MatrA::value_type val_t;
#endif
      typedef typename traits::type_traits<val_t>::real_type real_t;

      int const lw = gesdd_work (opt, jobz, a);
      traits::detail::array<val_t> w (lw);
      if (!w.valid()) return -101;

      int const lrw = gesdd_rwork (jobz, a);
      traits::detail::array<real_t> rw (lrw);
      if (!rw.valid()) return -102;

      int const liw = gesdd_iwork (a);
      traits::detail::array<int> iw (liw);
      if (!iw.valid()) return -103;

      int info;
      detail::gesdd (jobz, m, n,
                     traits::matrix_storage (a),
                     traits::leading_dimension (a),
                     traits::vector_storage (s),
                     traits::matrix_storage (u),
                     traits::leading_dimension (u),
                     traits::matrix_storage (vt),
                     traits::leading_dimension (vt),
                     traits::vector_storage (w),
                     lw, //traits::vector_size (w),
                     traits::vector_storage (rw),
                     traits::vector_storage (iw),
                     &info);
      return info;
    }

#ifndef BOOST_NUMERIC_BINDINGS_LAPACK_2

    template <typename MatrA, typename VecS, typename MatrU, typename MatrV>
    inline
    int gesdd (char const jobz, MatrA& a, VecS& s, MatrU& u, MatrV& vt) {
      return gesdd ('O', jobz, a, s, u, vt);
    }

    template <typename MatrA, typename VecS, typename MatrU, typename MatrV>
    inline
    int gesdd (MatrA& a, VecS& s, MatrU& u, MatrV& vt) {
      return gesdd ('O', 'S', a, s, u, vt);
    }

    template <typename MatrA, typename VecS>
    inline
    int gesdd (MatrA& a, VecS& s) {

      int const m = traits::matrix_size1 (a);
      int const n = traits::matrix_size2 (a);
      int const minmn = m < n ? m : n;

assert (minmn == traits::vector_size (s));

      int const lw = gesdd_work ('O', 'N', a);
      traits::detail::array<val_t> w (lw);
      if (!w.valid()) return -101;

      int const lrw = gesdd_rwork ('N', a);
      traits::detail::array<real_t> rw (lrw);
      if (!rw.valid()) return -102;

      int const liw = gesdd_iwork (a);
      traits::detail::array<int> iw (liw);
      if (!iw.valid()) return -103;

      int info;
      detail::gesdd ('N', m, n,
                     traits::matrix_storage (a),
                     traits::leading_dimension (a),
                     traits::vector_storage (s),
                     0, // traits::matrix_storage (u),
                     1, // traits::leading_dimension (u),
                     0, // traits::matrix_storage (vt),
                     1, // traits::leading_dimension (vt),
                     traits::vector_storage (w),
                     traits::vector_size (w),
                     traits::vector_storage (rw),
                     traits::vector_storage (iw),
                     &info);
      return info;
    }

#endif

} // namespace lapack

}}}

#endif

Next message: Thomas Lemaire: "Re: [ublas] minor bug in lapack bindings gesdd"
Previous message: Thomas Lemaire: "Re: [ublas] minor bug in lapack bindings gesdd"
In reply to: Thomas Lemaire: "Re: [ublas] minor bug in lapack bindings gesdd"
Next in thread: Thomas Lemaire: "Re: [ublas] minor bug in lapack bindings gesdd"
Reply: Thomas Lemaire: "Re: [ublas] minor bug in lapack bindings gesdd"

Date view	Thread view	Subject view	Author view