Mercurial > octave

view liboctave/numeric/svd.cc @ 22204:469c817eb256

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
svd: reduce code duplication with more use of template and macro. * liboctave/numeric/svd.cc, liboctave/numeric/svd.h: remove unused constructor with reference for int (info). This allows to move all of the constructor into a single template, so remove init(). Two new methods, gesvd and gesdd, are fully specialized but the main hunck of code are the long list of arguments. Scope type and drive enums to the svd class for clarity, and rename member names. Add a new member for the drive used. * libinterp/corefcn/svd.cc: fix typenames for the svd enums which are now scoped. * CMatrix.cc, dMatrix.cc, fCMatrix.cc, fMatrix.cc: fix typenames for the svd enums which are now scoped.
author Carnë Draug <carandraug@octave.org>
date Thu, 04 Aug 2016 20:20:27 +0100
parents 407c66ae1e20
children 4afe3705ea75
line wrap: on
line source

/*

Copyright (C) 2016 Carnë Draug
Copyright (C) 1994-2015 John W. Eaton

This file is part of Octave.

Octave is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

Octave is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING.  If not, see
<http://www.gnu.org/licenses/>.

*/

#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

#include "svd.h"

#include <cassert>
#include <algorithm>

#include "CMatrix.h"
#include "dDiagMatrix.h"
#include "fDiagMatrix.h"
#include "dMatrix.h"
#include "f77-fcn.h"
#include "fCMatrix.h"
#include "fMatrix.h"
#include "lo-error.h"
#include "oct-locbuf.h"

extern "C"
{
  F77_RET_T
  F77_FUNC (dgesvd, DGESVD) (F77_CONST_CHAR_ARG_DECL,
                             F77_CONST_CHAR_ARG_DECL,
                             const F77_INT&, const F77_INT&,
                             F77_DBLE*, const F77_INT&, F77_DBLE*,
                             F77_DBLE*, const F77_INT&, F77_DBLE*,
                             const F77_INT&, F77_DBLE*,
                             const F77_INT&, F77_INT&
                             F77_CHAR_ARG_LEN_DECL
                             F77_CHAR_ARG_LEN_DECL);

  F77_RET_T
  F77_FUNC (dgesdd, DGESDD) (F77_CONST_CHAR_ARG_DECL,
                             const F77_INT&, const F77_INT&,
                             F77_DBLE*, const F77_INT&, F77_DBLE*,
                             F77_DBLE*, const F77_INT&, F77_DBLE*,
                             const F77_INT&, F77_DBLE*,
                             const F77_INT&, F77_INT *,
                             F77_INT&
                             F77_CHAR_ARG_LEN_DECL);

  F77_RET_T
  F77_FUNC (sgesvd, SGESVD) (F77_CONST_CHAR_ARG_DECL,
                             F77_CONST_CHAR_ARG_DECL,
                             const F77_INT&, const F77_INT&,
                             F77_REAL*, const F77_INT&, F77_REAL*,
                             F77_REAL*, const F77_INT&, F77_REAL*,
                             const F77_INT&, F77_REAL*,
                             const F77_INT&, F77_INT&
                             F77_CHAR_ARG_LEN_DECL
                             F77_CHAR_ARG_LEN_DECL);

  F77_RET_T
  F77_FUNC (sgesdd, SGESDD) (F77_CONST_CHAR_ARG_DECL,
                             const F77_INT&, const F77_INT&,
                             F77_REAL*, const F77_INT&, F77_REAL*,
                             F77_REAL*, const F77_INT&, F77_REAL*,
                             const F77_INT&, F77_REAL*,
                             const F77_INT&, F77_INT *,
                             F77_INT&
                             F77_CHAR_ARG_LEN_DECL);

  F77_RET_T
  F77_FUNC (zgesvd, ZGESVD) (F77_CONST_CHAR_ARG_DECL,
                             F77_CONST_CHAR_ARG_DECL,
                             const F77_INT&, const F77_INT&,
                             F77_DBLE_CMPLX*, const F77_INT&,
                             F77_DBLE*, F77_DBLE_CMPLX*, const F77_INT&,
                             F77_DBLE_CMPLX*, const F77_INT&, F77_DBLE_CMPLX*,
                             const F77_INT&, F77_DBLE*, F77_INT&
                             F77_CHAR_ARG_LEN_DECL
                             F77_CHAR_ARG_LEN_DECL);

  F77_RET_T
  F77_FUNC (zgesdd, ZGESDD) (F77_CONST_CHAR_ARG_DECL,
                             const F77_INT&, const F77_INT&,
                             F77_DBLE_CMPLX*, const F77_INT&,
                             F77_DBLE*, F77_DBLE_CMPLX*, const F77_INT&,
                             F77_DBLE_CMPLX*, const F77_INT&, F77_DBLE_CMPLX*,
                             const F77_INT&, F77_DBLE*,
                             F77_INT *, F77_INT&
                             F77_CHAR_ARG_LEN_DECL);
  F77_RET_T
  F77_FUNC (cgesvd, CGESVD) (F77_CONST_CHAR_ARG_DECL,
                             F77_CONST_CHAR_ARG_DECL,
                             const F77_INT&, const F77_INT&,
                             F77_CMPLX*, const F77_INT&, F77_REAL*,
                             F77_CMPLX*, const F77_INT&,
                             F77_CMPLX*, const F77_INT&,
                             F77_CMPLX*, const F77_INT&,
                             F77_REAL*, F77_INT&
                             F77_CHAR_ARG_LEN_DECL
                             F77_CHAR_ARG_LEN_DECL);

  F77_RET_T
  F77_FUNC (cgesdd, CGESDD) (F77_CONST_CHAR_ARG_DECL,
                             const F77_INT&, const F77_INT&,
                             F77_CMPLX*, const F77_INT&, F77_REAL*,
                             F77_CMPLX*, const F77_INT&,
                             F77_CMPLX*, const F77_INT&,
                             F77_CMPLX*, const F77_INT&,
                             F77_REAL*, F77_INT *, F77_INT&
                             F77_CHAR_ARG_LEN_DECL);
}

