# HG changeset patch
# User Markus Mützel <markus.muetzel@gmx.de>
# Date 1514927542 -3600
#      Tue Jan 02 22:12:22 2018 +0100
# Node ID 527564cf3740222f34966386a6d39196008fc209
# Parent  0a567c6489ef7c59b47519aeda63c355056f768d
Use new LAPACK functions in gsvd if available (bug #50463).

* gsvd.cc: Check on run-time whether *ggsvd3 functions are available in used
LAPACK library. Use *ggsvd3 functions if available (using function pointers).

* gsvd.h (ggsvd): Add additional parameter.

* lo-lapack-proto.h: Add declarations for new functions.

diff -r 0a567c6489ef -r 527564cf3740 liboctave/numeric/gsvd.cc
--- a/liboctave/numeric/gsvd.cc	Fri Dec 29 22:25:01 2017 +0100
+++ b/liboctave/numeric/gsvd.cc	Tue Jan 02 22:12:22 2018 +0100
@@ -26,6 +26,10 @@
 #  include <config.h>
 #endif
 
+#if not defined (OCTAVE_USE_WINDOWS_API)
+#  include <dlfcn.h>
+#endif
+
 #include <vector>
 
 #include "CMatrix.h"
@@ -37,6 +41,224 @@
 #include "gsvd.h"
 #include "lo-error.h"
 #include "lo-lapack-proto.h"
+#include "oct-shlib.h"
+
+#if defined (OCTAVE_USE_WINDOWS_API)
+static octave::dynamic_library lib_lapack;
+#endif
+
+static bool have_DGGSVD3;
+static bool gsvd_initialized = false;
+
+void initialize_gsvd (void)
+{
+#if defined (OCTAVE_USE_WINDOWS_API)
+  if (! lib_lapack)
+    {
+      // FIXME: This will only work as long as the filename of "liblapack.dll"
+      // doesn't change.
+      octave::dynamic_library lib ("liblapack");
+      lib_lapack = lib;
+      if (! lib_lapack)
+        {
+          // FIXME: Should we really throw an error if the library was not
+          // found?
+          have_DGGSVD3 = false;
+          (*current_liboctave_error_handler)
+            ("gsvd: Library \"liblapack\" not found.");
+          return;
+        }
+    }
+  have_DGGSVD3 = (lib_lapack.search (STRINGIZE (F77_FUNC (dggsvd3, DGGSVD3)))
+                  != nullptr);
+#else
+  have_DGGSVD3 = (dlsym (RTLD_DEFAULT, STRINGIZE (F77_FUNC (dggsvd3, DGGSVD3)))
+                  != nullptr);
+#endif
+
+  gsvd_initialized = true;
+}
+
+typedef F77_RET_T (*fcn_ptr) ( ... );
+template<class T1>
+struct real_ggsvd_ptr
+{
+  typedef F77_RET_T (*type)
+    (F77_CONST_CHAR_ARG_DECL,   // JOBU
+     F77_CONST_CHAR_ARG_DECL,   // JOBV
+     F77_CONST_CHAR_ARG_DECL,   // JOBQ
+     const F77_INT&,            // M
+     const F77_INT&,            // N
+     const F77_INT&,            // P
+     F77_INT &,                 // K
+     F77_INT &,                 // L
+     T1*,                       // A(LDA,N)
+     const F77_INT&,            // LDA
+     T1*,                       // B(LDB,N)
+     const F77_INT&,            // LDB
+     T1*,                       // ALPHA(N)
+     T1*,                       // BETA(N)
+     T1*,                       // U(LDU,M)
+     const F77_INT&,            // LDU
+     T1*,                       // V(LDV,P)
+     const F77_INT&,            // LDV
+     T1*,                       // Q(LDQ,N)
+     const F77_INT&,            // LDQ
+     T1*,                       // WORK
+     F77_INT*,                  // IWORK(N)
+     F77_INT&                   // INFO
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL);
+};
+
+template<class T1>
+struct real_ggsvd3_ptr
+{
+  typedef F77_RET_T (*type)
+    (F77_CONST_CHAR_ARG_DECL,   // JOBU
+     F77_CONST_CHAR_ARG_DECL,   // JOBV
+     F77_CONST_CHAR_ARG_DECL,   // JOBQ
+     const F77_INT&,            // M
+     const F77_INT&,            // N
+     const F77_INT&,            // P
+     F77_INT &,                 // K
+     F77_INT &,                 // L
+     T1*,                       // A(LDA,N)
+     const F77_INT&,            // LDA
+     T1*,                       // B(LDB,N)
+     const F77_INT&,            // LDB
+     T1*,                       // ALPHA(N)
+     T1*,                       // BETA(N)
+     T1*,                       // U(LDU,M)
+     const F77_INT&,            // LDU
+     T1*,                       // V(LDV,P)
+     const F77_INT&,            // LDV
+     T1*,                       // Q(LDQ,N)
+     const F77_INT&,            // LDQ
+     T1*,                       // WORK
+     const F77_INT&,            // LWORK
+     F77_INT*,                  // IWORK(N)
+     F77_INT&                   // INFO
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL);
+};
+
+template<class T1, class T2>
+struct comp_ggsvd_ptr
+{
+  typedef F77_RET_T (*type)
+    (F77_CONST_CHAR_ARG_DECL,   // JOBU
+     F77_CONST_CHAR_ARG_DECL,   // JOBV
+     F77_CONST_CHAR_ARG_DECL,   // JOBQ
+     const F77_INT&,            // M
+     const F77_INT&,            // N
+     const F77_INT&,            // P
+     F77_INT &,                 // K
+     F77_INT &,                 // L
+     T1*,                       // A(LDA,N)
+     const F77_INT&,            // LDA
+     T1*,                       // B(LDB,N)
+     const F77_INT&,            // LDB
+     T2*,                       // ALPHA(N)
+     T2*,                       // BETA(N)
+     T1*,                       // U(LDU,M)
+     const F77_INT&,            // LDU
+     T1*,                       // V(LDV,P)
+     const F77_INT&,            // LDV
+     T1*,                       // Q(LDQ,N)
+     const F77_INT&,            // LDQ
+     T1*,                       // WORK
+     T2*,                       // RWORK
+     F77_INT*,                  // IWORK(N)
+     F77_INT&                   // INFO
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL);
+};
+
+template<class T1, class T2>
+struct comp_ggsvd3_ptr
+{
+  typedef F77_RET_T (*type)
+    (F77_CONST_CHAR_ARG_DECL,   // JOBU
+     F77_CONST_CHAR_ARG_DECL,   // JOBV
+     F77_CONST_CHAR_ARG_DECL,   // JOBQ
+     const F77_INT&,            // M
+     const F77_INT&,            // N
+     const F77_INT&,            // P
+     F77_INT &,                 // K
+     F77_INT &,                 // L
+     T1*,                       // A(LDA,N)
+     const F77_INT&,            // LDA
+     T1*,                       // B(LDB,N)
+     const F77_INT&,            // LDB
+     T2*,                       // ALPHA(N)
+     T2*,                       // BETA(N)
+     T1*,                       // U(LDU,M)
+     const F77_INT&,            // LDU
+     T1*,                       // V(LDV,P)
+     const F77_INT&,            // LDV
+     T1*,                       // Q(LDQ,N)
+     const F77_INT&,            // LDQ
+     T1*,                       // WORK
+     const F77_INT&,            // LWORK
+     T2*,                       // RWORK
+     F77_INT*,                  // IWORK(N)
+     F77_INT&                   // INFO
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL);
+};
+// template specializations
+typedef real_ggsvd3_ptr<F77_DBLE>::type dggsvd3_type;
+typedef real_ggsvd_ptr<F77_DBLE>::type dggsvd_type;
+typedef real_ggsvd3_ptr<F77_REAL>::type sggsvd3_type;
+typedef real_ggsvd_ptr<F77_REAL>::type sggsvd_type;
+typedef comp_ggsvd3_ptr<F77_DBLE_CMPLX, F77_DBLE>::type zggsvd3_type;
+typedef comp_ggsvd_ptr<F77_DBLE_CMPLX, F77_DBLE>::type zggsvd_type;
+typedef comp_ggsvd3_ptr<F77_CMPLX, F77_REAL>::type cggsvd3_type;
+typedef comp_ggsvd_ptr<F77_CMPLX, F77_REAL>::type cggsvd_type;
+
+#if defined (OCTAVE_USE_WINDOWS_API)
+#  define GET_FCN_POINTER(fcn_type, f, F)                                     \
+     fcn_type f_ptr = reinterpret_cast<fcn_type>                              \
+                      (lib_lapack.search(STRINGIZE (F77_FUNC (f, F))));
+#else
+#  define GET_FCN_POINTER(fcn_type, f, F)                                     \
+     fcn_type f_ptr = reinterpret_cast<fcn_type>                              \
+                      (dlsym (RTLD_DEFAULT, STRINGIZE (F77_FUNC (f, F))));
+#endif
+
+// This macro is very similar to F77_XFCN in "f77-fcn.h"
+#define F77_XFCN_GGSVD(fcn_type, f, F, ...)                                   \
+  do                                                                          \
+    {                                                                         \
+      octave_jmp_buf saved_context;                                           \
+      sig_atomic_t saved_octave_interrupt_immediately =                       \
+                                                octave_interrupt_immediately; \
+      f77_exception_encountered = 0;                                          \
+      octave_save_current_context (saved_context);                            \
+      if (octave_set_current_context)                                         \
+        {                                                                     \
+          octave_interrupt_immediately = saved_octave_interrupt_immediately;  \
+          octave_restore_current_context (saved_context);                     \
+          if (f77_exception_encountered)                                      \
+            F77_XFCN_ERROR (f, F);                                            \
+          else                                                                \
+            octave_rethrow_exception ();                                      \
+        }                                                                     \
+      else                                                                    \
+        {                                                                     \
+          GET_FCN_POINTER (fcn_type, f, F)                                    \
+          octave_interrupt_immediately++;                                     \
+          f_ptr (__VA_ARGS__);                                                \
+          octave_interrupt_immediately--;                                     \
+          octave_restore_current_context (saved_context);                     \
+        }                                                                     \
+    }                                                                         \
+  while (0)
 
 namespace octave
 {
@@ -49,20 +271,40 @@
                          double *tmp_dataA, F77_INT m1, double *tmp_dataB,
                          F77_INT p1, Matrix& alpha, Matrix& beta, double *u,
                          F77_INT nrow_u, double *v, F77_INT nrow_v, double *q,
-                         F77_INT nrow_q, Matrix& work, F77_INT *iwork,
-                         F77_INT& info)
+                         F77_INT nrow_q, Matrix& work, F77_INT lwork,
+                         F77_INT *iwork, F77_INT& info)
     {
-      F77_XFCN (dggsvd, DGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
-                                 F77_CONST_CHAR_ARG2 (&jobv, 1),
-                                 F77_CONST_CHAR_ARG2 (&jobq, 1),
-                                 m, n, p, k, l, tmp_dataA, m1,
-                                 tmp_dataB, p1, alpha.fortran_vec (),
-                                 beta.fortran_vec (), u, nrow_u,
-                                 v, nrow_v, q, nrow_q, work.fortran_vec (),
-                                 iwork, info
-                                 F77_CHAR_ARG_LEN (1)
-                                 F77_CHAR_ARG_LEN (1)
-                                 F77_CHAR_ARG_LEN (1)));
+      if (! gsvd_initialized)
+        initialize_gsvd ();
+
+      if (have_DGGSVD3)
+        {
+          F77_XFCN_GGSVD (dggsvd3_type, dggsvd3, DGGSVD3,
+                          F77_CONST_CHAR_ARG2 (&jobu, 1),
+                          F77_CONST_CHAR_ARG2 (&jobv, 1),
+                          F77_CONST_CHAR_ARG2 (&jobq, 1),
+                          m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
+                          alpha.fortran_vec (), beta.fortran_vec (),
+                          u, nrow_u, v, nrow_v, q, nrow_q,
+                          work.fortran_vec (), lwork, iwork, info
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1));
+        }
+      else
+        {
+          F77_XFCN_GGSVD (dggsvd_type, dggsvd, DGGSVD,
+                          F77_CONST_CHAR_ARG2 (&jobu, 1),
+                          F77_CONST_CHAR_ARG2 (&jobv, 1),
+                          F77_CONST_CHAR_ARG2 (&jobq, 1),
+                          m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
+                          alpha.fortran_vec (), beta.fortran_vec (),
+                          u, nrow_u, v, nrow_v, q, nrow_q,
+                          work.fortran_vec (), iwork, info
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1));
+        }
     }
 
     template <>
