# HG changeset patch
# User Carnë Draug <carandraug@octave.org>
# Date 1470338427 -3600
# Node ID 469c817eb2560e93208a10b267c35c4bccd46bf3
# Parent  2f301a30aeedcc36a3358d59f3090c5b7e071f15
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.

diff -r 2f301a30aeed -r 469c817eb256 libinterp/corefcn/svd.cc
--- a/libinterp/corefcn/svd.cc	Wed Aug 03 10:23:41 2016 -0700
+++ b/libinterp/corefcn/svd.cc	Thu Aug 04 20:20:27 2016 +0100
@@ -37,19 +37,19 @@
 static std::string Vsvd_driver = "gesvd";
 
 template <typename T>
-static typename svd<T>::type
+static typename svd<T>::Type
 svd_type (int nargin, int nargout)
 {
   return ((nargout == 0 || nargout == 1)
-          ? svd<T>::sigma_only
-          : (nargin == 2) ? svd<T>::economy : svd<T>::std);
+          ? svd<T>::Type::sigma_only
+          : (nargin == 2) ? svd<T>::Type::economy : svd<T>::Type::std);
 }
 
 template <typename T>
-static typename svd<T>::driver
+static typename svd<T>::Driver
 svd_driver (void)
 {
-  return Vsvd_driver == "gesvd" ? svd<T>::GESVD : svd<T>::GESDD;
+  return Vsvd_driver == "gesvd" ? svd<T>::Driver::GESVD : svd<T>::Driver::GESDD;
 }
 
 DEFUN (svd, args, nargout,
diff -r 2f301a30aeed -r 469c817eb256 liboctave/array/CMatrix.cc
--- a/liboctave/array/CMatrix.cc	Wed Aug 03 10:23:41 2016 -0700
+++ b/liboctave/array/CMatrix.cc	Thu Aug 04 20:20:27 2016 +0100
@@ -1139,7 +1139,7 @@
 {
   ComplexMatrix retval;
 
-  svd<ComplexMatrix> result (*this, svd<ComplexMatrix>::economy);
+  svd<ComplexMatrix> result (*this, svd<ComplexMatrix>::Type::economy);
 
   DiagMatrix S = result.singular_values ();
   ComplexMatrix U = result.left_singular_matrix ();
diff -r 2f301a30aeed -r 469c817eb256 liboctave/array/dMatrix.cc
--- a/liboctave/array/dMatrix.cc	Wed Aug 03 10:23:41 2016 -0700
+++ b/liboctave/array/dMatrix.cc	Thu Aug 04 20:20:27 2016 +0100
@@ -823,7 +823,7 @@
 Matrix
 Matrix::pseudo_inverse (double tol) const
 {
-  svd<Matrix> result (*this, svd<Matrix>::economy);
+  svd<Matrix> result (*this, svd<Matrix>::Type::economy);
 
   DiagMatrix S = result.singular_values ();
   Matrix U = result.left_singular_matrix ();
diff -r 2f301a30aeed -r 469c817eb256 liboctave/array/fCMatrix.cc
--- a/liboctave/array/fCMatrix.cc	Wed Aug 03 10:23:41 2016 -0700
+++ b/liboctave/array/fCMatrix.cc	Thu Aug 04 20:20:27 2016 +0100
@@ -1144,7 +1144,7 @@
 {
   FloatComplexMatrix retval;
 
-  svd<FloatComplexMatrix> result (*this, svd<FloatComplexMatrix>::economy);
+  svd<FloatComplexMatrix> result (*this, svd<FloatComplexMatrix>::Type::economy);
 
   FloatDiagMatrix S = result.singular_values ();
   FloatComplexMatrix U = result.left_singular_matrix ();
diff -r 2f301a30aeed -r 469c817eb256 liboctave/array/fMatrix.cc
--- a/liboctave/array/fMatrix.cc	Wed Aug 03 10:23:41 2016 -0700
+++ b/liboctave/array/fMatrix.cc	Thu Aug 04 20:20:27 2016 +0100
@@ -830,7 +830,7 @@
 FloatMatrix
 FloatMatrix::pseudo_inverse (float tol) const
 {
-  svd<FloatMatrix> result (*this, svd<FloatMatrix>::economy);
+  svd<FloatMatrix> result (*this, svd<FloatMatrix>::Type::economy);
 
   FloatDiagMatrix S = result.singular_values ();
   FloatMatrix U = result.left_singular_matrix ();
diff -r 2f301a30aeed -r 469c817eb256 liboctave/numeric/oct-norm.cc
--- a/liboctave/numeric/oct-norm.cc	Wed Aug 03 10:23:41 2016 -0700
+++ b/liboctave/numeric/oct-norm.cc	Thu Aug 04 20:20:27 2016 +0100
@@ -485,7 +485,7 @@
   R res = 0;
   if (p == 2)
     {
-      svd<MatrixT> fact (m, svd<MatrixT>::sigma_only);
+      svd<MatrixT> fact (m, svd<MatrixT>::Type::sigma_only);
       res = fact.singular_values () (0,0);
     }
   else if (p == 1)
diff -r 2f301a30aeed -r 469c817eb256 liboctave/numeric/svd.cc
--- a/liboctave/numeric/svd.cc	Wed Aug 03 10:23:41 2016 -0700
+++ b/liboctave/numeric/svd.cc	Thu Aug 04 20:20:27 2016 +0100
@@ -1,5 +1,6 @@
 /*
 
+Copyright (C) 2016 Carnë Draug
 Copyright (C) 1994-2015 John W. Eaton
 
 This file is part of Octave.
@@ -24,7 +25,10 @@
 #  include "config.h"
 #endif
 
+#include "svd.h"
+
 #include <cassert>
+#include <algorithm>
 
 #include "CMatrix.h"
 #include "dDiagMatrix.h"
@@ -35,7 +39,6 @@
 #include "fMatrix.h"
 #include "lo-error.h"
 #include "oct-locbuf.h"
-#include "svd.h"
 
 extern "C"
 {
@@ -128,7 +131,7 @@
 T
 svd<T>::left_singular_matrix (void) const
 {
-  if (type_computed == svd::sigma_only)
+  if (type == svd::Type::sigma_only)
     (*current_liboctave_error_handler)
       ("svd: U not computed because type == svd::sigma_only");
 
@@ -139,65 +142,270 @@
 T
 svd<T>::right_singular_matrix (void) const
 {
-  if (type_computed == svd::sigma_only)
+  if (type == svd::Type::sigma_only)
     (*current_liboctave_error_handler)
       ("svd: V not computed because type == svd::sigma_only");
 
   return right_sm;
 }
 
-template <typename T>
-octave_idx_type
-svd<T>::empty_init (octave_idx_type nr, octave_idx_type nc, svd::type svd_type)
+
+// 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)
 {
-  assert (nr == 0 || nc == 0);
+  GESVD_REAL_STEP (dgesvd, DGESVD);
 
-  static typename T::element_type zero (0);
-  static typename T::element_type one (1);
+  lwork = static_cast<octave_idx_type> (work(0));
+  work.resize (lwork, 1);
+
+  GESVD_REAL_STEP (dgesvd, DGESVD);
+}
 
-  switch (svd_type)
-    {
-    case svd::std:
-      left_sm = T (nr, nr, zero);
-      for (octave_idx_type i = 0; i < nr; i++)
-        left_sm.xelem (i, i) = one;
-      sigma = DM_T (nr, nc);
-      right_sm = T (nc, nc, zero);
-      for (octave_idx_type i = 0; i < nc; i++)
-        right_sm.xelem (i, i) = one;
-      break;
+// 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);
+}
 
-    case svd::economy:
-      left_sm = T (nr, 0, zero);
-      sigma = DM_T (0, 0);
-      right_sm = T (0, nc, zero);
-      break;
+// 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);
+}
 
-    case svd::sigma_only:
-    default:
-      sigma = DM_T (0, 1);
-      break;
-    }
+// 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);
 
-  return 0;
+  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);
 }
 
-// Specializations.
+#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);
+}
 
-template <>
-octave_idx_type
-svd<Matrix>::init (const Matrix& a, svd::type svd_type, svd::driver svd_driver)
+// 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 info = 0;
+
+  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)
-    return empty_init (m, n, svd_type);
+    {
+      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;
 
-  Matrix atmp = a;
-  double *tmp_data = atmp.fortran_vec ();
+        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;
 
@@ -209,139 +417,14 @@
   octave_idx_type nrow_s = m;
   octave_idx_type ncol_s = n;
 
-  switch (svd_type)
+  switch (type)
     {
-    case svd::economy:
+    case svd::Type::economy:
       jobu = jobv = 'S';
       ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
       break;
 
-    case svd::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;
-    }
-
-  type_computed = svd_type;
-
-  if (! (jobu == 'N' || jobu == 'O'))
-    left_sm.resize (m, ncol_u);
-
-  double *u = left_sm.fortran_vec ();
-
-  sigma.resize (nrow_s, ncol_s);
-  double *s_vec = sigma.fortran_vec ();
-
-  if (! (jobv == 'N' || jobv == 'O'))
-    right_sm.resize (nrow_vt, n);
-
-  double *vt = right_sm.fortran_vec ();
-
-  // Query DGESVD for the correct dimension of WORK.
-
-  octave_idx_type lwork = -1;
-
-  Array<double> work (dim_vector (1, 1));
-
-  octave_idx_type one = 1;
-  octave_idx_type m1 = std::max (m, one);
-  octave_idx_type nrow_vt1 = std::max (nrow_vt, one);
-
-  if (svd_driver == svd::GESVD)
-    {
-      F77_XFCN (dgesvd, DGESVD, (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)));
-
-      lwork = static_cast<octave_idx_type> (work(0));
-      work.resize (dim_vector (lwork, 1));
-
-      F77_XFCN (dgesvd, DGESVD, (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)));
-
-    }
-  else if (svd_driver == svd::GESDD)
-    {
-      assert (jobu == jobv);
-      char jobz = jobu;
-      OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn);
-
-      F77_XFCN (dgesdd, DGESDD, (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)));
-
-      lwork = static_cast<octave_idx_type> (work(0));
-      work.resize (dim_vector (lwork, 1));
-
-      F77_XFCN (dgesdd, DGESDD, (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)));
-
-    }
-  else
-    abort ();
-
-  if (! (jobv == 'N' || jobv == 'O'))
-    right_sm = right_sm.transpose ();
-
-  return info;
-}
-
-template <>
-octave_idx_type
-svd<FloatMatrix>::init (const FloatMatrix& a, svd::type svd_type,
-                        svd::driver svd_driver)
-{
-  octave_idx_type info;
-
-  octave_idx_type m = a.rows ();
-  octave_idx_type n = a.cols ();
-
-  if (m == 0 || n == 0)
-    return empty_init (m, n, svd_type);
-
-  FloatMatrix atmp = a;
-  float *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 (svd_type)
-    {
-    case svd::economy:
-      jobu = jobv = 'S';
-      ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
-      break;
-
-    case svd::sigma_only:
+    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
@@ -359,356 +442,46 @@
       break;
     }
 
-  type_computed = svd_type;
-
   if (! (jobu == 'N' || jobu == 'O'))
     left_sm.resize (m, ncol_u);
 
-  float *u = left_sm.fortran_vec ();
+  P* u = left_sm.fortran_vec ();
 
   sigma.resize (nrow_s, ncol_s);
-  float *s_vec = sigma.fortran_vec ();
+  DM_P* s_vec = sigma.fortran_vec ();
 
   if (! (jobv == 'N' || jobv == 'O'))
     right_sm.resize (nrow_vt, n);
 
-  float *vt = right_sm.fortran_vec ();
+  P* vt = right_sm.fortran_vec ();
 
-  // Query SGESVD for the correct dimension of WORK.
+  // Query _GESVD for the correct dimension of WORK.
 
   octave_idx_type lwork = -1;
 
-  Array<float> work (dim_vector (1, 1));
+  T work (1, 1);
 
-  octave_idx_type one = 1;
-  octave_idx_type m1 = std::max (m, one);
-  octave_idx_type nrow_vt1 = std::max (nrow_vt, one);
+  octave_idx_type m1 = std::max (m, 1);
+  octave_idx_type nrow_vt1 = std::max (nrow_vt, 1);
 
-  if (svd_driver == svd::GESVD)
-    {
-      F77_XFCN (sgesvd, SGESVD, (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)));
-
-      lwork = static_cast<octave_idx_type> (work(0));
-      work.resize (dim_vector (lwork, 1));
-
-      F77_XFCN (sgesvd, SGESVD, (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)));
-
-    }
-  else if (svd_driver == svd::GESDD)
+  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*min_mn);
 
-      F77_XFCN (sgesdd, SGESDD, (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)));
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8 * std::min (m, n));
 
-      lwork = static_cast<octave_idx_type> (work(0));
-      work.resize (dim_vector (lwork, 1));
-
-      F77_XFCN (sgesdd, SGESDD, (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)));
-
+      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 ();
-
-  return info;
-}
-
-template <>
-octave_idx_type
-svd<ComplexMatrix>::init (const ComplexMatrix& a, svd::type svd_type,
-                          svd::driver svd_driver)
-{
-  octave_idx_type info;
-
-  octave_idx_type m = a.rows ();
-  octave_idx_type n = a.cols ();
-
-  if (m == 0 || n == 0)
-    return empty_init (m, n, svd_type);
-
-  ComplexMatrix atmp = a;
-  Complex *tmp_data = atmp.fortran_vec ();
-
-  octave_idx_type min_mn = m < n ? m : n;
-  octave_idx_type max_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 (svd_type)
-    {
-    case svd::economy:
-      jobu = jobv = 'S';
-      ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
-      break;
-
-    case svd::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;
-    }
-
-  type_computed = svd_type;
-
-  if (! (jobu == 'N' || jobu == 'O'))
-    left_sm.resize (m, ncol_u);
-
-  Complex *u = left_sm.fortran_vec ();
-
-  sigma.resize (nrow_s, ncol_s);
-  double *s_vec = sigma.fortran_vec ();
-
-  if (! (jobv == 'N' || jobv == 'O'))
-    right_sm.resize (nrow_vt, n);
-
-  Complex *vt = right_sm.fortran_vec ();
-
-  // Query ZGESVD for the correct dimension of WORK.
-
-  octave_idx_type lwork = -1;
-
-  Array<Complex> work (dim_vector (1, 1));
-
-  octave_idx_type one = 1;
-  octave_idx_type m1 = std::max (m, one);
-  octave_idx_type nrow_vt1 = std::max (nrow_vt, one);
-
-  if (svd_driver == svd::GESVD)
-    {
-      octave_idx_type lrwork = 5*max_mn;
-      Array<double> rwork (dim_vector (lrwork, 1));
-
-      F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
-                                 F77_CONST_CHAR_ARG2 (&jobv, 1),
-                                 m, n, F77_DBLE_CMPLX_ARG (tmp_data), m1, s_vec, F77_DBLE_CMPLX_ARG (u), m1, F77_DBLE_CMPLX_ARG (vt),
-                                 nrow_vt1, F77_DBLE_CMPLX_ARG (work.fortran_vec ()), lwork,
-                                 rwork.fortran_vec (), info
-                                 F77_CHAR_ARG_LEN (1)
-                                 F77_CHAR_ARG_LEN (1)));
-
-      lwork = static_cast<octave_idx_type> (work(0).real ());
-      work.resize (dim_vector (lwork, 1));
-
-      F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
-                                 F77_CONST_CHAR_ARG2 (&jobv, 1),
-                                 m, n, F77_DBLE_CMPLX_ARG (tmp_data), m1, s_vec, F77_DBLE_CMPLX_ARG (u), m1, F77_DBLE_CMPLX_ARG (vt),
-                                 nrow_vt1, F77_DBLE_CMPLX_ARG (work.fortran_vec ()), lwork,
-                                 rwork.fortran_vec (), info
-                                 F77_CHAR_ARG_LEN (1)
-                                 F77_CHAR_ARG_LEN (1)));
-    }
-  else if (svd_driver == svd::GESDD)
-    {
-      assert (jobu == jobv);
-      char jobz = jobu;
-
-      octave_idx_type lrwork;
-      if (jobz == 'N')
-        lrwork = 7*min_mn;
-      else
-        lrwork = 5*min_mn*min_mn + 5*min_mn;
-      Array<double> rwork (dim_vector (lrwork, 1));
-
-      OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn);
-
-      F77_XFCN (zgesdd, ZGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
-                                 m, n, F77_DBLE_CMPLX_ARG (tmp_data), m1, s_vec, F77_DBLE_CMPLX_ARG (u), m1, F77_DBLE_CMPLX_ARG (vt),
-                                 nrow_vt1, F77_DBLE_CMPLX_ARG (work.fortran_vec ()), lwork,
-                                 rwork.fortran_vec (), iwork, info
-                                 F77_CHAR_ARG_LEN (1)));
-
-      lwork = static_cast<octave_idx_type> (work(0).real ());
-      work.resize (dim_vector (lwork, 1));
-
-      F77_XFCN (zgesdd, ZGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
-                                 m, n, F77_DBLE_CMPLX_ARG (tmp_data), m1, s_vec, F77_DBLE_CMPLX_ARG (u), m1, F77_DBLE_CMPLX_ARG (vt),
-                                 nrow_vt1, F77_DBLE_CMPLX_ARG (work.fortran_vec ()), lwork,
-                                 rwork.fortran_vec (), iwork, info
-                                 F77_CHAR_ARG_LEN (1)));
-    }
-  else
-    abort ();
-
-  if (! (jobv == 'N' || jobv == 'O'))
-    right_sm = right_sm.hermitian ();
-
-  return info;
-}
-
-template <>
-octave_idx_type
-svd<FloatComplexMatrix>::init (const FloatComplexMatrix& a, svd::type svd_type,
-                              svd::driver svd_driver)
-{
-  octave_idx_type info;
-
-  octave_idx_type m = a.rows ();
-  octave_idx_type n = a.cols ();
-
-  if (m == 0 || n == 0)
-    return empty_init (m, n, svd_type);
-
-  FloatComplexMatrix atmp = a;
-  FloatComplex *tmp_data = atmp.fortran_vec ();
-
-  octave_idx_type min_mn = m < n ? m : n;
-  octave_idx_type max_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 (svd_type)
-    {
-    case svd::economy:
-      jobu = jobv = 'S';
-      ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
-      break;
-
-    case svd::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;
-    }
-
-  type_computed = svd_type;
-
-  if (! (jobu == 'N' || jobu == 'O'))
-    left_sm.resize (m, ncol_u);
-
-  FloatComplex *u = left_sm.fortran_vec ();
-
-  sigma.resize (nrow_s, ncol_s);
-  float *s_vec = sigma.fortran_vec ();
-
-  if (! (jobv == 'N' || jobv == 'O'))
-    right_sm.resize (nrow_vt, n);
-
-  FloatComplex *vt = right_sm.fortran_vec ();
-
-  // Query CGESVD for the correct dimension of WORK.
-
-  octave_idx_type lwork = -1;
-
-  Array<FloatComplex> work (dim_vector (1, 1));
-
-  octave_idx_type one = 1;
-  octave_idx_type m1 = std::max (m, one);
-  octave_idx_type nrow_vt1 = std::max (nrow_vt, one);
-
-  if (svd_driver == svd::GESVD)
-    {
-      octave_idx_type lrwork = 5*max_mn;
-      Array<float> rwork (dim_vector (lrwork, 1));
-
-      F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
-                                 F77_CONST_CHAR_ARG2 (&jobv, 1),
-                                 m, n, F77_CMPLX_ARG (tmp_data), m1, s_vec, F77_CMPLX_ARG (u), m1, F77_CMPLX_ARG (vt),
-                                 nrow_vt1, F77_CMPLX_ARG (work.fortran_vec ()), lwork,
-                                 rwork.fortran_vec (), info
-                                 F77_CHAR_ARG_LEN (1)
-                                 F77_CHAR_ARG_LEN (1)));
-
-      lwork = static_cast<octave_idx_type> (work(0).real ());
-      work.resize (dim_vector (lwork, 1));
-
-      F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
-                                 F77_CONST_CHAR_ARG2 (&jobv, 1),
-                                 m, n, F77_CMPLX_ARG (tmp_data), m1, s_vec, F77_CMPLX_ARG (u), m1, F77_CMPLX_ARG (vt),
-                                 nrow_vt1, F77_CMPLX_ARG (work.fortran_vec ()), lwork,
-                                 rwork.fortran_vec (), info
-                                 F77_CHAR_ARG_LEN (1)
-                                 F77_CHAR_ARG_LEN (1)));
-    }
-  else if (svd_driver == svd::GESDD)
-    {
-      assert (jobu == jobv);
-      char jobz = jobu;
-
-      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);
-      Array<float> rwork (dim_vector (lrwork, 1));
-
-      OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn);
-
-      F77_XFCN (cgesdd, CGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
-                                 m, n, F77_CMPLX_ARG (tmp_data), m1, s_vec, F77_CMPLX_ARG (u), m1, F77_CMPLX_ARG (vt),
-                                 nrow_vt1, F77_CMPLX_ARG (work.fortran_vec ()), lwork,
-                                 rwork.fortran_vec (), iwork, info
-                                 F77_CHAR_ARG_LEN (1)));
-
-      lwork = static_cast<octave_idx_type> (work(0).real ());
-      work.resize (dim_vector (lwork, 1));
-
-      F77_XFCN (cgesdd, CGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
-                                 m, n, F77_CMPLX_ARG (tmp_data), m1, s_vec, F77_CMPLX_ARG (u), m1, F77_CMPLX_ARG (vt),
-                                 nrow_vt1, F77_CMPLX_ARG (work.fortran_vec ()), lwork,
-                                 rwork.fortran_vec (), iwork, info
-                                 F77_CHAR_ARG_LEN (1)));
-    }
-  else
-    abort ();
-
-  if (! (jobv == 'N' || jobv == 'O'))
-    right_sm = right_sm.hermitian ();
-
-  return info;
 }
 
 // Instantiations we need.
diff -r 2f301a30aeed -r 469c817eb256 liboctave/numeric/svd.h
--- a/liboctave/numeric/svd.h	Wed Aug 03 10:23:41 2016 -0700
+++ b/liboctave/numeric/svd.h	Thu Aug 04 20:20:27 2016 +0100
@@ -1,5 +1,6 @@
 /*
 
+Copyright (C) 2016 Carnë Draug
 Copyright (C) 1994-2015 John W. Eaton
 
 This file is part of Octave.
@@ -35,49 +36,40 @@
 
   typedef typename T::real_diag_matrix_type DM_T;
 
-  enum type
+  enum class Type
   {
     std,
     economy,
     sigma_only
   };
 
-  enum driver
+  enum class Driver
   {
     GESVD,
     GESDD
   };
 
   svd (void)
-    : type_computed (), left_sm (), sigma (), right_sm ()
+    : type (), driver (), left_sm (), sigma (), right_sm ()
   { }
 
-  svd (const T& a, type svd_type = svd::std, driver svd_driver = svd::GESVD)
-    : type_computed (), left_sm (), sigma (), right_sm ()
-  {
-    init (a, svd_type, svd_driver);
-  }
-
-  svd (const T& a, octave_idx_type& info, type svd_type = svd::std,
-       driver svd_driver = svd::GESVD)
-    : type_computed (), left_sm (), sigma (), right_sm ()
-  {
-    info = init (a, svd_type, svd_driver);
-  }
+  svd (const T& a, svd::Type type = svd::Type::std,
+       svd::Driver driver = svd::Driver::GESVD);
 
   svd (const svd& a)
-    : type_computed (a.type_computed), left_sm (a.left_sm),
-      sigma (a.sigma), right_sm (a.right_sm)
+    : type (a.type), driver (a.driver), left_sm (a.left_sm), sigma (a.sigma),
+      right_sm (a.right_sm)
   { }
 
   svd& operator = (const svd& a)
   {
     if (this != &a)
       {
-        type_computed = a.type_computed;
+        type = a.type;
         left_sm = a.left_sm;
         sigma = a.sigma;
         right_sm = a.right_sm;
+        driver = a.driver;
       }
 
     return *this;
@@ -93,17 +85,26 @@
 
 private:
 
-  svd::type type_computed;
+  typedef typename T::element_type P;
+  typedef typename DM_T::element_type DM_P;
+
+  svd::Type type;
+  svd::Driver driver;
 
   T left_sm;
   DM_T sigma;
   T right_sm;
 
-  octave_idx_type
-  init (const T& a, type svd_type, driver svd_driver);
+  void gesvd (char& jobu, char& jobv, octave_idx_type m, octave_idx_type n,
+              P* tmp_data, octave_idx_type m1, DM_P* s_vec, P* u, P* vt,
+              octave_idx_type nrow_vt1, T& work, octave_idx_type& lwork,
+              octave_idx_type& info);
 
-  octave_idx_type
-  empty_init (octave_idx_type nr, octave_idx_type nc, type svd_type);
+  void gesdd (char& jobz, octave_idx_type m, octave_idx_type n,
+              P* tmp_data, octave_idx_type m1, DM_P* s_vec, P* u, P* vt,
+              octave_idx_type nrow_vt1, T& work, octave_idx_type& lwork,
+              octave_idx_type* iwork, octave_idx_type& info);
+
 };
 
 #endif