template <typename T>
T
svd<T>::left_singular_matrix (void) const
{
  if (type == svd::Type::sigma_only)
    (*current_liboctave_error_handler)
      ("svd: U not computed because type == svd::sigma_only");

  return left_sm;
}

template <typename T>
T
svd<T>::right_singular_matrix (void) const
{
  if (type == svd::Type::sigma_only)
    (*current_liboctave_error_handler)
      ("svd: V not computed because type == svd::sigma_only");

  return right_sm;
}


// GESVD specializations

#define GESVD_REAL_STEP(f, F)                                 \
  F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobu, 1),            \
                   F77_CONST_CHAR_ARG2 (&jobv, 1),            \
                   m, n, tmp_data, m1, s_vec, u, m1, vt,      \
                   nrow_vt1, work.fortran_vec (), lwork, info \
                   F77_CHAR_ARG_LEN (1)                       \
                   F77_CHAR_ARG_LEN (1)))

#define GESVD_COMPLEX_STEP(f, F, CMPLX_ARG)           \
  F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobu, 1),    \
                   F77_CONST_CHAR_ARG2 (&jobv, 1),    \
                   m, n, CMPLX_ARG (tmp_data),        \
                   m1, s_vec, CMPLX_ARG (u), m1,      \
                   CMPLX_ARG (vt), nrow_vt1,          \
                   CMPLX_ARG (work.fortran_vec ()),   \
                   lwork, rwork.fortran_vec (), info  \
                   F77_CHAR_ARG_LEN (1)               \
                   F77_CHAR_ARG_LEN (1)))

// DGESVD
template<>
void
svd<Matrix>::gesvd (char& jobu, char& jobv, octave_idx_type m,
                    octave_idx_type n, double* tmp_data, octave_idx_type m1,
                    double* s_vec, double* u, double* vt,
                    octave_idx_type nrow_vt1, Matrix& work,
                    octave_idx_type& lwork, octave_idx_type& info)
{
  GESVD_REAL_STEP (dgesvd, DGESVD);

  lwork = static_cast<octave_idx_type> (work(0));
  work.resize (lwork, 1);

  GESVD_REAL_STEP (dgesvd, DGESVD);
}

// SGESVD
template<>
void
svd<FloatMatrix>::gesvd (char& jobu, char& jobv, octave_idx_type m,
                         octave_idx_type n, float* tmp_data,
                         octave_idx_type m1, float* s_vec, float* u, float* vt,
                         octave_idx_type nrow_vt1, FloatMatrix& work,
                         octave_idx_type& lwork, octave_idx_type& info)
{
  GESVD_REAL_STEP (sgesvd, SGESVD);

  lwork = static_cast<octave_idx_type> (work(0));
  work.resize (lwork, 1);

  GESVD_REAL_STEP (sgesvd, SGESVD);
}

