# HG changeset patch
# User John W. Eaton <jwe@octave.org>
# Date 1455651666 18000
# Node ID cbced1c0991678a321a129d236cfd477290604de
# Parent  987c1a79d33fa8ea99df82e9e797c7852afd6de2
better use of templates for svd classes

* liboctave/numeric/svd.h, liboctave/numeric/svd.cc: New files for svd
classes generated from CmplxSVD.cc, CmplxSVD.h, dbleSVD.cc, dbleSVD.h,
fCmplxSVD.cc, fCmplxSVD.h, floatSVD.cc, and floatSVD.h and converted
to templates.

* liboctave/numeric/module.mk: Update.

* __qp__.cc, svd.cc, CMatrix.cc, CMatrix.h, dDiagMatrix.h, dMatrix.cc,
dMatrix.h, fCMatrix.cc, fCMatrix.h, fDiagMatrix.h, fMatrix.cc,
fMatrix.h, oct-norm.cc, mx-defs.h, mx-ext.h: Use new classes.

diff -r 987c1a79d33f -r cbced1c09916 libinterp/corefcn/__qp__.cc
--- a/libinterp/corefcn/__qp__.cc	Tue Feb 16 14:39:52 2016 -0500
+++ b/libinterp/corefcn/__qp__.cc	Tue Feb 16 14:41:06 2016 -0500
@@ -27,7 +27,7 @@
 #include <cfloat>
 
 #include "chol.h"
-#include "dbleSVD.h"
+#include "svd.h"
 #include "mx-m-dm.h"
 #include "EIG.h"
 