@@ -74,19 +316,39 @@
                               FloatMatrix& beta, float *u, F77_INT nrow_u,
                               float *v, F77_INT nrow_v, float *q,
                               F77_INT nrow_q, FloatMatrix& work,
-                              F77_INT *iwork, F77_INT& info)
+                              F77_INT lwork, F77_INT *iwork, F77_INT& info)
     {
-      F77_XFCN (sggsvd, SGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
-                                 F77_CONST_CHAR_ARG2 (&jobv, 1),
-                                 F77_CONST_CHAR_ARG2 (&jobq, 1),
-                                 m, n, p, k, l, tmp_dataA, m1,
-                                 tmp_dataB, p1, alpha.fortran_vec (),
-                                 beta.fortran_vec (), u, nrow_u,
-                                 v, nrow_v, q, nrow_q, work.fortran_vec (),
-                                 iwork, info
-                                 F77_CHAR_ARG_LEN (1)
-                                 F77_CHAR_ARG_LEN (1)
-                                 F77_CHAR_ARG_LEN (1)));
+      if (! gsvd_initialized)
+        initialize_gsvd ();
+
+      if (have_DGGSVD3)
+        {
+          F77_XFCN_GGSVD (sggsvd3_type, sggsvd3, SGGSVD3,
+                          F77_CONST_CHAR_ARG2 (&jobu, 1),
+                          F77_CONST_CHAR_ARG2 (&jobv, 1),
+                          F77_CONST_CHAR_ARG2 (&jobq, 1),
+                          m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
+                          alpha.fortran_vec (), beta.fortran_vec (),
+                          u, nrow_u, v, nrow_v, q, nrow_q,
+                          work.fortran_vec (), lwork, iwork, info
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1));
+        }
+      else
+        {
+          F77_XFCN_GGSVD (sggsvd_type, sggsvd, SGGSVD,
+                          F77_CONST_CHAR_ARG2 (&jobu, 1),
+                          F77_CONST_CHAR_ARG2 (&jobv, 1),
+                          F77_CONST_CHAR_ARG2 (&jobq, 1),
+                          m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
+                          alpha.fortran_vec (), beta.fortran_vec (),
+                          u, nrow_u, v, nrow_v, q, nrow_q,
+                          work.fortran_vec (), iwork, info
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1));
+        }
     }
 
     template <>