// ZGESVD
template<>
void
svd<ComplexMatrix>::gesvd (char& jobu, char& jobv, octave_idx_type m,
                           octave_idx_type n, Complex* tmp_data,
                           octave_idx_type m1, double* s_vec, Complex* u,
                           Complex* vt, octave_idx_type nrow_vt1,
                           ComplexMatrix& work,
                           octave_idx_type& lwork, octave_idx_type& info)
{
  Matrix rwork (5 * std::max (m, n), 1);

  GESVD_COMPLEX_STEP (zgesvd, ZGESVD, F77_DBLE_CMPLX_ARG);

  lwork = static_cast<octave_idx_type> (work(0).real ());
  work.resize (lwork, 1);

  GESVD_COMPLEX_STEP (zgesvd, ZGESVD, F77_DBLE_CMPLX_ARG);
}

// CGESVD
template<>
void
svd<FloatComplexMatrix>::gesvd (char& jobu, char& jobv,
                                octave_idx_type m, octave_idx_type n,
                                FloatComplex* tmp_data, octave_idx_type m1,
                                float* s_vec, FloatComplex* u,
                                FloatComplex* vt, octave_idx_type nrow_vt1,
                                FloatComplexMatrix& work,
                                octave_idx_type& lwork, octave_idx_type& info)
{
  FloatMatrix rwork (5 * std::max (m, n), 1);

  GESVD_COMPLEX_STEP (cgesvd, CGESVD, F77_CMPLX_ARG);

  lwork = static_cast<octave_idx_type> (work(0).real ());
  work.resize (lwork, 1);

  GESVD_COMPLEX_STEP (cgesvd, CGESVD, F77_CMPLX_ARG);
}

#undef GESVD_REAL_STEP
#undef GESVD_COMPLEX_STEP


// GESDD specializations

#define GESDD_REAL_STEP(f, F)                                       \
  F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobz, 1),                  \
                   m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1,  \
                   work.fortran_vec (), lwork, iwork, info          \
                   F77_CHAR_ARG_LEN (1)))

#define GESDD_COMPLEX_STEP(f, F, CMPLX_ARG)                 \
  F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobz, 1), m, n,    \
                   CMPLX_ARG (tmp_data), m1,                \
                   s_vec, CMPLX_ARG (u), m1,                \
                   CMPLX_ARG (vt), nrow_vt1,                \
                   CMPLX_ARG (work.fortran_vec ()), lwork,  \
                   rwork.fortran_vec (), iwork, info        \
                   F77_CHAR_ARG_LEN (1)))

// DGESDD
template<>
void
svd<Matrix>::gesdd (char& jobz, octave_idx_type m, octave_idx_type n,
                    double* tmp_data, octave_idx_type m1,
                    double* s_vec, double* u,
                    double* vt, octave_idx_type nrow_vt1,
                    Matrix& work, octave_idx_type& lwork,
                    octave_idx_type* iwork, octave_idx_type& info)
{
  GESDD_REAL_STEP (dgesdd, DGESDD);

  lwork = static_cast<octave_idx_type> (work(0));
  work.resize (lwork, 1);

  GESDD_REAL_STEP (dgesdd, DGESDD);
}

// SGESDD
template<>
void
svd<FloatMatrix>::gesdd (char& jobz, octave_idx_type m, octave_idx_type n,
                         float* tmp_data, octave_idx_type m1,
                         float* s_vec, float* u,
                         float* vt, octave_idx_type nrow_vt1,
                         FloatMatrix& work, octave_idx_type& lwork,
                         octave_idx_type* iwork, octave_idx_type& info)
{
  GESDD_REAL_STEP (sgesdd, SGESDD);

  lwork = static_cast<octave_idx_type> (work(0));
  work.resize (lwork, 1);

  GESDD_REAL_STEP (sgesdd, SGESDD);
}

// ZGESDD
template<>
void
svd<ComplexMatrix>::gesdd (char& jobz, octave_idx_type m, octave_idx_type n,
                           Complex* tmp_data, octave_idx_type m1,
                           double* s_vec, Complex* u,
                           Complex* vt, octave_idx_type nrow_vt1,
                           ComplexMatrix& work, octave_idx_type& lwork,
                           octave_idx_type* iwork, octave_idx_type& info)
{

  octave_idx_type min_mn = std::min (m, n);

  octave_idx_type lrwork;
  if (jobz == 'N')
    lrwork = 7*min_mn;
  else
    lrwork = 5*min_mn*min_mn + 5*min_mn;
  Matrix rwork (lrwork, 1);


  GESDD_COMPLEX_STEP (zgesdd, ZGESDD, F77_DBLE_CMPLX_ARG);

  lwork = static_cast<octave_idx_type> (work(0).real ());
  work.resize (lwork, 1);

  GESDD_COMPLEX_STEP (zgesdd, ZGESDD, F77_DBLE_CMPLX_ARG);
}

