--- a/liboctave/numeric/svd.cc	Wed Aug 17 10:37:57 2016 -0400
+++ b/liboctave/numeric/svd.cc	Wed Aug 17 10:55:38 2016 -0400
@@ -42,376 +42,374 @@
 
 namespace octave
 {
-namespace math
-{
-
-template <typename T>
-T
-svd<T>::left_singular_matrix (void) const
-{
-  if (m_type == svd::Type::sigma_only)
-    (*current_liboctave_error_handler)
-      ("svd: U not computed because type == svd::sigma_only");
+  namespace math
+  {
+    template <typename T>
+    T
+    svd<T>::left_singular_matrix (void) const
+    {
+      if (m_type == svd::Type::sigma_only)
+        (*current_liboctave_error_handler)
+          ("svd: U not computed because type == svd::sigma_only");
 
-  return left_sm;
-}
+      return left_sm;
+    }
 
-template <typename T>
-T
-svd<T>::right_singular_matrix (void) const
-{
-  if (m_type == svd::Type::sigma_only)
-    (*current_liboctave_error_handler)
-      ("svd: V not computed because type == svd::sigma_only");
+    template <typename T>
+    T
+    svd<T>::right_singular_matrix (void) const
+    {
+      if (m_type == svd::Type::sigma_only)
+        (*current_liboctave_error_handler)
+          ("svd: V not computed because type == svd::sigma_only");
 
-  return right_sm;
-}
+      return right_sm;
+    }
 
 
-// GESVD specializations
+    // 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.data (), lwork, info \
-                   F77_CHAR_ARG_LEN (1)                       \
-                   F77_CHAR_ARG_LEN (1)))
+#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.data (), 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.data ()),   \
-                   lwork, rwork.data (), 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.data ()),          \
+                     lwork, rwork.data (), 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, std::vector<double>& work,
-                    octave_idx_type& lwork, octave_idx_type& info)
-{
-  GESVD_REAL_STEP (dgesvd, DGESVD);
+    // 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, std::vector<double>& work,
+                        octave_idx_type& lwork, octave_idx_type& info)
+    {
+      GESVD_REAL_STEP (dgesvd, DGESVD);
 
-  lwork = work[0];
-  work.reserve (lwork);
+      lwork = work[0];
+      work.reserve (lwork);
 
-  GESVD_REAL_STEP (dgesvd, DGESVD);
-}
+      GESVD_REAL_STEP (dgesvd, DGESVD);
+    }
 
-// SGESVD
-template<>
-void
-svd<FloatMatrix>::gesvd (char& jobu, char& jobv, octave_idx_type m,
-                         octave_idx_type n, float* tmp_data,
-                         octave_idx_type m1, float* s_vec, float* u, float* vt,
-                         octave_idx_type nrow_vt1, std::vector<float>& work,
-                         octave_idx_type& lwork, octave_idx_type& info)
-{
-  GESVD_REAL_STEP (sgesvd, SGESVD);
+    // 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, std::vector<float>& work,
+                             octave_idx_type& lwork, octave_idx_type& info)
+    {
+      GESVD_REAL_STEP (sgesvd, SGESVD);
 
-  lwork = work[0];
-  work.reserve (lwork);
+      lwork = work[0];
+      work.reserve (lwork);
 
-  GESVD_REAL_STEP (sgesvd, SGESVD);
-}
+      GESVD_REAL_STEP (sgesvd, SGESVD);
+    }
 
-// ZGESVD
-template<>
-void
-svd<ComplexMatrix>::gesvd (char& jobu, char& jobv, octave_idx_type m,
-                           octave_idx_type n, Complex* tmp_data,
-                           octave_idx_type m1, double* s_vec, Complex* u,
-                           Complex* vt, octave_idx_type nrow_vt1,
-                           std::vector<Complex>& work,
-                           octave_idx_type& lwork, octave_idx_type& info)
-{
-  std::vector<double> rwork (5 * std::max (m, n));
+    // 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,
+                               std::vector<Complex>& work,
+                               octave_idx_type& lwork, octave_idx_type& info)
+    {
+      std::vector<double> rwork (5 * std::max (m, n));
 
-  GESVD_COMPLEX_STEP (zgesvd, ZGESVD, F77_DBLE_CMPLX_ARG);
+      GESVD_COMPLEX_STEP (zgesvd, ZGESVD, F77_DBLE_CMPLX_ARG);
 
-  lwork = work[0].real ();
-  work.reserve (lwork);
+      lwork = work[0].real ();
+      work.reserve (lwork);
 
-  GESVD_COMPLEX_STEP (zgesvd, ZGESVD, F77_DBLE_CMPLX_ARG);
-}
+      GESVD_COMPLEX_STEP (zgesvd, ZGESVD, F77_DBLE_CMPLX_ARG);
+    }
 