@@ -98,23 +360,50 @@
                                 Matrix& beta, Complex *u, F77_INT nrow_u,
                                 Complex *v, F77_INT nrow_v, Complex *q,
                                 F77_INT nrow_q, ComplexMatrix& work,
-                                F77_INT *iwork, F77_INT& info)
+                                F77_INT lwork, F77_INT *iwork, F77_INT& info)
     {
+      if (! gsvd_initialized)
+        initialize_gsvd ();
+
       Matrix rwork(2*n, 1);
-      F77_XFCN (zggsvd, ZGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
-                                 F77_CONST_CHAR_ARG2 (&jobv, 1),
-                                 F77_CONST_CHAR_ARG2 (&jobq, 1),
-                                 m, n, p, k, l, F77_DBLE_CMPLX_ARG (tmp_dataA),
-                                 m1, F77_DBLE_CMPLX_ARG (tmp_dataB), p1,
-                                 alpha.fortran_vec (), beta.fortran_vec (),
-                                 F77_DBLE_CMPLX_ARG (u), nrow_u,
-                                 F77_DBLE_CMPLX_ARG (v), nrow_v,
-                                 F77_DBLE_CMPLX_ARG (q), nrow_q,
-                                 F77_DBLE_CMPLX_ARG (work.fortran_vec ()),
-                                 rwork.fortran_vec (), iwork, info
-                                 F77_CHAR_ARG_LEN (1)
-                                 F77_CHAR_ARG_LEN (1)
-                                 F77_CHAR_ARG_LEN (1)));
+      if (have_DGGSVD3)
+        {
+          F77_XFCN_GGSVD (zggsvd3_type, zggsvd3, ZGGSVD3,
+                          F77_CONST_CHAR_ARG2 (&jobu, 1),
+                          F77_CONST_CHAR_ARG2 (&jobv, 1),
+                          F77_CONST_CHAR_ARG2 (&jobq, 1),
+                          m, n, p, k, l,
+                          F77_DBLE_CMPLX_ARG (tmp_dataA), m1,
+                          F77_DBLE_CMPLX_ARG (tmp_dataB), p1,
+                          alpha.fortran_vec (), beta.fortran_vec (),
+                          F77_DBLE_CMPLX_ARG (u), nrow_u,
+                          F77_DBLE_CMPLX_ARG (v), nrow_v,
+                          F77_DBLE_CMPLX_ARG (q), nrow_q,
+                          F77_DBLE_CMPLX_ARG (work.fortran_vec ()),
+                          lwork, rwork.fortran_vec (), iwork, info
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1));
+        }
+      else
+        {
+          F77_XFCN_GGSVD (zggsvd_type, zggsvd, ZGGSVD,
+                          F77_CONST_CHAR_ARG2 (&jobu, 1),
+                          F77_CONST_CHAR_ARG2 (&jobv, 1),
+                          F77_CONST_CHAR_ARG2 (&jobq, 1),
+                          m, n, p, k, l,
+                          F77_DBLE_CMPLX_ARG (tmp_dataA), m1,
+                          F77_DBLE_CMPLX_ARG (tmp_dataB), p1,
+                          alpha.fortran_vec (), beta.fortran_vec (),
+                          F77_DBLE_CMPLX_ARG (u), nrow_u,
+                          F77_DBLE_CMPLX_ARG (v), nrow_v,
+                          F77_DBLE_CMPLX_ARG (q), nrow_q,
+                          F77_DBLE_CMPLX_ARG (work.fortran_vec ()),
+                          rwork.fortran_vec (), iwork, info
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1));
+        }
     }
 
     template <>