// CGESDD
template<>
void
svd<FloatComplexMatrix>::gesdd (char& jobz, octave_idx_type m,
                                octave_idx_type n,
                                FloatComplex* tmp_data, octave_idx_type m1,
                                float* s_vec, FloatComplex* u,
                                FloatComplex* vt, octave_idx_type nrow_vt1,
                                FloatComplexMatrix& work,
                                octave_idx_type& lwork, octave_idx_type* iwork,
                                octave_idx_type& info)
{
  octave_idx_type min_mn = std::min (m, n);
  octave_idx_type max_mn = std::max (m, n);

  octave_idx_type lrwork;
  if (jobz == 'N')
    lrwork = 5*min_mn;
  else
    lrwork = min_mn * std::max (5*min_mn+7, 2*max_mn+2*min_mn+1);
  FloatMatrix rwork (lrwork, 1);

  GESDD_COMPLEX_STEP (cgesdd, CGESDD, F77_CMPLX_ARG);

  lwork = static_cast<octave_idx_type> (work(0).real ());
  work.resize (lwork, 1);

  GESDD_COMPLEX_STEP (cgesdd, CGESDD, F77_CMPLX_ARG);
}

#undef GESDD_REAL_STEP
#undef GESDD_COMPLEX_STEP


template<typename T>
svd<T>::svd (const T& a, svd::Type type,
             svd::Driver driver)
  : type (type), driver (driver), left_sm (), sigma (), right_sm ()
{
  octave_idx_type info;

  octave_idx_type m = a.rows ();
  octave_idx_type n = a.cols ();

  if (m == 0 || n == 0)
    {
      switch (type)
        {
        case svd::Type::std:
          left_sm = T (m, m, 0);
          for (octave_idx_type i = 0; i < m; i++)
            left_sm.xelem (i, i) = 1;
          sigma = DM_T (m, n);
          right_sm = T (n, n, 0);
          for (octave_idx_type i = 0; i < n; i++)
            right_sm.xelem (i, i) = 1;
          break;

        case svd::Type::economy:
          left_sm = T (m, 0, 0);
          sigma = DM_T (0, 0);
          right_sm = T (0, n, 0);
          break;

        case svd::Type::sigma_only:
        default:
          sigma = DM_T (0, 1);
          break;
        }
      return;
    }

  T atmp = a;
  P* tmp_data = atmp.fortran_vec ();

  octave_idx_type min_mn = m < n ? m : n;

  char jobu = 'A';
  char jobv = 'A';

  octave_idx_type ncol_u = m;
  octave_idx_type nrow_vt = n;
  octave_idx_type nrow_s = m;
  octave_idx_type ncol_s = n;

  switch (type)
    {
    case svd::Type::economy:
      jobu = jobv = 'S';
      ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
      break;

    case svd::Type::sigma_only:

      // Note:  for this case, both jobu and jobv should be 'N', but
      // there seems to be a bug in dgesvd from Lapack V2.0.  To
      // demonstrate the bug, set both jobu and jobv to 'N' and find
      // the singular values of [eye(3), eye(3)].  The result is
      // [-sqrt(2), -sqrt(2), -sqrt(2)].
      //
      // For Lapack 3.0, this problem seems to be fixed.

      jobu = jobv = 'N';
      ncol_u = nrow_vt = 1;
      break;

    default:
      break;
    }

  if (! (jobu == 'N' || jobu == 'O'))
    left_sm.resize (m, ncol_u);

  P* u = left_sm.fortran_vec ();

  sigma.resize (nrow_s, ncol_s);
  DM_P* s_vec = sigma.fortran_vec ();

  if (! (jobv == 'N' || jobv == 'O'))
    right_sm.resize (nrow_vt, n);

  P* vt = right_sm.fortran_vec ();

  // Query _GESVD for the correct dimension of WORK.

  octave_idx_type lwork = -1;

  T work (1, 1);

  octave_idx_type m1 = std::max (m, 1);
  octave_idx_type nrow_vt1 = std::max (nrow_vt, 1);

  if (driver == svd::Driver::GESVD)
    gesvd (jobu, jobv, m, n, tmp_data, m1, s_vec, u, vt, nrow_vt1,
           work, lwork, info);
  else if (driver == svd::Driver::GESDD)
    {
      assert (jobu == jobv);
      char jobz = jobu;

      OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8 * std::min (m, n));

      gesdd (jobz, m, n, tmp_data, m1, s_vec, u, vt, nrow_vt1,
             work, lwork, iwork, info);
    }
  else
    abort ();

  if (! (jobv == 'N' || jobv == 'O'))
    right_sm = right_sm.transpose ();
}

// Instantiations we need.

template class svd<Matrix>;

template class svd<FloatMatrix>;

template class svd<ComplexMatrix>;

template class svd<FloatComplexMatrix>;