-// CGESVD
-template<>
-void
-svd<FloatComplexMatrix>::gesvd (char& jobu, char& jobv,
-                                octave_idx_type m, octave_idx_type n,
-                                FloatComplex* tmp_data, octave_idx_type m1,
-                                float* s_vec, FloatComplex* u,
-                                FloatComplex* vt, octave_idx_type nrow_vt1,
-                                std::vector<FloatComplex>& work,
-                                octave_idx_type& lwork, octave_idx_type& info)
-{
-  std::vector<float> rwork (5 * std::max (m, n));
+    // 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,
+                                    std::vector<FloatComplex>& work,
+                                    octave_idx_type& lwork, octave_idx_type& info)
+    {
+      std::vector<float> rwork (5 * std::max (m, n));
 
-  GESVD_COMPLEX_STEP (cgesvd, CGESVD, F77_CMPLX_ARG);
+      GESVD_COMPLEX_STEP (cgesvd, CGESVD, F77_CMPLX_ARG);
 
-  lwork = work[0].real ();
-  work.reserve (lwork);
+      lwork = work[0].real ();
+      work.reserve (lwork);
 
-  GESVD_COMPLEX_STEP (cgesvd, CGESVD, F77_CMPLX_ARG);
-}
+      GESVD_COMPLEX_STEP (cgesvd, CGESVD, F77_CMPLX_ARG);
+    }
 
 #undef GESVD_REAL_STEP
 #undef GESVD_COMPLEX_STEP
 
 
-// GESDD specializations
+    // 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.data (), lwork, iwork, info                 \
-                   F77_CHAR_ARG_LEN (1)))
+#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.data (), 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.data ()), lwork,         \
-                   rwork.data (), 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.data ()), lwork,           \
+                     rwork.data (), 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,
-                    std::vector<double>& work, octave_idx_type& lwork,
-                    octave_idx_type* iwork, octave_idx_type& info)
-{
-  GESDD_REAL_STEP (dgesdd, DGESDD);
+    // 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,
+                        std::vector<double>& work, octave_idx_type& lwork,
+                        octave_idx_type* iwork, octave_idx_type& info)
+    {
+      GESDD_REAL_STEP (dgesdd, DGESDD);
 
-  lwork = work[0];
-  work.reserve (lwork);
+      lwork = work[0];
+      work.reserve (lwork);
 
-  GESDD_REAL_STEP (dgesdd, DGESDD);
-}
+      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,
-                         std::vector<float>& work, octave_idx_type& lwork,
-                         octave_idx_type* iwork, octave_idx_type& info)
-{
-  GESDD_REAL_STEP (sgesdd, SGESDD);
+    // 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,
+                             std::vector<float>& work, octave_idx_type& lwork,
+                             octave_idx_type* iwork, octave_idx_type& info)
+    {
+      GESDD_REAL_STEP (sgesdd, SGESDD);
 
-  lwork = work[0];
-  work.reserve (lwork);
+      lwork = work[0];
+      work.reserve (lwork);
 
-  GESDD_REAL_STEP (sgesdd, SGESDD);
-}
+      GESDD_REAL_STEP (sgesdd, SGESDD);
+    }
 