@@ -128,24 +417,51 @@
                                      FloatComplex *u, F77_INT nrow_u,
                                      FloatComplex *v, F77_INT nrow_v,
                                      FloatComplex *q, F77_INT nrow_q,
-                                     FloatComplexMatrix& work,
+                                     FloatComplexMatrix& work, F77_INT lwork,
                                      F77_INT *iwork, F77_INT& info)
     {
+      if (! gsvd_initialized)
+        initialize_gsvd ();
+
       FloatMatrix rwork(2*n, 1);
-      F77_XFCN (cggsvd, CGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
-                                 F77_CONST_CHAR_ARG2 (&jobv, 1),
-                                 F77_CONST_CHAR_ARG2 (&jobq, 1),
-                                 m, n, p, k, l, F77_CMPLX_ARG (tmp_dataA), m1,
-                                 F77_CMPLX_ARG (tmp_dataB), p1,
-                                 alpha.fortran_vec (), beta.fortran_vec (),
-                                 F77_CMPLX_ARG (u), nrow_u,
-                                 F77_CMPLX_ARG (v), nrow_v,
-                                 F77_CMPLX_ARG (q), nrow_q,
-                                 F77_CMPLX_ARG (work.fortran_vec ()),
-                                 rwork.fortran_vec (), iwork, info
-                                 F77_CHAR_ARG_LEN (1)
-                                 F77_CHAR_ARG_LEN (1)
-                                 F77_CHAR_ARG_LEN (1)));
+      if (have_DGGSVD3)
+        {
+          F77_XFCN_GGSVD (cggsvd3_type, cggsvd3, CGGSVD3,
+                          F77_CONST_CHAR_ARG2 (&jobu, 1),
+                          F77_CONST_CHAR_ARG2 (&jobv, 1),
+                          F77_CONST_CHAR_ARG2 (&jobq, 1),
+                          m, n, p, k, l,
+                          F77_CMPLX_ARG (tmp_dataA), m1,
+                          F77_CMPLX_ARG (tmp_dataB), p1,
+                          alpha.fortran_vec (), beta.fortran_vec (),
+                          F77_CMPLX_ARG (u), nrow_u,
+                          F77_CMPLX_ARG (v), nrow_v,
+                          F77_CMPLX_ARG (q), nrow_q,
+                          F77_CMPLX_ARG (work.fortran_vec ()), lwork,
+                          rwork.fortran_vec (), iwork, info
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1));
+        }
+      else
+        {
+          F77_XFCN_GGSVD (cggsvd_type, cggsvd, CGGSVD,
+                          F77_CONST_CHAR_ARG2 (&jobu, 1),
+                          F77_CONST_CHAR_ARG2 (&jobv, 1),
+                          F77_CONST_CHAR_ARG2 (&jobq, 1),
+                          m, n, p, k, l,
+                          F77_CMPLX_ARG (tmp_dataA), m1,
+                          F77_CMPLX_ARG (tmp_dataB), p1,
+                          alpha.fortran_vec (), beta.fortran_vec (),
+                          F77_CMPLX_ARG (u), nrow_u,
+                          F77_CMPLX_ARG (v), nrow_v,
+                          F77_CMPLX_ARG (q), nrow_q,
+                          F77_CMPLX_ARG (work.fortran_vec ()),
+                          rwork.fortran_vec (), iwork, info
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1)
+                          F77_CHAR_ARG_LEN (1));
+        }
     }
 
     template <typename T>