@@ -47,7 +47,7 @@
 
   if (! A.is_empty ())
     {
-      SVD A_svd (A);
+      svd<Matrix> A_svd (A);
 
       DiagMatrix S = A_svd.singular_values ();
 
diff -r 987c1a79d33f -r cbced1c09916 libinterp/corefcn/svd.cc
--- a/libinterp/corefcn/svd.cc	Tue Feb 16 14:39:52 2016 -0500
+++ b/libinterp/corefcn/svd.cc	Tue Feb 16 14:41:06 2016 -0500
@@ -24,10 +24,7 @@
 #  include <config.h>
 #endif
 
-#include "CmplxSVD.h"
-#include "dbleSVD.h"
-#include "fCmplxSVD.h"
-#include "floatSVD.h"
+#include "svd.h"
 
 #include "defun.h"
 #include "error.h"
@@ -37,7 +34,23 @@
 #include "utils.h"
 #include "variables.h"
 
-static int Vsvd_driver = SVD::GESVD;
+static std::string Vsvd_driver = "gesvd";
+
+template <typename T>
+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);
+}
+
+template <typename T>
+static typename svd<T>::driver
+svd_driver (void)
+{
+  return Vsvd_driver == "gesvd" ? svd<T>::GESVD : svd<T>::GESDD;
+}
 
 DEFUN (svd, args, nargout,
        "-*- texinfo -*-\n\
@@ -138,142 +151,93 @@
 
   bool isfloat = arg.is_single_type ();
 
-  SVD::type type = ((nargout == 0 || nargout == 1)
-                    ? SVD::sigma_only
-                    : (nargin == 2) ? SVD::economy : SVD::std);
-
-  octave_idx_type nr = arg.rows ();
-  octave_idx_type nc = arg.columns ();
-
-  SVD::driver driver = static_cast<SVD::driver> (Vsvd_driver);
+  if (isfloat)
+    {
+      if (arg.is_real_type ())
+        {
+          FloatMatrix tmp = arg.float_matrix_value ();
 
-  if (nr == 0 || nc == 0)
-    {
-      if (isfloat)
-        {
-          switch (type)
-            {
-            case SVD::std:
-              retval = ovl (FloatDiagMatrix (nr, nr, 1.0f),
-                            FloatMatrix (nr, nc),
-                            FloatDiagMatrix (nc, nc, 1.0f));
-              break;
+          if (tmp.any_element_is_inf_or_nan ())
+            error ("svd: cannot take SVD of matrix containing Inf or NaN values");
+
+          svd<FloatMatrix> result (tmp,
+                                   svd_type<FloatMatrix> (nargin, nargout),
+                                   svd_driver<FloatMatrix> ());
+
+          FloatDiagMatrix sigma = result.singular_values ();
 
-            case SVD::economy:
-              retval = ovl (FloatDiagMatrix (nr, 0, 1.0f),
-                            FloatMatrix (0, 0),
-                            FloatDiagMatrix (0, nc, 1.0f));
-              break;
-
-            case SVD::sigma_only: default:
-              retval(0) = FloatMatrix (0, 1);
-              break;
-            }
+          if (nargout == 0 || nargout == 1)
+            retval(0) = sigma.extract_diag ();
+          else
+            retval = ovl (result.left_singular_matrix (),
+                          sigma,
+                          result.right_singular_matrix ());
         }
-      else
+      else if (arg.is_complex_type ())
         {
-          switch (type)
-            {
-            case SVD::std:
-              retval = ovl (DiagMatrix (nr, nr, 1.0),
-                            Matrix (nr, nc),
-                            DiagMatrix (nc, nc, 1.0));
-              break;
+          FloatComplexMatrix ctmp = arg.float_complex_matrix_value ();
+
+          if (ctmp.any_element_is_inf_or_nan ())
+            error ("svd: cannot take SVD of matrix containing Inf or NaN values");
+
+          svd<FloatComplexMatrix> result (ctmp,
+                                          svd_type<FloatComplexMatrix> (nargin, nargout),
+                                          svd_driver<FloatComplexMatrix> ());
 
-            case SVD::economy:
-              retval = ovl (DiagMatrix (nr, 0, 1.0),
-                            Matrix (0, 0),
-                            DiagMatrix (0, nc, 1.0));
-              break;
+          FloatDiagMatrix sigma = result.singular_values ();
 
-            case SVD::sigma_only: default:
-              retval(0) = Matrix (0, 1);
-              break;
-            }
+          if (nargout == 0 || nargout == 1)
+            retval(0) = sigma.extract_diag ();
+          else
+            retval = ovl (result.left_singular_matrix (),
+                          sigma,
+                          result.right_singular_matrix ());
         }
     }
   else
     {
-      if (isfloat)
+      if (arg.is_real_type ())
         {
-          if (arg.is_real_type ())
-            {
-              FloatMatrix tmp = arg.float_matrix_value ();
+          Matrix tmp = arg.matrix_value ();
 
-              if (tmp.any_element_is_inf_or_nan ())
-                error ("svd: cannot take SVD of matrix containing Inf or NaN values");
+          if (tmp.any_element_is_inf_or_nan ())
+            error ("svd: cannot take SVD of matrix containing Inf or NaN values");
 
-              FloatSVD result (tmp, type, driver);
+          svd<Matrix> result (tmp,
+                              svd_type<Matrix> (nargin, nargout),
+                              svd_driver<Matrix> ());
 
-              FloatDiagMatrix sigma = result.singular_values ();
+          DiagMatrix sigma = result.singular_values ();
 
-              if (nargout == 0 || nargout == 1)
-                retval(0) = sigma.extract_diag ();
-              else
-                retval = ovl (result.left_singular_matrix (),
-                              sigma,
-                              result.right_singular_matrix ());
-            }
-          else if (arg.is_complex_type ())
-            {
-              FloatComplexMatrix ctmp = arg.float_complex_matrix_value ();
+          if (nargout == 0 || nargout == 1)
+            retval(0) = sigma.extract_diag ();
+          else
+            retval = ovl (result.left_singular_matrix (),
+                          sigma,
+                          result.right_singular_matrix ());
+        }
+      else if (arg.is_complex_type ())
+        {
+          ComplexMatrix ctmp = arg.complex_matrix_value ();
 
-              if (ctmp.any_element_is_inf_or_nan ())
-                error ("svd: cannot take SVD of matrix containing Inf or NaN values");
+          if (ctmp.any_element_is_inf_or_nan ())
+            error ("svd: cannot take SVD of matrix containing Inf or NaN values");
 
-              FloatComplexSVD result (ctmp, type, driver);
-
-              FloatDiagMatrix sigma = result.singular_values ();
+          svd<ComplexMatrix> result (ctmp,
+                                     svd_type<ComplexMatrix> (nargin, nargout),
+                                     svd_driver<ComplexMatrix> ());
 
-              if (nargout == 0 || nargout == 1)
-                retval(0) = sigma.extract_diag ();
-              else
-                retval = ovl (result.left_singular_matrix (),
-                              sigma,
-                              result.right_singular_matrix ());
-            }
+          DiagMatrix sigma = result.singular_values ();
+
+          if (nargout == 0 || nargout == 1)
+            retval(0) = sigma.extract_diag ();
+          else
+            retval = ovl (result.left_singular_matrix (),
+                          sigma,
+                          result.right_singular_matrix ());
         }
       else
-        {
-          if (arg.is_real_type ())
-            {
-              Matrix tmp = arg.matrix_value ();
-
-              if (tmp.any_element_is_inf_or_nan ())
-                error ("svd: cannot take SVD of matrix containing Inf or NaN values");
-
-              SVD result (tmp, type, driver);
-
-              DiagMatrix sigma = result.singular_values ();
-
-              if (nargout == 0 || nargout == 1)
-                retval(0) = sigma.extract_diag ();
-              else
-                retval = ovl (result.left_singular_matrix (),
-                              sigma,
-                              result.right_singular_matrix ());
-            }
-          else if (arg.is_complex_type ())
-            {
-              ComplexMatrix ctmp = arg.complex_matrix_value ();
-
-              if (ctmp.any_element_is_inf_or_nan ())
-                error ("svd: cannot take SVD of matrix containing Inf or NaN values");
-
-              ComplexSVD result (ctmp, type, driver);
-
-              DiagMatrix sigma = result.singular_values ();
-
-              if (nargout == 0 || nargout == 1)
-                retval(0) = sigma.extract_diag ();
-              else
-                retval = ovl (result.left_singular_matrix (),
-                              sigma,
-                              result.right_singular_matrix ());
-            }
-          else
-            err_wrong_type_arg ("svd", arg);
-        }
+        err_wrong_type_arg ("svd", arg);
     }
 
   return retval;
diff -r 987c1a79d33f -r cbced1c09916 liboctave/array/CMatrix.cc
--- a/liboctave/array/CMatrix.cc	Tue Feb 16 14:39:52 2016 -0500
+++ b/liboctave/array/CMatrix.cc	Tue Feb 16 14:41:06 2016 -0500
@@ -47,7 +47,7 @@
 #include "CDiagMatrix.h"
 #include "dDiagMatrix.h"
 #include "schur.h"
-#include "CmplxSVD.h"
+#include "svd.h"
 #include "DET.h"
 #include "f77-fcn.h"
 #include "functor.h"
@@ -1137,7 +1137,7 @@
 {
   ComplexMatrix retval;
 
-  ComplexSVD result (*this, SVD::economy);
+  svd<ComplexMatrix> result (*this, svd<ComplexMatrix>::economy);
 
   DiagMatrix S = result.singular_values ();
   ComplexMatrix U = result.left_singular_matrix ();
diff -r 987c1a79d33f -r cbced1c09916 liboctave/array/CMatrix.h
--- a/liboctave/array/CMatrix.h	Tue Feb 16 14:39:52 2016 -0500
+++ b/liboctave/array/CMatrix.h	Tue Feb 16 14:41:06 2016 -0500
@@ -50,6 +50,9 @@
   typedef Matrix real_matrix_type;
   typedef ComplexMatrix complex_matrix_type;
 
+  typedef DiagMatrix real_diag_matrix_type;
+  typedef ComplexDiagMatrix complex_diag_matrix_type;
+
   typedef double real_elt_type;
   typedef Complex complex_elt_type;
 
diff -r 987c1a79d33f -r cbced1c09916 liboctave/array/dDiagMatrix.h
--- a/liboctave/array/dDiagMatrix.h	Tue Feb 16 14:39:52 2016 -0500
+++ b/liboctave/array/dDiagMatrix.h	Tue Feb 16 14:41:06 2016 -0500
@@ -37,9 +37,6 @@
 OCTAVE_API
 DiagMatrix : public MDiagArray2<double>
 {
-  friend class SVD;
-  friend class ComplexSVD;
-
 public:
 
   DiagMatrix (void) : MDiagArray2<double> () { }
diff -r 987c1a79d33f -r cbced1c09916 liboctave/array/dMatrix.cc
--- a/liboctave/array/dMatrix.cc	Tue Feb 16 14:39:52 2016 -0500
+++ b/liboctave/array/dMatrix.cc	Tue Feb 16 14:41:06 2016 -0500
@@ -46,7 +46,7 @@
 #include "PermMatrix.h"
 #include "DET.h"
 #include "schur.h"
-#include "dbleSVD.h"
+#include "svd.h"
 #include "f77-fcn.h"
 #include "functor.h"
 #include "lo-error.h"
@@ -825,7 +825,7 @@
 Matrix
 Matrix::pseudo_inverse (double tol) const
 {
-  SVD result (*this, SVD::economy);
+  svd<Matrix> result (*this, svd<Matrix>::economy);
 
   DiagMatrix S = result.singular_values ();
   Matrix U = result.left_singular_matrix ();
diff -r 987c1a79d33f -r cbced1c09916 liboctave/array/dMatrix.h
--- a/liboctave/array/dMatrix.h	Tue Feb 16 14:39:52 2016 -0500
+++ b/liboctave/array/dMatrix.h	Tue Feb 16 14:41:06 2016 -0500
@@ -49,6 +49,9 @@
   typedef Matrix real_matrix_type;
   typedef ComplexMatrix complex_matrix_type;
 
+  typedef DiagMatrix real_diag_matrix_type;
+  typedef ComplexDiagMatrix complex_diag_matrix_type;
+
   typedef double real_elt_type;
   typedef Complex complex_elt_type;
 
diff -r 987c1a79d33f -r cbced1c09916 liboctave/array/fCMatrix.cc
--- a/liboctave/array/fCMatrix.cc	Tue Feb 16 14:39:52 2016 -0500
+++ b/liboctave/array/fCMatrix.cc	Tue Feb 16 14:41:06 2016 -0500
@@ -47,7 +47,7 @@
 #include "fCColVector.h"
 #include "fCRowVector.h"
 #include "schur.h"
-#include "fCmplxSVD.h"
+#include "svd.h"
 #include "functor.h"
 #include "lo-error.h"
 #include "lo-ieee.h"
@@ -1141,7 +1141,7 @@
 {
   FloatComplexMatrix retval;
 
-  FloatComplexSVD result (*this, SVD::economy);
+  svd<FloatComplexMatrix> result (*this, svd<FloatComplexMatrix>::economy);
 
   FloatDiagMatrix S = result.singular_values ();
   FloatComplexMatrix U = result.left_singular_matrix ();
diff -r 987c1a79d33f -r cbced1c09916 liboctave/array/fCMatrix.h
--- a/liboctave/array/fCMatrix.h	Tue Feb 16 14:39:52 2016 -0500
+++ b/liboctave/array/fCMatrix.h	Tue Feb 16 14:41:06 2016 -0500
@@ -50,6 +50,9 @@
   typedef FloatMatrix real_matrix_type;
   typedef FloatComplexMatrix complex_matrix_type;
 
+  typedef FloatDiagMatrix real_diag_matrix_type;
+  typedef FloatComplexDiagMatrix complex_diag_matrix_type;
+
   typedef float real_elt_type;
   typedef FloatComplex complex_elt_type;
 
diff -r 987c1a79d33f -r cbced1c09916 liboctave/array/fDiagMatrix.h
--- a/liboctave/array/fDiagMatrix.h	Tue Feb 16 14:39:52 2016 -0500
+++ b/liboctave/array/fDiagMatrix.h	Tue Feb 16 14:41:06 2016 -0500
@@ -37,9 +37,6 @@
 OCTAVE_API
 FloatDiagMatrix : public MDiagArray2<float>
 {
-  friend class FloatSVD;
-  friend class FloatComplexSVD;
-
 public:
 
   FloatDiagMatrix (void) : MDiagArray2<float> () { }
diff -r 987c1a79d33f -r cbced1c09916 liboctave/array/fMatrix.cc
--- a/liboctave/array/fMatrix.cc	Tue Feb 16 14:39:52 2016 -0500
+++ b/liboctave/array/fMatrix.cc	Tue Feb 16 14:41:06 2016 -0500
@@ -49,7 +49,7 @@
 #include "f77-fcn.h"
 #include "fMatrix.h"
 #include "schur.h"
-#include "floatSVD.h"
+#include "svd.h"
 #include "functor.h"
 #include "lo-error.h"
 #include "lo-ieee.h"
@@ -832,7 +832,7 @@
 FloatMatrix
 FloatMatrix::pseudo_inverse (float tol) const
 {
-  FloatSVD result (*this, SVD::economy);
+  svd<FloatMatrix> result (*this, svd<FloatMatrix>::economy);
 
   FloatDiagMatrix S = result.singular_values ();
   FloatMatrix U = result.left_singular_matrix ();
diff -r 987c1a79d33f -r cbced1c09916 liboctave/array/fMatrix.h
--- a/liboctave/array/fMatrix.h	Tue Feb 16 14:39:52 2016 -0500
+++ b/liboctave/array/fMatrix.h	Tue Feb 16 14:41:06 2016 -0500
@@ -49,6 +49,9 @@
   typedef FloatMatrix real_matrix_type;
   typedef FloatComplexMatrix complex_matrix_type;
 
+  typedef FloatDiagMatrix real_diag_matrix_type;
+  typedef FloatComplexDiagMatrix complex_diag_matrix_type;
+
   typedef float real_elt_type;
   typedef FloatComplex complex_elt_type;
 
diff -r 987c1a79d33f -r cbced1c09916 liboctave/numeric/CmplxSVD.cc
--- a/liboctave/numeric/CmplxSVD.cc	Tue Feb 16 14:39:52 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,209 +0,0 @@
-/*
-
-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/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#  include <config.h>
-#endif
-
-#include "CmplxSVD.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-#include "oct-locbuf.h"
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (zgesvd, ZGESVD) (F77_CONST_CHAR_ARG_DECL,
-                             F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, const octave_idx_type&,
-                             Complex*, const octave_idx_type&,
-                             double*, Complex*, const octave_idx_type&,
-                             Complex*, const octave_idx_type&, Complex*,
-                             const octave_idx_type&, double*, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL
-                             F77_CHAR_ARG_LEN_DECL);
-
-  F77_RET_T
-  F77_FUNC (zgesdd, ZGESDD) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, const octave_idx_type&,
-                             Complex*, const octave_idx_type&,
-                             double*, Complex*, const octave_idx_type&,
-                             Complex*, const octave_idx_type&, Complex*,
-                             const octave_idx_type&, double*,
-                             octave_idx_type *, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-}
-
-ComplexMatrix
-ComplexSVD::left_singular_matrix (void) const
-{
-  if (type_computed == SVD::sigma_only)
-    (*current_liboctave_error_handler)
-      ("ComplexSVD: U not computed because type == SVD::sigma_only");
-
-  return left_sm;
-}
-
-ComplexMatrix
-ComplexSVD::right_singular_matrix (void) const
-{
-  if (type_computed == SVD::sigma_only)
-    (*current_liboctave_error_handler)
-      ("ComplexSVD: V not computed because type == SVD::sigma_only");
-
-  return right_sm;
-}
-
-octave_idx_type
-ComplexSVD::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 ();
-
-  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, tmp_data, m1, s_vec, u, m1, vt,
-                                 nrow_vt1, 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, tmp_data, m1, s_vec, u, m1, vt,
-                                 nrow_vt1, 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, tmp_data, m1, s_vec, u, m1, vt,
-                                 nrow_vt1, 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, tmp_data, m1, s_vec, u, m1, vt,
-                                 nrow_vt1, work.fortran_vec (), lwork,
-                                 rwork.fortran_vec (), iwork, info
-                                 F77_CHAR_ARG_LEN (1)));
-    }
-  else
-    assert (0); // impossible
-
-  if (! (jobv == 'N' || jobv == 'O'))
-    right_sm = right_sm.hermitian ();
-
-  return info;
-}
diff -r 987c1a79d33f -r cbced1c09916 liboctave/numeric/CmplxSVD.h
--- a/liboctave/numeric/CmplxSVD.h	Tue Feb 16 14:39:52 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,99 +0,0 @@
-/*
-
-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 (octave_CmplxSVD_h)
-#define octave_CmplxSVD_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "dDiagMatrix.h"
-#include "CMatrix.h"
-#include "dbleSVD.h"
-
-class
-OCTAVE_API
-ComplexSVD
-{
-public:
-
-  ComplexSVD (void)
-    : type_computed (), sigma (), left_sm (), right_sm ()
-  { }
-
-  ComplexSVD (const ComplexMatrix& a, SVD::type svd_type = SVD::std,
-              SVD::driver svd_driver = SVD::GESVD)
-    : type_computed (), sigma (), left_sm (), right_sm ()
-  {
-    init (a, svd_type, svd_driver);
-  }
-
-  ComplexSVD (const ComplexMatrix& a, octave_idx_type& info,
-              SVD::type svd_type = SVD::std,
-              SVD::driver svd_driver = SVD::GESVD)
-    : type_computed (), sigma (), left_sm (), right_sm ()
-  {
-    info = init (a, svd_type, svd_driver);
-  }
-
-  ComplexSVD (const ComplexSVD& a)
-    : type_computed (a.type_computed), sigma (a.sigma),
-      left_sm (a.left_sm), right_sm (a.right_sm)
-  { }
-
-  ComplexSVD& operator = (const ComplexSVD& a)
-  {
-    if (this != &a)
-      {
-        type_computed = a.type_computed;
-        sigma = a.sigma;
-        left_sm = a.left_sm;
-        right_sm = a.right_sm;
-      }
-    return *this;
-  }
-
-  ~ComplexSVD (void) { }
-
-  DiagMatrix singular_values (void) const { return sigma; }
-
-  ComplexMatrix left_singular_matrix (void) const;
-
-  ComplexMatrix right_singular_matrix (void) const;
-
-  friend std::ostream&  operator << (std::ostream& os, const ComplexSVD& a);
-
-private:
-
-  SVD::type type_computed;
-
-  DiagMatrix sigma;
-  ComplexMatrix left_sm;
-  ComplexMatrix right_sm;
-
-  octave_idx_type init (const ComplexMatrix& a,
-                        SVD::type svd_type = SVD::std,
-                        SVD::driver svd_driver = SVD::GESVD);
-};
-
-#endif
diff -r 987c1a79d33f -r cbced1c09916 liboctave/numeric/dbleSVD.cc
--- a/liboctave/numeric/dbleSVD.cc	Tue Feb 16 14:39:52 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,205 +0,0 @@
-/*
-
-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/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#  include <config.h>
-#endif
-
-#include <iostream>
-
-#include "dbleSVD.h"
-#include "f77-fcn.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 octave_idx_type&, const octave_idx_type&,
-                             double*, const octave_idx_type&, double*,
-                             double*, const octave_idx_type&, double*,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL
-                             F77_CHAR_ARG_LEN_DECL);
-
-  F77_RET_T
-  F77_FUNC (dgesdd, DGESDD) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, const octave_idx_type&,
-                             double*, const octave_idx_type&, double*,
-                             double*, const octave_idx_type&, double*,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, octave_idx_type *,
-                             octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-}
-
-Matrix
-SVD::left_singular_matrix (void) const
-{
-  if (type_computed == SVD::sigma_only)
-    (*current_liboctave_error_handler)
-      ("SVD: U not computed because type == SVD::sigma_only");
-
-  return left_sm;
-}
-
-Matrix
-SVD::right_singular_matrix (void) const
-{
-  if (type_computed == SVD::sigma_only)
-    (*current_liboctave_error_handler)
-      ("SVD: V not computed because type == SVD::sigma_only");
-
-  return right_sm;
-}
-
-octave_idx_type
-SVD::init (const Matrix& 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 ();
-
-  Matrix atmp = a;
-  double *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:
-
-      // 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
-    assert (0); // impossible
-
-  if (! (jobv == 'N' || jobv == 'O'))
-    right_sm = right_sm.transpose ();
-
-  return info;
-}
-
-std::ostream&
-operator << (std::ostream& os, const SVD& a)
-{
-  os << a.left_singular_matrix () << "\n";
-  os << a.singular_values () << "\n";
-  os << a.right_singular_matrix () << "\n";
-
-  return os;
-}
diff -r 987c1a79d33f -r cbced1c09916 liboctave/numeric/dbleSVD.h
--- a/liboctave/numeric/dbleSVD.h	Tue Feb 16 14:39:52 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,108 +0,0 @@
-/*
-
-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 (octave_dbleSVD_h)
-#define octave_dbleSVD_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "dDiagMatrix.h"
-#include "dMatrix.h"
-
-class
-OCTAVE_API
-SVD
-{
-public:
-
-  enum type
-  {
-    std,
-    economy,
-    sigma_only
-  };
-
-  enum driver
-  {
-    GESVD,
-    GESDD
-  };
-
-  SVD (void) : type_computed (), sigma (), left_sm (), right_sm () { }
-
-  SVD (const Matrix& a,
-       type svd_type = SVD::std, driver svd_driver = SVD::GESVD)
-    : type_computed (), sigma (), left_sm (), right_sm ()
-  {
-    init (a, svd_type, svd_driver);
-  }
-
-  SVD (const Matrix& a, octave_idx_type& info,
-       type svd_type = SVD::std, driver svd_driver = SVD::GESVD)
-    : type_computed (), sigma (), left_sm (), right_sm ()
-  {
-    info = init (a, svd_type, svd_driver);
-  }
-
-  SVD (const SVD& a)
-    : type_computed (a.type_computed), sigma (a.sigma),
-      left_sm (a.left_sm), right_sm (a.right_sm)
-  { }
-
-  SVD& operator = (const SVD& a)
-  {
-    if (this != &a)
-      {
-        type_computed = a.type_computed;
-        sigma = a.sigma;
-        left_sm = a.left_sm;
-        right_sm = a.right_sm;
-      }
-
-    return *this;
-  }
-
-  ~SVD (void) { }
-
-  DiagMatrix singular_values (void) const { return sigma; }
-
-  Matrix left_singular_matrix (void) const;
-
-  Matrix right_singular_matrix (void) const;
-
-  friend std::ostream&  operator << (std::ostream& os, const SVD& a);
-
-private:
-
-  SVD::type type_computed;
-
-  DiagMatrix sigma;
-  Matrix left_sm;
-  Matrix right_sm;
-
-  octave_idx_type init (const Matrix& a,
-                        type svd_type = std, driver svd_driver = GESVD);
-};
-
-#endif
diff -r 987c1a79d33f -r cbced1c09916 liboctave/numeric/fCmplxSVD.cc
--- a/liboctave/numeric/fCmplxSVD.cc	Tue Feb 16 14:39:52 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,216 +0,0 @@
-/*
-
-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/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#  include <config.h>
-#endif
-
-#include "fCmplxSVD.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-#include "oct-locbuf.h"
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (cgesvd, CGESVD) (F77_CONST_CHAR_ARG_DECL,
-                             F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, const octave_idx_type&,
-                             FloatComplex*, const octave_idx_type&, float*,
-                             FloatComplex*, const octave_idx_type&,
-                             FloatComplex*, const octave_idx_type&,
-                             FloatComplex*, const octave_idx_type&,
-                             float*, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL
-                             F77_CHAR_ARG_LEN_DECL);
-
-  F77_RET_T
-  F77_FUNC (cgesdd, CGESDD) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, const octave_idx_type&,
-                             FloatComplex*, const octave_idx_type&, float*,
-                             FloatComplex*, const octave_idx_type&,
-                             FloatComplex*, const octave_idx_type&,
-                             FloatComplex*, const octave_idx_type&,
-                             float*, octave_idx_type *, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-}
-
-FloatComplexMatrix
-FloatComplexSVD::left_singular_matrix (void) const
-{
-  if (type_computed == SVD::sigma_only)
-    {
-      (*current_liboctave_error_handler)
-        ("FloatComplexSVD: U not computed because type == SVD::sigma_only");
-      return FloatComplexMatrix ();
-    }
-  else
-    return left_sm;
-}
-
-FloatComplexMatrix
-FloatComplexSVD::right_singular_matrix (void) const
-{
-  if (type_computed == SVD::sigma_only)
-    {
-      (*current_liboctave_error_handler)
-        ("FloatComplexSVD: V not computed because type == SVD::sigma_only");
-      return FloatComplexMatrix ();
-    }
-  else
-    return right_sm;
-}
-
-octave_idx_type
-FloatComplexSVD::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 ();
-
-  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, tmp_data, m1, s_vec, u, m1, vt,
-                                 nrow_vt1, 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, tmp_data, m1, s_vec, u, m1, vt,
-                                 nrow_vt1, 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, tmp_data, m1, s_vec, u, m1, vt,
-                                 nrow_vt1, 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, tmp_data, m1, s_vec, u, m1, vt,
-                                 nrow_vt1, work.fortran_vec (), lwork,
-                                 rwork.fortran_vec (), iwork, info
-                                 F77_CHAR_ARG_LEN (1)));
-    }
-  else
-    assert (0); // impossible
-
-  if (! (jobv == 'N' || jobv == 'O'))
-    right_sm = right_sm.hermitian ();
-
-  return info;
-}
diff -r 987c1a79d33f -r cbced1c09916 liboctave/numeric/fCmplxSVD.h
--- a/liboctave/numeric/fCmplxSVD.h	Tue Feb 16 14:39:52 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,101 +0,0 @@
-/*
-
-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 (octave_fCmplxSVD_h)
-#define octave_fCmplxSVD_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "fDiagMatrix.h"
-#include "fCMatrix.h"
-#include "dbleSVD.h"
-
-class
-OCTAVE_API
-FloatComplexSVD
-{
-public:
-
-  FloatComplexSVD (void)
-    : type_computed (), sigma (), left_sm (), right_sm ()
-  { }
-
-  FloatComplexSVD (const FloatComplexMatrix& a,
-                   SVD::type svd_type = SVD::std,
-                   SVD::driver svd_driver = SVD::GESVD)
-    : type_computed (), sigma (), left_sm (), right_sm ()
-  {
-    init (a, svd_type, svd_driver);
-  }
-
-  FloatComplexSVD (const FloatComplexMatrix& a, octave_idx_type& info,
-                   SVD::type svd_type = SVD::std,
-                   SVD::driver svd_driver = SVD::GESVD)
-    : type_computed (), sigma (), left_sm (), right_sm ()
-  {
-    info = init (a, svd_type, svd_driver);
-  }
-
-  FloatComplexSVD (const FloatComplexSVD& a)
-    : type_computed (a.type_computed), sigma (a.sigma),
-      left_sm (a.left_sm), right_sm (a.right_sm)
-  { }
-
-  FloatComplexSVD& operator = (const FloatComplexSVD& a)
-  {
-    if (this != &a)
-      {
-        type_computed = a.type_computed;
-        sigma = a.sigma;
-        left_sm = a.left_sm;
-        right_sm = a.right_sm;
-      }
-    return *this;
-  }
-
-  ~FloatComplexSVD (void) { }
-
-  FloatDiagMatrix singular_values (void) const { return sigma; }
-
-  FloatComplexMatrix left_singular_matrix (void) const;
-
-  FloatComplexMatrix right_singular_matrix (void) const;
-
-  friend std::ostream&  operator << (std::ostream& os,
-                                     const FloatComplexSVD& a);
-
-private:
-
-  SVD::type type_computed;
-
-  FloatDiagMatrix sigma;
-  FloatComplexMatrix left_sm;
-  FloatComplexMatrix right_sm;
-
-  octave_idx_type init (const FloatComplexMatrix& a,
-                        SVD::type svd_type = SVD::std,
-                        SVD::driver svd_driver = SVD::GESVD);
-};
-
-#endif
diff -r 987c1a79d33f -r cbced1c09916 liboctave/numeric/floatSVD.cc
--- a/liboctave/numeric/floatSVD.cc	Tue Feb 16 14:39:52 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,206 +0,0 @@
-/*
-
-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/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#  include <config.h>
-#endif
-
-#include <iostream>
-
-#include "floatSVD.h"
-#include "f77-fcn.h"
-#include "oct-locbuf.h"
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (sgesvd, SGESVD) (F77_CONST_CHAR_ARG_DECL,
-                             F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, const octave_idx_type&,
-                             float*, const octave_idx_type&, float*,
-                             float*, const octave_idx_type&, float*,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL
-                             F77_CHAR_ARG_LEN_DECL);
-
-  F77_RET_T
-  F77_FUNC (sgesdd, SGESDD) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, const octave_idx_type&,
-                             float*, const octave_idx_type&, float*,
-                             float*, const octave_idx_type&, float*,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, octave_idx_type *,
-                             octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-}
-
-FloatMatrix
-FloatSVD::left_singular_matrix (void) const
-{
-  if (type_computed == SVD::sigma_only)
-    (*current_liboctave_error_handler)
-      ("FloatSVD: U not computed because type == SVD::sigma_only");
-
-  return left_sm;
-}
-
-FloatMatrix
-FloatSVD::right_singular_matrix (void) const
-{
-  if (type_computed == SVD::sigma_only)
-    (*current_liboctave_error_handler)
-      ("FloatSVD: V not computed because type == SVD::sigma_only");
-
-  return right_sm;
-}
-
-octave_idx_type
-FloatSVD::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 ();
-
-  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:
-
-      // 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);
-
-  float *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);
-
-  float *vt = right_sm.fortran_vec ();
-
-  // Query SGESVD for the correct dimension of WORK.
-
-  octave_idx_type lwork = -1;
-
-  Array<float> 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 (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)
-    {
-      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)));
-
-      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)));
-
-    }
-  else
-    assert (0); // impossible
-
-  if (! (jobv == 'N' || jobv == 'O'))
-    right_sm = right_sm.transpose ();
-
-  return info;
-}
-
-std::ostream&
-operator << (std::ostream& os, const FloatSVD& a)
-{
-  os << a.left_singular_matrix () << "\n";
-  os << a.singular_values () << "\n";
-  os << a.right_singular_matrix () << "\n";
-
-  return os;
-}
diff -r 987c1a79d33f -r cbced1c09916 liboctave/numeric/floatSVD.h
--- a/liboctave/numeric/floatSVD.h	Tue Feb 16 14:39:52 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,98 +0,0 @@
-/*
-
-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 (octave_floatSVD_h)
-#define octave_floatSVD_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "fDiagMatrix.h"
-#include "fMatrix.h"
-#include "dbleSVD.h"
-
-class
-OCTAVE_API
-FloatSVD
-{
-public:
-
-  FloatSVD (void) : type_computed (), sigma (), left_sm (), right_sm () { }
-
-  FloatSVD (const FloatMatrix& a,
-            SVD::type svd_type = SVD::std, SVD::driver svd_driver = SVD::GESVD)
-    : type_computed (), sigma (), left_sm (), right_sm ()
-  {
-    init (a, svd_type, svd_driver);
-  }
-
-  FloatSVD (const FloatMatrix& a, octave_idx_type& info,
-            SVD::type svd_type = SVD::std,
-            SVD::driver svd_driver = SVD::GESVD)
-    : type_computed (), sigma (), left_sm (), right_sm ()
-  {
-    info = init (a, svd_type, svd_driver);
-  }
-
-  FloatSVD (const FloatSVD& a)
-    : type_computed (a.type_computed), sigma (a.sigma),
-      left_sm (a.left_sm), right_sm (a.right_sm)
-  { }
-
-  FloatSVD& operator = (const FloatSVD& a)
-  {
-    if (this != &a)
-      {
-        type_computed = a.type_computed;
-        sigma = a.sigma;
-        left_sm = a.left_sm;
-        right_sm = a.right_sm;
-      }
-
-    return *this;
-  }
-
-  ~FloatSVD (void) { }
-
-  FloatDiagMatrix singular_values (void) const { return sigma; }
-
-  FloatMatrix left_singular_matrix (void) const;
-
-  FloatMatrix right_singular_matrix (void) const;
-
-  friend std::ostream&  operator << (std::ostream& os, const FloatSVD& a);
-
-private:
-
-  SVD::type type_computed;
-
-  FloatDiagMatrix sigma;
-  FloatMatrix left_sm;
-  FloatMatrix right_sm;
-
-  octave_idx_type init (const FloatMatrix& a,
-                        SVD::type svd_type = SVD::std,
-                        SVD::driver svd_driver = SVD::GESVD);
-};
-
-#endif
diff -r 987c1a79d33f -r cbced1c09916 liboctave/numeric/module.mk
--- a/liboctave/numeric/module.mk	Tue Feb 16 14:39:52 2016 -0500
+++ b/liboctave/numeric/module.mk	Tue Feb 16 14:41:06 2016 -0500
@@ -10,7 +10,6 @@
 NUMERIC_INC = \
   liboctave/numeric/CmplxQR.h \
   liboctave/numeric/CmplxQRP.h \
-  liboctave/numeric/CmplxSVD.h \
   liboctave/numeric/CollocWt.h \
   liboctave/numeric/DAE.h \
   liboctave/numeric/DAEFunc.h \
@@ -37,15 +36,12 @@
   liboctave/numeric/chol.h \
   liboctave/numeric/dbleQR.h \
   liboctave/numeric/dbleQRP.h \
-  liboctave/numeric/dbleSVD.h \
   liboctave/numeric/eigs-base.h \
   liboctave/numeric/fCmplxQR.h \
   liboctave/numeric/fCmplxQRP.h \
-  liboctave/numeric/fCmplxSVD.h \
   liboctave/numeric/fEIG.h \
   liboctave/numeric/floatQR.h \
   liboctave/numeric/floatQRP.h \
-  liboctave/numeric/floatSVD.h \
   liboctave/numeric/gepbalance.h \
   liboctave/numeric/hess.h \
   liboctave/numeric/lo-mappers.h \
@@ -63,7 +59,8 @@
   liboctave/numeric/sparse-chol.h \
   liboctave/numeric/sparse-dmsolve.h \
   liboctave/numeric/sparse-lu.h \
-  liboctave/numeric/sparse-qr.h
+  liboctave/numeric/sparse-qr.h \
+  liboctave/numeric/svd.h
 
 NUMERIC_C_SRC = \
   liboctave/numeric/randgamma.c \
@@ -73,7 +70,6 @@
 NUMERIC_SRC = \
   liboctave/numeric/CmplxQR.cc \
   liboctave/numeric/CmplxQRP.cc \
-  liboctave/numeric/CmplxSVD.cc \
   liboctave/numeric/CollocWt.cc \
   liboctave/numeric/DASPK.cc \
   liboctave/numeric/DASRT.cc \
@@ -86,15 +82,12 @@
   liboctave/numeric/chol.cc \
   liboctave/numeric/dbleQR.cc \
   liboctave/numeric/dbleQRP.cc \
-  liboctave/numeric/dbleSVD.cc \
   liboctave/numeric/eigs-base.cc \
   liboctave/numeric/fCmplxQR.cc \
   liboctave/numeric/fCmplxQRP.cc \
-  liboctave/numeric/fCmplxSVD.cc \
   liboctave/numeric/fEIG.cc \
   liboctave/numeric/floatQR.cc \
   liboctave/numeric/floatQRP.cc \
-  liboctave/numeric/floatSVD.cc \
   liboctave/numeric/gepbalance.cc \
   liboctave/numeric/hess.cc \
   liboctave/numeric/lo-mappers.cc \
@@ -110,6 +103,7 @@
   liboctave/numeric/sparse-dmsolve.cc \
   liboctave/numeric/sparse-lu.cc \
   liboctave/numeric/sparse-qr.cc \
+  liboctave/numeric/svd.cc \
   $(NUMERIC_C_SRC)
 
 LIBOCTAVE_TEMPLATE_SRC += \
diff -r 987c1a79d33f -r cbced1c09916 liboctave/numeric/oct-norm.cc
--- a/liboctave/numeric/oct-norm.cc	Tue Feb 16 14:39:52 2016 -0500
+++ b/liboctave/numeric/oct-norm.cc	Tue Feb 16 14:41:06 2016 -0500
@@ -33,33 +33,32 @@
 #include <iostream>
 #include <vector>
 
-#include "oct-cmplx.h"
+#include "Array-util.h"
+#include "Array.h"
+#include "CColVector.h"
+#include "CMatrix.h"
+#include "CRowVector.h"
+#include "CSparse.h"
+#include "dColVector.h"
+#include "dDiagMatrix.h"
+#include "dMatrix.h"
+#include "dRowVector.h"
+#include "dSparse.h"
+#include "fCColVector.h"
+#include "fCMatrix.h"
+#include "fCRowVector.h"
+#include "fColVector.h"
+#include "fDiagMatrix.h"
+#include "fMatrix.h"
+#include "fRowVector.h"
 #include "lo-error.h"
 #include "lo-ieee.h"
 #include "mx-cm-s.h"
-#include "mx-s-cm.h"
 #include "mx-fcm-fs.h"
 #include "mx-fs-fcm.h"
-#include "Array.h"
-#include "Array-util.h"
-#include "CMatrix.h"
-#include "dMatrix.h"
-#include "fCMatrix.h"
-#include "fMatrix.h"
-#include "CColVector.h"
-#include "dColVector.h"
-#include "CRowVector.h"
-#include "dRowVector.h"
-#include "fCColVector.h"
-#include "fColVector.h"
-#include "fCRowVector.h"
-#include "fRowVector.h"
-#include "CSparse.h"
-#include "dSparse.h"
-#include "dbleSVD.h"
-#include "CmplxSVD.h"
-#include "floatSVD.h"
-#include "fCmplxSVD.h"
+#include "mx-s-cm.h"
+#include "oct-cmplx.h"
+#include "svd.h"
 
 // Theory: norm accumulator is an object that has an accum method able
 // to handle both real and complex element, and a cast operator
@@ -475,14 +474,14 @@
 static int max_norm_iter = 100;
 
 // version with SVD for dense matrices
-template <typename MatrixT, typename VectorT, typename SVDT, typename R>
-R matrix_norm (const MatrixT& m, R p, VectorT, SVDT)
+template <typename MatrixT, typename VectorT, typename R>
+R svd_matrix_norm (const MatrixT& m, R p, VectorT)
 {
   R res = 0;
   if (p == 2)
     {
-      SVDT svd (m, SVD::sigma_only);
-      res = svd.singular_values () (0,0);
+      svd<MatrixT> fact (m, svd<MatrixT>::sigma_only);
+      res = fact.singular_values () (0,0);
     }
   else if (p == 1)
     res = xcolnorms (m, 1).max ();
@@ -529,7 +528,7 @@
   OCTAVE_API RTYPE xnorm (const PREFIX##RowVector& x, RTYPE p) \
   { return vector_norm (x, p); } \
   OCTAVE_API RTYPE xnorm (const PREFIX##Matrix& x, RTYPE p) \
-  { return matrix_norm (x, p, PREFIX##Matrix (), PREFIX##SVD ()); } \
+  { return svd_matrix_norm (x, p, PREFIX##Matrix ()); } \
   OCTAVE_API RTYPE xfrobnorm (const PREFIX##Matrix& x) \
   { return vector_norm (x, static_cast<RTYPE> (2)); }
 
diff -r 987c1a79d33f -r cbced1c09916 liboctave/numeric/svd.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/svd.cc	Tue Feb 16 14:41:06 2016 -0500
@@ -0,0 +1,722 @@
+/*
+
+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/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#  include <config.h>
+#endif
+
+#include <cassert>
+
+#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"
+#include "svd.h"
+
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC (dgesvd, DGESVD) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             double*, const octave_idx_type&, double*,
+                             double*, const octave_idx_type&, double*,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dgesdd, DGESDD) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             double*, const octave_idx_type&, double*,
+                             double*, const octave_idx_type&, double*,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, octave_idx_type *,
+                             octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (sgesvd, SGESVD) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             float*, const octave_idx_type&, float*,
+                             float*, const octave_idx_type&, float*,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (sgesdd, SGESDD) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             float*, const octave_idx_type&, float*,
+                             float*, const octave_idx_type&, float*,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, octave_idx_type *,
+                             octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (zgesvd, ZGESVD) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             Complex*, const octave_idx_type&,
+                             double*, Complex*, const octave_idx_type&,
+                             Complex*, const octave_idx_type&, Complex*,
+                             const octave_idx_type&, double*, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (zgesdd, ZGESDD) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             Complex*, const octave_idx_type&,
+                             double*, Complex*, const octave_idx_type&,
+                             Complex*, const octave_idx_type&, Complex*,
+                             const octave_idx_type&, double*,
+                             octave_idx_type *, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+  F77_RET_T
+  F77_FUNC (cgesvd, CGESVD) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             FloatComplex*, const octave_idx_type&, float*,
+                             FloatComplex*, const octave_idx_type&,
+                             FloatComplex*, const octave_idx_type&,
+                             FloatComplex*, const octave_idx_type&,
+                             float*, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (cgesdd, CGESDD) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             FloatComplex*, const octave_idx_type&, float*,
+                             FloatComplex*, const octave_idx_type&,
+                             FloatComplex*, const octave_idx_type&,
+                             FloatComplex*, const octave_idx_type&,
+                             float*, octave_idx_type *, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+}
+
+template <typename T>
+T
+svd<T>::left_singular_matrix (void) const
+{
+  if (type_computed == svd::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_computed == svd::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)
+{
+  assert (nr == 0 || nc == 0);
+
+  static typename T::element_type zero (0);
+  static typename T::element_type one (1);
+
+  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;
+
+    case svd::economy:
+      left_sm = T (nr, 0, zero);
+      sigma = DM_T (0, 0);
+      right_sm = T (0, nc, zero);
+      break;
+
+    case svd::sigma_only:
+    default:
+      sigma = DM_T (0, 1);
+      break;
+    }
+
+  return 0;
+}
+
+// Specializations.
+
+template <>
+octave_idx_type
+svd<Matrix>::init (const Matrix& a, svd::type svd_type, svd::driver svd_driver)
+{
+  octave_idx_type info = 0;
+
+  octave_idx_type m = a.rows ();
+  octave_idx_type n = a.cols ();
+
+  if (m == 0 || n == 0)
+    return empty_init (m, n, svd_type);
+
+  Matrix atmp = a;
+  double *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:
+
+      // 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:
+
+      // 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);
+
+  float *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);
+
+  float *vt = right_sm.fortran_vec ();
+
+  // Query SGESVD for the correct dimension of WORK.
+
+  octave_idx_type lwork = -1;
+
+  Array<float> 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 (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)
+    {
+      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)));
+
+      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)));
+
+    }
+  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, tmp_data, m1, s_vec, u, m1, vt,
+                                 nrow_vt1, 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, tmp_data, m1, s_vec, u, m1, vt,
+                                 nrow_vt1, 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, tmp_data, m1, s_vec, u, m1, vt,
+                                 nrow_vt1, 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, tmp_data, m1, s_vec, u, m1, vt,
+                                 nrow_vt1, 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, tmp_data, m1, s_vec, u, m1, vt,
+                                 nrow_vt1, 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, tmp_data, m1, s_vec, u, m1, vt,
+                                 nrow_vt1, 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, tmp_data, m1, s_vec, u, m1, vt,
+                                 nrow_vt1, 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, tmp_data, m1, s_vec, u, m1, vt,
+                                 nrow_vt1, 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.
+
+template class svd<Matrix>;
+
+template class svd<FloatMatrix>;
+
+template class svd<ComplexMatrix>;
+
+template class svd<FloatComplexMatrix>;
diff -r 987c1a79d33f -r cbced1c09916 liboctave/numeric/svd.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/svd.h	Tue Feb 16 14:41:06 2016 -0500
@@ -0,0 +1,109 @@
+/*
+
+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 (octave_svd_h)
+#define octave_svd_h 1
+
+#include "octave-config.h"
+
+#include <iosfwd>
+
+template <typename T>
+class
+svd
+{
+public:
+
+  typedef typename T::real_diag_matrix_type DM_T;
+
+  enum type
+  {
+    std,
+    economy,
+    sigma_only
+  };
+
+  enum driver
+  {
+    GESVD,
+    GESDD
+  };
+
+  svd (void)
+    : type_computed (), 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 svd& a)
+    : type_computed (a.type_computed), 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;
+        left_sm = a.left_sm;
+        sigma = a.sigma;
+        right_sm = a.right_sm;
+      }
+
+    return *this;
+  }
+
+  ~svd (void) { }
+
+  T left_singular_matrix (void) const;
+
+  DM_T singular_values (void) const { return sigma; }
+
+  T right_singular_matrix (void) const;
+
+private:
+
+  svd::type type_computed;
+
+  T left_sm;
+  DM_T sigma;
+  T right_sm;
+
+  octave_idx_type
+  init (const T& a, type svd_type, driver svd_driver);
+
+  octave_idx_type
+  empty_init (octave_idx_type nr, octave_idx_type nc, type svd_type);
+};
+
+#endif
diff -r 987c1a79d33f -r cbced1c09916 liboctave/operators/mx-defs.h
--- a/liboctave/operators/mx-defs.h	Tue Feb 16 14:39:52 2016 -0500
+++ b/liboctave/operators/mx-defs.h	Tue Feb 16 14:41:06 2016 -0500
@@ -70,10 +70,7 @@
 
 template <typename T> class schur;
 
-class SVD;
-class ComplexSVD;
-class FloatSVD;
-class FloatComplexSVD;
+template <typename T> class svd;
 
 template <typename T> class lu;
 
diff -r 987c1a79d33f -r cbced1c09916 liboctave/operators/mx-ext.h
--- a/liboctave/operators/mx-ext.h	Tue Feb 16 14:39:52 2016 -0500
+++ b/liboctave/operators/mx-ext.h	Tue Feb 16 14:41:06 2016 -0500
@@ -51,10 +51,7 @@
 
 // Result of a Singular Value Decomposition.
 
-#include "dbleSVD.h"
-#include "CmplxSVD.h"
-#include "floatSVD.h"
-#include "fCmplxSVD.h"
+#include "svd.h"
 
 // Result of an Eigenvalue computation.