-// ZGESDD
-template<>
-void
-svd<ComplexMatrix>::gesdd (char& jobz, octave_idx_type m, octave_idx_type n,
-                           Complex* tmp_data, octave_idx_type m1,
-                           double* s_vec, Complex* u,
-                           Complex* vt, octave_idx_type nrow_vt1,
-                           std::vector<Complex>& work, octave_idx_type& lwork,
-                           octave_idx_type* iwork, octave_idx_type& info)
-{
+    // 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,
+                               std::vector<Complex>& 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 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;
+      octave_idx_type lrwork;
+      if (jobz == 'N')
+        lrwork = 7*min_mn;
+      else
+        lrwork = 5*min_mn*min_mn + 5*min_mn;
 
-  std::vector<double> rwork (lrwork);
+      std::vector<double> rwork (lrwork);
 
-  GESDD_COMPLEX_STEP (zgesdd, ZGESDD, F77_DBLE_CMPLX_ARG);
+      GESDD_COMPLEX_STEP (zgesdd, ZGESDD, F77_DBLE_CMPLX_ARG);
 
-  lwork = work[0].real ();
-  work.reserve (lwork);
+      lwork = work[0].real ();
+      work.reserve (lwork);
 
-  GESDD_COMPLEX_STEP (zgesdd, ZGESDD, F77_DBLE_CMPLX_ARG);
-}
+      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,
-                                std::vector<FloatComplex>& 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);
+    // 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,
+                                    std::vector<FloatComplex>& 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);
-  std::vector<float> rwork (lrwork);
+      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);
+      std::vector<float> rwork (lrwork);
 
-  GESDD_COMPLEX_STEP (cgesdd, CGESDD, F77_CMPLX_ARG);
+      GESDD_COMPLEX_STEP (cgesdd, CGESDD, F77_CMPLX_ARG);
 
-  lwork = work[0].real ();
-  work.reserve (lwork);
+      lwork = work[0].real ();
+      work.reserve (lwork);
 