@@ -268,19 +584,41 @@
 
       P *q = right_sm.fortran_vec ();
 
-      F77_INT lwork = 3*n;
-      lwork = (lwork > m ? lwork : m);
-      lwork = (lwork > p ? lwork : p) + n;
-
-      T work (lwork, 1);
       real_matrix alpha (n, 1);
       real_matrix beta (n, 1);
 
       std::vector<F77_INT> iwork (n);
 
+      if (! gsvd_initialized)
+        initialize_gsvd ();
+
+      F77_INT lwork;
+      if (have_DGGSVD3)
+        {
+          lwork = -1;
+          T work_tmp (1, 1);
+
+          gsvd<T>::ggsvd (jobu, jobv, jobq, m, n, p, k, l,
+                          tmp_dataA, m, tmp_dataB, p, alpha, beta, u,
+                          nrow_u, v, nrow_v, q, nrow_q, work_tmp, lwork,
+                          iwork.data (), info);
+
+          lwork = static_cast<F77_INT> (std::abs (work_tmp(0, 0)));
+        }
+      else
+        {
+          lwork = 3*n;
+          lwork = (lwork > m ? lwork : m);
+          lwork = (lwork > p ? lwork : p) + n;
+        }
+      info = 0;
+
+      T work (lwork, 1);
+
       gsvd<T>::ggsvd (jobu, jobv, jobq, m, n, p, k, l,
                       tmp_dataA, m, tmp_dataB, p, alpha, beta, u,
-                      nrow_u, v, nrow_v, q, nrow_q, work, iwork.data (), info);
+                      nrow_u, v, nrow_v, q, nrow_q, work, lwork, iwork.data (),
+                      info);
 
       if (f77_exception_encountered)
         (*current_liboctave_error_handler) ("unrecoverable error in *ggsvd");
