# HG changeset patch
# User Markus Mützel <markus.muetzel@gmx.de>
# Date 1515358084 -3600
# Node ID f0de21a6a4268f07668b89647bcac76e8b7b5c95
# Parent  ef60416c4686bba785ee2121d4d0e9050b1ae578
Use new LAPACK functions in gsvd if available (bug #50463).

* oct-shlib.cc: Search all loaded functions if function name is empty.

* oct-shlib.h (destructor): Prevent double free on close.

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

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

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

* configure.ac: Add check for psapi.

diff -r ef60416c4686 -r f0de21a6a426 configure.ac
--- a/configure.ac	Mon Jan 29 16:22:30 2018 -0500
+++ b/configure.ac	Sun Jan 07 21:48:04 2018 +0100
@@ -1978,6 +1978,9 @@
   AC_SUBST(X11_LIBS)
 fi
 
+## FIXME: This check should probably be for Windows only.
+AC_CHECK_LIB([psapi], [EnumProcessModules])
+
 ### Check for the Carbon framework on macOS systems
 
 OCTAVE_HAVE_FRAMEWORK([Carbon],
diff -r ef60416c4686 -r f0de21a6a426 liboctave/numeric/gsvd.cc
--- a/liboctave/numeric/gsvd.cc	Mon Jan 29 16:22:30 2018 -0500
+++ b/liboctave/numeric/gsvd.cc	Sun Jan 07 21:48:04 2018 +0100
@@ -37,6 +37,176 @@
 #include "gsvd.h"
 #include "lo-error.h"
 #include "lo-lapack-proto.h"
+#include "oct-shlib.h"
+
+static octave::dynamic_library dyn_libs;
+
+static bool have_DGGSVD3 = false;
+static bool gsvd_initialized = false;
+
+/* Hack to stringize macro results. */
+#define xSTRINGIZE(x) #x
+#define STRINGIZE(x) xSTRINGIZE(x)
+
+void initialize_gsvd (void)
+{
+  if (! dyn_libs)
+    {
+      dyn_libs = octave::dynamic_library ("");
+      if (! dyn_libs)
+        {
+          // FIXME: Should we throw an error if we cannot check the libraries?
+          have_DGGSVD3 = false;
+          return;
+        }
+    }
+  have_DGGSVD3 = (dyn_libs.search (STRINGIZE (F77_FUNC (dggsvd3, DGGSVD3)))
+                  != nullptr);
+
+  gsvd_initialized = true;
+}
+
+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;
 
 namespace octave
 {
@@ -49,20 +219,42 @@
                          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)
+        {
+          dggsvd3_type f_ptr = reinterpret_cast<dggsvd3_type>
+            (dyn_libs.search (STRINGIZE (F77_FUNC (dggsvd3, DGGSVD3))));
+          f_ptr (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
+        {
+          dggsvd_type f_ptr = reinterpret_cast<dggsvd_type>
+            (dyn_libs.search (STRINGIZE (F77_FUNC (dggsvd, DGGSVD))));
+          f_ptr (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 +266,41 @@
                               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)
+        {
+          sggsvd3_type f_ptr = reinterpret_cast<sggsvd3_type>
+            (dyn_libs.search (STRINGIZE (F77_FUNC (sggsvd3, SGGSVD3))));
+          f_ptr (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
+        {
+          sggsvd_type f_ptr = reinterpret_cast<sggsvd_type>
+            (dyn_libs.search (STRINGIZE (F77_FUNC (sggsvd, SGGSVD))));
+          f_ptr (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 +312,52 @@
                                 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)
+        {
+          zggsvd3_type f_ptr = reinterpret_cast<zggsvd3_type>
+            (dyn_libs.search (STRINGIZE (F77_FUNC (zggsvd3, ZGGSVD3))));
+          f_ptr (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
+        {
+          zggsvd_type f_ptr = reinterpret_cast<zggsvd_type>
+            (dyn_libs.search (STRINGIZE (F77_FUNC (zggsvd, ZGGSVD))));
+          f_ptr (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 +371,53 @@
                                      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)
+        {
+          cggsvd3_type f_ptr = reinterpret_cast<cggsvd3_type>
+            (dyn_libs.search (STRINGIZE (F77_FUNC (cggsvd3, CGGSVD3))));
+          f_ptr (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
+        {
+          cggsvd_type f_ptr = reinterpret_cast<cggsvd_type>
+            (dyn_libs.search (STRINGIZE (F77_FUNC (cggsvd, CGGSVD))));
+          f_ptr (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 +540,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 (info < 0)
         (*current_liboctave_error_handler) ("*ggsvd.f: argument %d illegal",
diff -r ef60416c4686 -r f0de21a6a426 liboctave/numeric/gsvd.h
--- a/liboctave/numeric/gsvd.h	Mon Jan 29 16:22:30 2018 -0500
+++ b/liboctave/numeric/gsvd.h	Sun Jan 07 21:48:04 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 ef60416c4686 -r f0de21a6a426 liboctave/numeric/lo-lapack-proto.h
--- a/liboctave/numeric/lo-lapack-proto.h	Mon Jan 29 16:22:30 2018 -0500
+++ b/liboctave/numeric/lo-lapack-proto.h	Sun Jan 07 21:48:04 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
diff -r ef60416c4686 -r f0de21a6a426 liboctave/util/oct-shlib.cc
--- a/liboctave/util/oct-shlib.cc	Mon Jan 29 16:22:30 2018 -0500
+++ b/liboctave/util/oct-shlib.cc	Sun Jan 07 21:48:04 2018 +0100
@@ -52,6 +52,7 @@
 #elif defined (HAVE_LOADLIBRARY_API)
 #  define WIN32_LEAN_AND_MEAN 1
 #  include <windows.h>
+#  include <psapi.h>
 #endif
 }
 
@@ -186,10 +187,13 @@
   private:
 
     void *library;
+
+    bool search_all_loaded;
   };
 
   octave_dlopen_shlib::octave_dlopen_shlib (const std::string& f)
-    : dynamic_library::dynlib_rep (f), library (nullptr)
+    : dynamic_library::dynlib_rep (f), library (nullptr),
+      search_all_loaded (false)
   {
     int flags = 0;
 
@@ -207,6 +211,12 @@
     flags |= RTLD_GLOBAL;
 #  endif
 
+    if (file.empty())
+      {
+        search_all_loaded = true;
+        return;
+      }
+
     library = dlopen (file.c_str (), flags);
 
     if (! library)
@@ -234,7 +244,7 @@
   {
     void *function = nullptr;
 
-    if (! is_open ())
+    if (! search_all_loaded && ! is_open ())
       (*current_liboctave_error_handler)
         ("shared library %s is not open", file.c_str ());
 
@@ -243,7 +253,10 @@
     if (mangler)
       sym_name = mangler (name);
 
-    function = dlsym (library, sym_name.c_str ());
+    if (search_all_loaded)
+      function = dlsym (RTLD_DEFAULT, sym_name.c_str ());
+    else
+      function = dlsym (library, sym_name.c_str ());
 
     return function;
   }
@@ -273,13 +286,21 @@
   private:
 
     shl_t library;
+
+    bool search_all_loaded;
   };
 
   octave_shl_load_shlib::octave_shl_load_shlib (const std::string& f)
-    : dynamic_library::dynlib_rep (f), library (0)
+    : dynamic_library::dynlib_rep (f), library (0), search_all_loaded (false)
   {
     file = f;
 
+    if (file.empty())
+      {
+        search_all_loaded = true;
+        return;
+      }
+
     library = shl_load (file.c_str (), BIND_IMMEDIATE, 0L);
 
     if (! library)
@@ -301,7 +322,7 @@
   {
     void *function = nullptr;
 
-    if (! is_open ())
+    if (! search_all_loaded && ! is_open ())
       (*current_liboctave_error_handler)
         ("shared library %s is not open", file.c_str ());
 
@@ -310,8 +331,12 @@
     if (mangler)
       sym_name = mangler (name);
 
-    int status = shl_findsym (&library, sym_name.c_str (),
-                              TYPE_UNDEFINED, &function);
+    if (search_all_loaded)
+      int status = shl_findsym (nullptr, sym_name.c_str (),
+                                TYPE_UNDEFINED, &function);
+    else
+      int status = shl_findsym (&library, sym_name.c_str (),
+                                TYPE_UNDEFINED, &function);
 
     return function;
   }
@@ -336,11 +361,15 @@
     void * search (const std::string& name,
                    dynamic_library::name_mangler mangler = 0);
 
+    void * global_search (const std::string& sym_name);
+
     bool is_open (void) const { return (handle != 0); }
 
   private:
 
     HINSTANCE handle;
+
+    bool search_all_loaded;
   };
 
   static void
@@ -350,8 +379,14 @@
   }
 
   octave_w32_shlib::octave_w32_shlib (const std::string& f)
-    : dynamic_library::dynlib_rep (f), handle (0)
+    : dynamic_library::dynlib_rep (f), handle (0), search_all_loaded (false)
   {
+    if (f.empty())
+      {
+        search_all_loaded = true;
+        return;
+      }
+
     std::string dir = sys::file_ops::dirname (f);
 
     set_dll_directory (dir);
@@ -395,12 +430,61 @@
   }
 
   void *
+  octave_w32_shlib::global_search (const std::string& sym_name)
+  {
+    void *function = nullptr;
+
+    HANDLE proc = GetCurrentProcess ();
+
+    if (! proc)
+      (*current_liboctave_error_handler)
+        ("Unable to get handle to own process.");
+
+    size_t lib_num = 64;
+    size_t size_lib = sizeof (HMODULE);
+    HMODULE *h_libs;
+    DWORD bytes_all_libs;
+    bool got_libs;
+
+    // Get a list of all the libraries in own process.
+    h_libs = static_cast<HMODULE *> (malloc (size_lib*lib_num));
+    got_libs = EnumProcessModules (proc, h_libs, size_lib*lib_num,
+                                   &bytes_all_libs);
+    int ii = 0;
+    while (((size_lib*lib_num) < bytes_all_libs) && ii++ < 3)
+      {
+        lib_num = bytes_all_libs / size_lib;
+        h_libs = static_cast<HMODULE *> (realloc (h_libs, bytes_all_libs));
+        got_libs = EnumProcessModules (proc, h_libs, bytes_all_libs,
+                                       &bytes_all_libs);
+      }
+
+     if (got_libs)
+      {
+        for (int i = 0; i < (bytes_all_libs / size_lib); i++)
+          {
+              // Check for function in library.
+              function = reinterpret_cast<void *>
+                           (GetProcAddress (h_libs[i], sym_name.c_str ()));
+
+              if (function)
+                break;
+          }
+      }
+
+    // Release the handle to the process.
+    CloseHandle (proc);
+
+    return function;
+  }
+
+  void *
   octave_w32_shlib::search (const std::string& name,
                             dynamic_library::name_mangler mangler)
   {
     void *function = nullptr;
 
-    if (! is_open ())
+    if (! search_all_loaded && ! is_open ())
       (*current_liboctave_error_handler)
         ("shared library %s is not open", file.c_str ());
 
@@ -409,8 +493,11 @@
     if (mangler)
       sym_name = mangler (name);
 
-    function = reinterpret_cast<void *> (GetProcAddress (handle,
-                                                         sym_name.c_str ()));
+    if (search_all_loaded)
+      function = global_search (sym_name);
+    else
+      function = reinterpret_cast<void *> (GetProcAddress (handle,
+                                                           sym_name.c_str ()));
 
     return function;
   }
@@ -445,11 +532,17 @@
 
     NSObjectFileImage img;
     NSModule handle;
+
+    bool search_all_loaded;
   };
 
   octave_dyld_shlib::octave_dyld_shlib (const std::string& f)
     : dynamic_library::dynlib_rep (f), handle (0)
   {
+    if (f.empty())
+      (*current_liboctave_error_handler)
+        ("global search is not implemented for DYLD_API");
+
     int returnCode = NSCreateObjectFileImageFromFile (file.c_str (), &img);
 
     if (NSObjectFileImageSuccess != returnCode)
diff -r ef60416c4686 -r f0de21a6a426 liboctave/util/oct-shlib.h
--- a/liboctave/util/oct-shlib.h	Mon Jan 29 16:22:30 2018 -0500
+++ b/liboctave/util/oct-shlib.h	Sun Jan 07 21:48:04 2018 +0100
@@ -123,8 +123,7 @@
 
     ~dynamic_library (void)
     {
-      if (--rep->count == 0)
-        delete rep;
+      rep->count--;
     }
 
     dynamic_library (const dynamic_library& sl)