-  GESDD_COMPLEX_STEP (cgesdd, CGESDD, F77_CMPLX_ARG);
-}
+      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)
-  : m_type (type), m_driver (driver), left_sm (), sigma (), right_sm ()
-{
-  octave_idx_type info;
+    template<typename T>
+    svd<T>::svd (const T& a, svd::Type type,
+                 svd::Driver driver)
+      : m_type (type), m_driver (driver), left_sm (), sigma (), right_sm ()
+    {
+      octave_idx_type info;
+
+      octave_idx_type m = a.rows ();
+      octave_idx_type n = a.cols ();
+
+      if (m == 0 || n == 0)
+        {
+          switch (m_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;
 
-  octave_idx_type m = a.rows ();
-  octave_idx_type n = a.cols ();
+            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;
+        }
 
-  if (m == 0 || n == 0)
-    {
+      T atmp = a;
+      P* tmp_data = atmp.fortran_vec ();
+
+      octave_idx_type min_mn = m < n ? m : n;
+
+      char jobu = 'A';
+      char jobv = 'A';
+
+      octave_idx_type ncol_u = m;
+      octave_idx_type nrow_vt = n;
+      octave_idx_type nrow_s = m;
+      octave_idx_type ncol_s = n;
+
       switch (m_type)
         {
-        case svd::Type::std:
-          left_sm = T (m, m, 0);
-          for (octave_idx_type i = 0; i < m; i++)
-            left_sm.xelem (i, i) = 1;
-          sigma = DM_T (m, n);
-          right_sm = T (n, n, 0);
-          for (octave_idx_type i = 0; i < n; i++)
-            right_sm.xelem (i, i) = 1;
-          break;
-
         case svd::Type::economy:
-          left_sm = T (m, 0, 0);
-          sigma = DM_T (0, 0);
-          right_sm = T (0, n, 0);
+          jobu = jobv = 'S';
+          ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
           break;
 
         case svd::Type::sigma_only:
+
+          // Note:  for this case, both jobu and jobv should be 'N', but
+          // there seems to be a bug in dgesvd from Lapack V2.0.  To
+          // demonstrate the bug, set both jobu and jobv to 'N' and find
+          // the singular values of [eye(3), eye(3)].  The result is
+          // [-sqrt(2), -sqrt(2), -sqrt(2)].
+          //
+          // For Lapack 3.0, this problem seems to be fixed.
+
+          jobu = jobv = 'N';
+          ncol_u = nrow_vt = 1;
+          break;
+
         default:
-          sigma = DM_T (0, 1);
           break;
         }
-      return;
-    }
+
+      if (! (jobu == 'N' || jobu == 'O'))
+        left_sm.resize (m, ncol_u);
 
-  T atmp = a;
-  P* tmp_data = atmp.fortran_vec ();
+      P* u = left_sm.fortran_vec ();
 
-  octave_idx_type min_mn = m < n ? m : n;
+      sigma.resize (nrow_s, ncol_s);
+      DM_P* s_vec = sigma.fortran_vec ();
 
-  char jobu = 'A';
-  char jobv = 'A';
+      if (! (jobv == 'N' || jobv == 'O'))
+        right_sm.resize (nrow_vt, n);
+
+      P* vt = right_sm.fortran_vec ();
 
-  octave_idx_type ncol_u = m;
-  octave_idx_type nrow_vt = n;
-  octave_idx_type nrow_s = m;
-  octave_idx_type ncol_s = n;
+      // Query _GESVD for the correct dimension of WORK.
+
+      octave_idx_type lwork = -1;
+
+      std::vector<P> work (1);
 
-  switch (m_type)
-    {
-    case svd::Type::economy:
-      jobu = jobv = 'S';
-      ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
-      break;
-
-    case svd::Type::sigma_only:
+      octave_idx_type m1 = std::max (m, static_cast<octave_idx_type> (1));
+      octave_idx_type nrow_vt1 = std::max (nrow_vt,
+                                           static_cast<octave_idx_type> (1));
 
-      // 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.
+      if (m_driver == svd::Driver::GESVD)
+        gesvd (jobu, jobv, m, n, tmp_data, m1, s_vec, u, vt, nrow_vt1,
+               work, lwork, info);
+      else if (m_driver == svd::Driver::GESDD)
+        {
+          assert (jobu == jobv);
+          char jobz = jobu;
 
-      jobu = jobv = 'N';
-      ncol_u = nrow_vt = 1;
-      break;
+          std::vector<octave_idx_type> iwork (8 * std::min (m, n));
 
-    default:
-      break;
+          gesdd (jobz, m, n, tmp_data, m1, s_vec, u, vt, nrow_vt1,
+                 work, lwork, iwork.data (), info);
+        }
+      else
+        abort ();
+
+      if (! (jobv == 'N' || jobv == 'O'))
+        right_sm = right_sm.transpose ();
     }
 
-  if (! (jobu == 'N' || jobu == 'O'))
-    left_sm.resize (m, ncol_u);
-
-  P* u = left_sm.fortran_vec ();
-
-  sigma.resize (nrow_s, ncol_s);
-  DM_P* s_vec = sigma.fortran_vec ();
-
-  if (! (jobv == 'N' || jobv == 'O'))
-    right_sm.resize (nrow_vt, n);
+    // Instantiations we need.
 
-  P* vt = right_sm.fortran_vec ();
-
-  // Query _GESVD for the correct dimension of WORK.
-
-  octave_idx_type lwork = -1;
+    template class svd<Matrix>;
 
-  std::vector<P> work (1);
-
-  octave_idx_type m1 = std::max (m, static_cast<octave_idx_type> (1));
-  octave_idx_type nrow_vt1 = std::max (nrow_vt,
-    static_cast<octave_idx_type> (1));
+    template class svd<FloatMatrix>;
 
-  if (m_driver == svd::Driver::GESVD)
-    gesvd (jobu, jobv, m, n, tmp_data, m1, s_vec, u, vt, nrow_vt1,
-           work, lwork, info);
-  else if (m_driver == svd::Driver::GESDD)
-    {
-      assert (jobu == jobv);
-      char jobz = jobu;
-
-      std::vector<octave_idx_type> iwork (8 * std::min (m, n));
-
-      gesdd (jobz, m, n, tmp_data, m1, s_vec, u, vt, nrow_vt1,
-             work, lwork, iwork.data (), info);
-    }
-  else
-    abort ();
+    template class svd<ComplexMatrix>;
 
-  if (! (jobv == 'N' || jobv == 'O'))
-    right_sm = right_sm.transpose ();
+    template class svd<FloatComplexMatrix>;
+  }
 }
-
-// Instantiations we need.
-
-template class svd<Matrix>;
-
-template class svd<FloatMatrix>;
-
-template class svd<ComplexMatrix>;
-
-template class svd<FloatComplexMatrix>;
-
-}
-}