diff -r 0a567c6489ef -r 527564cf3740 liboctave/numeric/gsvd.h
--- a/liboctave/numeric/gsvd.h	Fri Dec 29 22:25:01 2017 +0100
+++ b/liboctave/numeric/gsvd.h	Tue Jan 02 22:12:22 2018 +0100
@@ -104,7 +104,8 @@
                   P *u, octave_f77_int_type nrow_u,
                   P *v, octave_f77_int_type nrow_v,
                   P *q, octave_f77_int_type nrow_q,
-                  T& work, octave_f77_int_type *iwork,
+                  T& work, octave_f77_int_type lwork,
+                  octave_f77_int_type *iwork,
                   octave_f77_int_type& info);
     };
   }
diff -r 0a567c6489ef -r 527564cf3740 liboctave/numeric/lo-lapack-proto.h
--- a/liboctave/numeric/lo-lapack-proto.h	Fri Dec 29 22:25:01 2017 +0100
+++ b/liboctave/numeric/lo-lapack-proto.h	Tue Jan 02 22:12:22 2018 +0100
@@ -839,7 +839,7 @@
 
   F77_RET_T
   F77_FUNC (sggsvd, SGGSVD)
-   (F77_CONST_CHAR_ARG_DECL,    // JOBU
+    (F77_CONST_CHAR_ARG_DECL,   // JOBU
      F77_CONST_CHAR_ARG_DECL,   // JOBV
      F77_CONST_CHAR_ARG_DECL,   // JOBQ
      const F77_INT&,            // M
@@ -898,7 +898,7 @@
 
   F77_RET_T
   F77_FUNC (cggsvd, CGGSVD)
-   (F77_CONST_CHAR_ARG_DECL,    // JOBU
+    (F77_CONST_CHAR_ARG_DECL,   // JOBU
      F77_CONST_CHAR_ARG_DECL,   // JOBV
      F77_CONST_CHAR_ARG_DECL,   // JOBQ
      const F77_INT&,            // M
@@ -926,6 +926,130 @@
      F77_CHAR_ARG_LEN_DECL
      F77_CHAR_ARG_LEN_DECL);
 
+  // GGSVD3
+
+  F77_RET_T
+  F77_FUNC (dggsvd3, DGGSVD3)
+    (F77_CONST_CHAR_ARG_DECL,   // JOBU
+     F77_CONST_CHAR_ARG_DECL,   // JOBV
+     F77_CONST_CHAR_ARG_DECL,   // JOBQ
+     const F77_INT&,            // M
+     const F77_INT&,            // N
+     const F77_INT&,            // P
+     F77_INT &,                 // K
+     F77_INT &,                 // L
+     F77_DBLE*,                 // A(LDA,N)
+     const F77_INT&,            // LDA
+     F77_DBLE*,                 // B(LDB,N)
+     const F77_INT&,            // LDB
+     F77_DBLE*,                 // ALPHA(N)
+     F77_DBLE*,                 // BETA(N)
+     F77_DBLE*,                 // U(LDU,M)
+     const F77_INT&,            // LDU
+     F77_DBLE*,                 // V(LDV,P)
+     const F77_INT&,            // LDV
+     F77_DBLE*,                 // Q(LDQ,N)
+     const F77_INT&,            // LDQ
+     F77_DBLE*,                 // WORK
+     const F77_INT&,            // LWORK
+     F77_INT*,                  // IWORK(N)
+     F77_INT&                   // INFO
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (sggsvd3, SGGSVD3)
+    (F77_CONST_CHAR_ARG_DECL,   // JOBU
+     F77_CONST_CHAR_ARG_DECL,   // JOBV
+     F77_CONST_CHAR_ARG_DECL,   // JOBQ
+     const F77_INT&,            // M
+     const F77_INT&,            // N
+     const F77_INT&,            // P
+     F77_INT &,                 // K
+     F77_INT &,                 // L
+     F77_REAL*,                 // A
+     const F77_INT&,            // LDA
+     F77_REAL*,                 // B
+     const F77_INT&,            // LDB
+     F77_REAL*,                 // ALPHA
+     F77_REAL*,                 // BETA
+     F77_REAL*,                 // U
+     const F77_INT&,            // LDU
+     F77_REAL*,                 // V
+     const F77_INT&,            // LDV
+     F77_REAL*,                 // Q
+     const F77_INT&,            // LDQ
+     F77_REAL*,                 // WORK
+     const F77_INT&,            // LWORK
+     F77_INT*,                  // IWORK
+     F77_INT&                   // INFO
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (zggsvd3, ZGGSVD3)
+    (F77_CONST_CHAR_ARG_DECL,   // JOBU
+     F77_CONST_CHAR_ARG_DECL,   // JOBV
+     F77_CONST_CHAR_ARG_DECL,   // JOBQ
+     const F77_INT&,            // M
+     const F77_INT&,            // N
+     const F77_INT&,            // P
+     F77_INT &,                 // K
+     F77_INT &,                 // L
+     F77_DBLE_CMPLX*,           // A(LDA,N)
+     const F77_INT&,            // LDA
+     F77_DBLE_CMPLX*,           // B(LDB,N)
+     const F77_INT&,            // LDB
+     F77_DBLE*,                 // ALPHA(N)
+     F77_DBLE*,                 // BETA(N)
+     F77_DBLE_CMPLX*,           // U(LDU,M)
+     const F77_INT&,            // LDU
+     F77_DBLE_CMPLX*,           // V(LDV,P)
+     const F77_INT&,            // LDV
+     F77_DBLE_CMPLX*,           // Q(LDQ,N)
+     const F77_INT&,            // LDQ
+     F77_DBLE_CMPLX*,           // WORK
+     const F77_INT&,            // LWORK
+     F77_DBLE*,                 // RWORK
+     F77_INT*,                  // IWORK(N)
+     F77_INT&                   // INFO
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (cggsvd3, CGGSVD3)
+    (F77_CONST_CHAR_ARG_DECL,   // JOBU
+     F77_CONST_CHAR_ARG_DECL,   // JOBV
+     F77_CONST_CHAR_ARG_DECL,   // JOBQ
+     const F77_INT&,            // M
+     const F77_INT&,            // N
+     const F77_INT&,            // P
+     F77_INT &,                 // K
+     F77_INT &,                 // L
+     F77_CMPLX*,                // A
+     const F77_INT&,            // LDA
+     F77_CMPLX*,                // B
+     const F77_INT&,            // LDB
+     F77_REAL*,                 // ALPHA
+     F77_REAL*,                 // BETA
+     F77_CMPLX*,                // U
+     const F77_INT&,            // LDU
+     F77_CMPLX*,                // V
+     const F77_INT&,            // LDV
+     F77_CMPLX*,                // Q
+     const F77_INT&,            // LDQ
+     F77_CMPLX*,                // WORK
+     const F77_INT&,            // LWORK
+     F77_REAL*,                 // RWORK
+     F77_INT*,                  // IWORK
+     F77_INT&                   // INFO
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL);
+
   // GTSV
 
   F77_RET_T
