# HG changeset patch
# User Eddy Xiao <bewantbe@gmail.com>
# Date 1550452096 -28800
#      Mon Feb 18 09:08:16 2019 +0800
# Node ID 5438e38e5cb057ed934d6c9fb360ff11c6ecb56f
# Parent  f85d491d232270ccb5fad8cf015d7a1bb5050e7f
Add GEJSV routines as additional "svd drivers" (bug #55727).

* liboctave/numeric/svd.cc
    (class svd): add a member function gejsv;
    (class svd::Driver): add enum GEJSV.
* liboctave/numeric/svd.h
    (class svd): add a member function gejsv;
    (svd<T>::svd): call gejsv if it is selected svd driver;
    (class Get_jsv_lwork): add utilize class for computing workspace size of gejsv.
* liboctave/numeric/lo-lapack-proto.h: Declare LAPACK functions: GEJSV for singular value decomposition, and GELQF, ORMLQ, ORMQR for workspace size query.
* libinterp/corefcn/svd.cc
    (svd_driver): add GEJSV.
    (Fsvd): Document GEJSV; Add a test for GEJSV.
    (Fsvd_driver): Document GEJSV; Add driver name "gejsv"; Add a test for GEJSV.

diff -r f85d491d2322 -r 5438e38e5cb0 libinterp/corefcn/svd.cc
--- a/libinterp/corefcn/svd.cc	Fri Feb 15 08:29:56 2019 -0800
+++ b/libinterp/corefcn/svd.cc	Mon Feb 18 09:08:16 2019 +0800
@@ -59,9 +59,12 @@
 static typename octave::math::svd<T>::Driver
 svd_driver (void)
 {
-  return (Vsvd_driver == "gesvd"
-          ? octave::math::svd<T>::Driver::GESVD
-          : octave::math::svd<T>::Driver::GESDD);
+  if (Vsvd_driver == "gejsv")
+      return octave::math::svd<T>::Driver::GEJSV;
+  else if (Vsvd_driver == "gesdd")
+      return octave::math::svd<T>::Driver::GESDD;
+  else
+      return octave::math::svd<T>::Driver::GESVD;  // default
 }
 
 DEFUN (svd, args, nargout,
@@ -161,7 +164,8 @@
 singular matrices in addition to singular values) there is a choice of two
 routines in @sc{lapack}.  The default routine used by Octave is @code{gesvd}.
 The alternative is @code{gesdd} which is 5X faster, but may use more memory
-and may be inaccurate for some input matrices.  See the documentation for
+and may be inaccurate for some input matrices. In contrast, for best accuracy
+of tiny singular values, try @code{gejsv}. See the documentation for
 @code{svd_driver} for more information on choosing a driver.
 @seealso{svd_driver, svds, eig, lu, chol, hess, qr, qz}
 @end deftypefn */)
@@ -374,6 +378,15 @@
 
 %!assert <*55710> (1 / svd (-0), Inf)
 
+%!test
+%! old_driver = svd_driver ("gejsv");
+%! s0 = [1 1e-10 1e-20];
+%! q = sqrt (0.5);
+%! a = fliplr (s0)' .* [q 0 -q; -0.5 q -0.5; 0.5 q 0.5];
+%! s1 = svd (a);
+%! svd_driver (old_driver);
+%! assert (abs(s1 - s0) ./ s0 < sqrt (eps));
+
 %!error svd ()
 %!error svd ([1, 2; 4, 5], 2, 3)
 */
@@ -385,14 +398,14 @@
 @deftypefnx {} {} svd_driver (@var{new_val}, "local")
 Query or set the underlying @sc{lapack} driver used by @code{svd}.
 
-Currently recognized values are @qcode{"gesdd"} and @qcode{"gesvd"}.
+Currently recognized values are @qcode{"gesdd"}, @qcode{"gesvd"} and @qcode{"gejsv"}.
 The default is @qcode{"gesvd"}.
 
 When called from inside a function with the @qcode{"local"} option, the
 variable is changed locally for the function and any subroutines it calls.
 The original variable value is restored when exiting the function.
 
-Algorithm Notes: The @sc{lapack} library provides two routines for calculating
+Algorithm Notes: The @sc{lapack} library provides three routines for calculating
 the full singular value decomposition (left and right singular matrices as
 well as singular values).  When calculating just the singular values the
 following discussion is not relevant.
@@ -403,6 +416,11 @@
 For an @nospell{MxN} input matrix the memory usage is of order O(min(M,N) ^ 2),
 whereas the alternative is of order O(max(M,N)).
 
+The @code{gejsv} uses a preconditioned Jacobi SVD algorithm. It can sometimes 
+compute tiny singular values and their singular vectors much more accurately 
+than other SVD drivers, at the cost of speed (mostly single thread) and memory 
+(O(min(M,N) ^ 2 + M + N)).
+
 Beyond speed and memory issues, there have been instances where some input
 matrices were not accurately decomposed by @code{gesdd}.  See currently active
 bug @url{https://savannah.gnu.org/bugs/?55564}.  Until these accuracy issues
@@ -412,7 +430,7 @@
 @seealso{svd}
 @end deftypefn */)
 {
-  static const char *driver_names[] = { "gesvd", "gesdd", nullptr };
+  static const char *driver_names[] = { "gesvd", "gesdd", "gejsv", nullptr };
 
   return SET_INTERNAL_VARIABLE_CHOICES (svd_driver, driver_names);
 }
@@ -424,8 +442,14 @@
 %! [U1, S1, V1] = svd (A);
 %! svd_driver ("gesdd");
 %! [U2, S2, V2] = svd (A);
+%! svd_driver ("gejsv");
+%! [U3, S3, V3] = svd (A);
+%! assert(svd_driver (), "gejsv")
 %! svd_driver (old_driver);
 %! assert (U1, U2, 5*eps);
 %! assert (S1, S2, 5*eps);
 %! assert (V1, V2, 5*eps);
+%! assert (U1, U3, 5*eps);
+%! assert (S1, S3, 5*eps);
+%! assert (V1, V3, 5*eps);
 */
diff -r f85d491d2322 -r 5438e38e5cb0 liboctave/numeric/lo-lapack-proto.h
--- a/liboctave/numeric/lo-lapack-proto.h	Fri Feb 15 08:29:56 2019 -0800
+++ b/liboctave/numeric/lo-lapack-proto.h	Mon Feb 18 09:08:16 2019 +0800
@@ -263,6 +263,124 @@
                              F77_DBLE_CMPLX*, F77_DBLE_CMPLX*,
                              const F77_INT&, F77_INT&);
 
+  // GELQF
+
+  F77_RET_T
+  F77_FUNC (cgelqf, CGELQF) (const F77_INT&, const F77_INT&,
+                             F77_CMPLX*, const F77_INT&,
+                             F77_CMPLX*, F77_CMPLX*,
+                             const F77_INT&, F77_INT&);
+
+  F77_RET_T
+  F77_FUNC (dgelqf, DGELQF) (const F77_INT&, const F77_INT&,
+                             F77_DBLE*, const F77_INT&,
+                             F77_DBLE*, F77_DBLE*,
+                             const F77_INT&, F77_INT&);
+
+  F77_RET_T
+  F77_FUNC (sgelqf, SGELQF) (const F77_INT&, const F77_INT&,
+                             F77_REAL*, const F77_INT&,
+                             F77_REAL*, F77_REAL*,
+                             const F77_INT&, F77_INT&);
+
+  F77_RET_T
+  F77_FUNC (zgelqf, ZGELQF) (const F77_INT&, const F77_INT&,
+                             F77_DBLE_CMPLX*, const F77_INT&,
+                             F77_DBLE_CMPLX*, F77_DBLE_CMPLX*,
+                             const F77_INT&, F77_INT&);
+
+  // ORMLQ
+  
+  F77_RET_T
+  F77_FUNC (cormlq, CORMLQ) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const F77_INT&, const F77_INT&, const F77_INT&,
+                             F77_CMPLX*, const F77_INT&,
+                             F77_CMPLX*, F77_CMPLX*,
+                             const F77_INT&, F77_CMPLX*,
+                             const F77_INT&, F77_INT&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dormlq, DORMLQ) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const F77_INT&, const F77_INT&, const F77_INT&,
+                             F77_DBLE*, const F77_INT&,
+                             F77_DBLE*, F77_DBLE*,
+                             const F77_INT&, F77_DBLE*,
+                             const F77_INT&, F77_INT&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (sormlq, SORMLQ) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const F77_INT&, const F77_INT&, const F77_INT&,
+                             F77_REAL*, const F77_INT&,
+                             F77_REAL*, F77_REAL*,
+                             const F77_INT&, F77_REAL*,
+                             const F77_INT&, F77_INT&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (zormlq, ZORMLQ) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const F77_INT&, const F77_INT&, const F77_INT&,
+                             F77_DBLE_CMPLX*, const F77_INT&,
+                             F77_DBLE_CMPLX*, F77_DBLE_CMPLX*,
+                             const F77_INT&, F77_DBLE_CMPLX*,
+                             const F77_INT&, F77_INT&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  // ORMQR
+
+  F77_RET_T
+  F77_FUNC (cormqr, CORMQR) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const F77_INT&, const F77_INT&, const F77_INT&,
+                             F77_CMPLX*, const F77_INT&,
+                             F77_CMPLX*, F77_CMPLX*,
+                             const F77_INT&, F77_CMPLX*,
+                             const F77_INT&, F77_INT&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dormqr, DORMQR) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const F77_INT&, const F77_INT&, const F77_INT&,
+                             F77_DBLE*, const F77_INT&,
+                             F77_DBLE*, F77_DBLE*,
+                             const F77_INT&, F77_DBLE*,
+                             const F77_INT&, F77_INT&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (sormqr, SORMQR) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const F77_INT&, const F77_INT&, const F77_INT&,
+                             F77_REAL*, const F77_INT&,
+                             F77_REAL*, F77_REAL*,
+                             const F77_INT&, F77_REAL*,
+                             const F77_INT&, F77_INT&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (zormqr, ZORMQR) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const F77_INT&, const F77_INT&, const F77_INT&,
+                             F77_DBLE_CMPLX*, const F77_INT&,
+                             F77_DBLE_CMPLX*, F77_DBLE_CMPLX*,
+                             const F77_INT&, F77_DBLE_CMPLX*,
+                             const F77_INT&, F77_INT&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
   // GESDD
 
   F77_RET_T
@@ -351,6 +469,90 @@
                              F77_CHAR_ARG_LEN_DECL
                              F77_CHAR_ARG_LEN_DECL);
 
+  // GEJSV
+
+  F77_RET_T
+  F77_FUNC (cgejsv, CGEJSV) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const F77_INT&, const F77_INT&,
+                             F77_CMPLX*, const F77_INT&, F77_REAL*,
+                             F77_CMPLX*, const F77_INT&,
+                             F77_CMPLX*, const F77_INT&,
+                             F77_CMPLX*, const F77_INT&,
+                             F77_REAL*,  const F77_INT&,
+                             F77_INT *, F77_INT&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dgejsv, DGEJSV) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const F77_INT&, const F77_INT&,
+                             F77_DBLE*, const F77_INT&, F77_DBLE*,
+                             F77_DBLE*, const F77_INT&,
+                             F77_DBLE*, const F77_INT&,
+                             F77_DBLE*, const F77_INT&,
+                             F77_INT *, F77_INT&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (sgejsv, SGEJSV) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const F77_INT&, const F77_INT&,
+                             F77_REAL*, const F77_INT&, F77_REAL*,
+                             F77_REAL*, const F77_INT&,
+                             F77_REAL*, const F77_INT&,
+                             F77_REAL*, const F77_INT&,
+                             F77_INT *, F77_INT&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (zgejsv, ZGEJSV) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const F77_INT&, const F77_INT&,
+                             F77_DBLE_CMPLX*, const F77_INT&, F77_DBLE*,
+                             F77_DBLE_CMPLX*, const F77_INT&,
+                             F77_DBLE_CMPLX*, const F77_INT&,
+                             F77_DBLE_CMPLX*, const F77_INT&,
+                             F77_DBLE*,       const F77_INT&,
+                             F77_INT *, F77_INT&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
   // GEESX
 
   typedef F77_INT (*double_selector) (const F77_DBLE&, const F77_DBLE&);
diff -r f85d491d2322 -r 5438e38e5cb0 liboctave/numeric/svd.cc
--- a/liboctave/numeric/svd.cc	Fri Feb 15 08:29:56 2019 -0800
+++ b/liboctave/numeric/svd.cc	Mon Feb 18 09:08:16 2019 +0800
@@ -39,6 +39,266 @@
 #include "lo-lapack-proto.h"
 #include "svd.h"
 
+template<typename T>
+class
+Get_jsv_lwork
+{
+public:
+  Get_jsv_lwork() = delete;
+
+  // Unfortrainately, dgejsv and sgejsv do not provide estimation of 'lwork'.
+  // Thus we have to estimate it according to corresponding LAPACK document
+  // and related source code (e.g. cgejsv).
+  // In LAPACKE (C interface to LAPACK), the memory handling code in 
+  // LAPACKE_dgejsv() (lapacke_dgejsv.c, last visit 2019-02-17) uses 
+  // the minimum required working space. In contracst, here the optimal 
+  // working space size is computed, at the cost of longer code.
+
+  static F77_INT optimal(char& joba, char& jobu, char& jobv,
+                         F77_INT m, F77_INT n);
+
+private:
+  typedef typename T::element_type P;
+  
+  static F77_INT opt_lwork_geqp3(F77_INT m, F77_INT n,
+                                 P *a, F77_INT lda,
+                                 F77_INT* jpvt, P *tau, P *work,
+                                 F77_INT lwork, F77_INT& info);
+  
+  static F77_INT opt_lwork_geqrf(F77_INT m, F77_INT n,
+                                 P *a, F77_INT lda,
+                                 P *tau, P *work,
+                                 F77_INT lwork, F77_INT& info);
+  
+  static F77_INT opt_lwork_gelqf(F77_INT m, F77_INT n,
+                                 P *a, F77_INT lda,
+                                 P *tau, P *work,
+                                 F77_INT lwork, F77_INT& info);
+  
+  static F77_INT opt_lwork_ormlq(char& side, char& trans,
+                                 F77_INT m, F77_INT n, F77_INT k,
+                                 P *a, F77_INT lda,
+                                 P *tau, P *c, F77_INT ldc,
+                                 P *work, F77_INT lwork, F77_INT& info);
+  
+  static F77_INT opt_lwork_ormqr(char& side, char& trans,
+                                 F77_INT m, F77_INT n, F77_INT k,
+                                 P *a, F77_INT lda,
+                                 P *tau, P *c, F77_INT ldc,
+                                 P *work, F77_INT lwork, F77_INT& info);
+};
+
+#define GEJSV_REAL_QP3_LWORK(f, F)                          \
+    F77_XFCN (f, F, (m, n, a, lda, jpvt,                    \
+                     tau, work, lwork, info))
+
+#define GEJSV_REAL_QR_LWORK(f, F)                           \
+    F77_XFCN (f, F, (m, n, a, lda,                          \
+                     tau, work, lwork, info))
+
+#define GEJSV_REAL_ORM_LWORK(f, F)                          \
+    F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&side,  1),       \
+                     F77_CONST_CHAR_ARG2 (&trans, 1),       \
+                     m, n, k, a, lda, tau,                  \
+                     c, ldc, work, lwork, info              \
+                     F77_CHAR_ARG_LEN (1)                   \
+                     F77_CHAR_ARG_LEN (1)))
+
+// For Matrix
+template<>
+F77_INT
+Get_jsv_lwork<Matrix>::opt_lwork_geqp3(F77_INT m, F77_INT n,
+                                       P *a, F77_INT lda,
+                                       F77_INT* jpvt, P *tau, P *work,
+                                       F77_INT lwork, F77_INT& info)
+{
+  GEJSV_REAL_QP3_LWORK(dgeqp3, DGEQP3);
+  return static_cast<F77_INT>(work[0]);
+}
+
+template<>
+F77_INT
+Get_jsv_lwork<Matrix>::opt_lwork_geqrf(F77_INT m, F77_INT n,
+                                       P *a, F77_INT lda,
+                                       P *tau, P *work,
+                                       F77_INT lwork, F77_INT& info)
+{
+  GEJSV_REAL_QR_LWORK(dgeqrf, DGEQRF);
+  return static_cast<F77_INT>(work[0]);
+}
+
+template<>
+F77_INT
+Get_jsv_lwork<Matrix>::opt_lwork_gelqf(F77_INT m, F77_INT n,
+                                       P *a, F77_INT lda,
+                                       P *tau, P *work,
+                                       F77_INT lwork, F77_INT& info)
+{
+  GEJSV_REAL_QR_LWORK(dgelqf, DGELQF);
+  return static_cast<F77_INT>(work[0]);
+}
+
+template<>
+F77_INT
+Get_jsv_lwork<Matrix>::opt_lwork_ormlq(char& side, char& trans,
+                                       F77_INT m, F77_INT n, F77_INT k,
+                                       P *a, F77_INT lda,
+                                       P *tau, P *c, F77_INT ldc,
+                                       P *work, F77_INT lwork, F77_INT& info)
+{
+  GEJSV_REAL_ORM_LWORK(dormlq, DORMLQ);
+  return static_cast<F77_INT>(work[0]);
+}
+
+template<>
+F77_INT
+Get_jsv_lwork<Matrix>::opt_lwork_ormqr(char& side, char& trans,
+                                       F77_INT m, F77_INT n, F77_INT k,
+                                       P *a, F77_INT lda,
+                                       P *tau, P *c, F77_INT ldc,
+                                       P *work, F77_INT lwork, F77_INT& info)
+{
+  GEJSV_REAL_ORM_LWORK(dormqr, DORMQR);
+  return static_cast<F77_INT>(work[0]);
+}
+
+// For FloatMatrix
+template<>
+F77_INT
+Get_jsv_lwork<FloatMatrix>::opt_lwork_geqp3(F77_INT m, F77_INT n,
+                                            P *a, F77_INT lda,
+                                            F77_INT* jpvt, P *tau, P *work,
+                                            F77_INT lwork, F77_INT& info)
+{
+  GEJSV_REAL_QP3_LWORK(sgeqp3, SGEQP3);
+  return static_cast<F77_INT>(work[0]);
+}
+
+template<>
+F77_INT
+Get_jsv_lwork<FloatMatrix>::opt_lwork_geqrf(F77_INT m, F77_INT n,
+                                            P *a, F77_INT lda,
+                                            P *tau, P *work,
+                                            F77_INT lwork, F77_INT& info)
+{
+  GEJSV_REAL_QR_LWORK(sgeqrf, SGEQRF);
+  return static_cast<F77_INT>(work[0]);
+}
+
+template<>
+F77_INT
+Get_jsv_lwork<FloatMatrix>::opt_lwork_gelqf(F77_INT m, F77_INT n,
+                                            P *a, F77_INT lda,
+                                            P *tau, P *work,
+                                            F77_INT lwork, F77_INT& info)
+{
+  GEJSV_REAL_QR_LWORK(sgelqf, SGELQF);
+  return static_cast<F77_INT>(work[0]);
+}
+
+template<>
+F77_INT
+Get_jsv_lwork<FloatMatrix>::opt_lwork_ormlq(char& side, char& trans,
+                                            F77_INT m, F77_INT n, F77_INT k,
+                                            P *a, F77_INT lda,
+                                            P *tau, P *c, F77_INT ldc,
+                                            P *work, F77_INT lwork, F77_INT& info)
+{
+  GEJSV_REAL_ORM_LWORK(sormlq, SORMLQ);
+  return static_cast<F77_INT>(work[0]);
+}
+
+template<>
+F77_INT
+Get_jsv_lwork<FloatMatrix>::opt_lwork_ormqr(char& side, char& trans,
+                                            F77_INT m, F77_INT n, F77_INT k,
+                                            P *a, F77_INT lda,
+                                            P *tau, P *c, F77_INT ldc,
+                                            P *work, F77_INT lwork, F77_INT& info)
+{
+  GEJSV_REAL_ORM_LWORK(sormqr, SORMQR);
+  return static_cast<F77_INT>(work[0]);
+}
+
+#undef GEJSV_REAL_QP3_LWORK
+#undef GEJSV_REAL_QR_LWORK
+#undef GEJSV_REAL_ORM_LWORK
+
+template<typename T>
+F77_INT
+Get_jsv_lwork<T>::optimal(char& joba, char& jobu, char& jobv,
+                          F77_INT m, F77_INT n)
+{
+  using std::max;
+  
+  F77_INT lwork = -1;             // work space size query
+  std::vector<P> work(2);
+  
+  // variables that mimic running environment of gejsv
+  F77_INT lda  = max (m, 1);
+  F77_INT ierr = 0;
+  char side  = 'L';
+  char trans = 'N';
+  std::vector<P> mat_a(1);
+  P *a = mat_a.data ();           // no real data is required.
+  std::vector<F77_INT> vec_jpvt = {0};
+  P *tau = work.data ();
+  P *u   = work.data ();
+  P *v   = work.data ();
+
+  bool need_lsvec = jobu == 'U' || jobu == 'F';
+  bool need_rsvec = jobv == 'V' || jobv == 'J';
+
+  F77_INT lwork_pocon = 3 * n;  // for [s,d]pocon
+  F77_INT lwork_geqp3 = opt_lwork_geqp3(m, n, a, lda, vec_jpvt.data (),
+                                        tau, work.data (), -1, ierr);
+  F77_INT lwork_geqrf = opt_lwork_geqrf(m, n, a, lda,
+                                        tau, work.data (), -1, ierr);
+  
+  if (! (need_lsvec || need_rsvec) )
+    {
+      // only SIGMA is needed
+      if (! (joba == 'E' || joba == 'G') )
+        lwork = max<int>({2*m + n, n + lwork_geqp3, n + lwork_geqrf, 7});
+      else
+        lwork = max<int>({2*m + n, n + lwork_geqp3, n + lwork_geqrf,
+                          n + n*n + lwork_pocon, 7});
+    }
+  else if (need_rsvec && ! need_lsvec)
+    {
+      // SIGMA and the right singular vectors are needed
+      F77_INT lwork_gelqf = opt_lwork_gelqf(n, n, a, lda, tau, work.data (), -1, ierr);
+      trans = 'T';
+      F77_INT lwork_ormlq = opt_lwork_ormlq(side, trans, n, n, n, a, lda,
+                                            tau, v, n, work.data (), -1, ierr);
+      lwork = max<int>({2*m + n, n + lwork_geqp3, n + lwork_pocon,
+                        n + lwork_gelqf, 2*n + lwork_geqrf, n + lwork_ormlq});
+    }
+  else if (need_lsvec && ! need_rsvec)
+    {
+      // SIGMA and the left singular vectors are needed
+      F77_INT n1 = (jobu == 'U') ? n : m;  // size of U is m x n1
+      F77_INT lwork_ormqr = opt_lwork_ormqr(side, trans, m, n1, n, a, lda,
+                                            tau, u, m, work.data (), -1, ierr);
+      lwork = max<int>({2*m + n, n + lwork_geqp3, n + lwork_pocon,
+                        2*n + lwork_geqrf, n + lwork_ormqr});
+    }
+  else  // full SVD is needed
+    {
+      if (jobv == 'V')
+        lwork = max<int>({2*m + n, 6*n + 2*n*n});
+      else if (jobv == 'J')
+        lwork = max<int>({2*m + n, 4*n + n*n, 2*n + n*n + 6});
+
+      F77_INT n1 = (jobu == 'U') ? n : m;  // size of U is m x n1
+      F77_INT lwork_ormqr = opt_lwork_ormqr(side, trans, m, n1, n, a, lda,
+                                            tau, u, m, work.data (), -1, ierr);
+      lwork = max<int>({lwork, n + lwork_ormqr});
+    }
+
+  return lwork;
+}
+
 namespace octave
 {
   namespace math
@@ -272,6 +532,139 @@
 #undef GESDD_REAL_STEP
 #undef GESDD_COMPLEX_STEP
 
+
+    // GEJSV specializations
+
+#define GEJSV_REAL_STEP(f, F)                                           \
+    F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&joba, 1),                    \
+                     F77_CONST_CHAR_ARG2 (&jobu, 1),                    \
+                     F77_CONST_CHAR_ARG2 (&jobv, 1),                    \
+                     F77_CONST_CHAR_ARG2 (&jobr, 1),                    \
+                     F77_CONST_CHAR_ARG2 (&jobt, 1),                    \
+                     F77_CONST_CHAR_ARG2 (&jobp, 1),                    \
+                     m, n, tmp_data, m1, s_vec, u, m1, v, nrow_v1,      \
+                     work.data (), lwork, iwork.data (), info           \
+                     F77_CHAR_ARG_LEN (1)                               \
+                     F77_CHAR_ARG_LEN (1)                               \
+                     F77_CHAR_ARG_LEN (1)                               \
+                     F77_CHAR_ARG_LEN (1)                               \
+                     F77_CHAR_ARG_LEN (1)                               \
+                     F77_CHAR_ARG_LEN (1)))
+
+#define GEJSV_COMPLEX_STEP(f, F, CMPLX_ARG)                       \
+    F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&joba, 1),              \
+                     F77_CONST_CHAR_ARG2 (&jobu, 1),              \
+                     F77_CONST_CHAR_ARG2 (&jobv, 1),              \
+                     F77_CONST_CHAR_ARG2 (&jobr, 1),              \
+                     F77_CONST_CHAR_ARG2 (&jobt, 1),              \
+                     F77_CONST_CHAR_ARG2 (&jobp, 1),              \
+                     m, n, CMPLX_ARG (tmp_data), m1,              \
+                     s_vec, CMPLX_ARG (u), m1,                    \
+                     CMPLX_ARG (v), nrow_v1,                      \
+                     CMPLX_ARG (work.data ()), lwork,             \
+                     rwork.data (), lrwork, iwork.data (), info   \
+                     F77_CHAR_ARG_LEN (1)                         \
+                     F77_CHAR_ARG_LEN (1)                         \
+                     F77_CHAR_ARG_LEN (1)                         \
+                     F77_CHAR_ARG_LEN (1)                         \
+                     F77_CHAR_ARG_LEN (1)                         \
+                     F77_CHAR_ARG_LEN (1)))
+
+    // DGEJSV
+    template<>
+    void
+    svd<Matrix>::gejsv (
+      char& joba, char& jobu, char& jobv, char& jobr, char& jobt,
+      char& jobp, F77_INT m, F77_INT n,
+      P *tmp_data, F77_INT m1, DM_P *s_vec, P *u,
+      P *v, F77_INT nrow_v1, std::vector<P>& work,
+      F77_INT& lwork, std::vector<F77_INT>& iwork,
+      F77_INT& info)
+    {
+      lwork = Get_jsv_lwork<Matrix>::optimal(joba, jobu, jobv, m, n);
+      work.reserve (lwork);
+      
+      GEJSV_REAL_STEP (dgejsv, DGEJSV);
+    }
+
+    // SGEJSV
+    template<>
+    void
+    svd<FloatMatrix>::gejsv (
+      char& joba, char& jobu, char& jobv, char& jobr, char& jobt,
+      char& jobp, F77_INT m, F77_INT n,
+      P *tmp_data, F77_INT m1, DM_P *s_vec, P *u,
+      P *v, F77_INT nrow_v1, std::vector<P>& work,
+      F77_INT& lwork, std::vector<F77_INT>& iwork,
+      F77_INT& info)
+    {
+      lwork = Get_jsv_lwork<FloatMatrix>::optimal(joba, jobu, jobv, m, n);
+      work.reserve (lwork);
+
+      GEJSV_REAL_STEP (sgejsv, SGEJSV);
+    }
+
+    // ZGEJSV
+    template<>
+    void
+    svd<ComplexMatrix>::gejsv (
+      char& joba, char& jobu, char& jobv, char& jobr, char& jobt,
+      char& jobp, F77_INT m, F77_INT n,
+      P *tmp_data, F77_INT m1, DM_P *s_vec, P *u,
+      P *v, F77_INT nrow_v1, std::vector<P>& work,
+      F77_INT& lwork, std::vector<F77_INT>& iwork,
+      F77_INT& info)
+    {
+      F77_INT lrwork = -1;          // work space size query
+      std::vector<double> rwork (1);
+      work.reserve (2);
+
+      GEJSV_COMPLEX_STEP (zgejsv, ZGEJSV, F77_DBLE_CMPLX_ARG);
+      
+      lwork = static_cast<F77_INT> (work[0].real ());
+      work.reserve (lwork);
+      
+      lrwork = static_cast<F77_INT> (rwork[0]);
+      rwork.reserve (lrwork);
+      
+      F77_INT liwork = static_cast<F77_INT> (iwork[0]);
+      iwork.reserve (liwork);
+
+      GEJSV_COMPLEX_STEP (zgejsv, ZGEJSV, F77_DBLE_CMPLX_ARG);
+    }
+
+    // CGEJSV
+    template<>
+    void
+    svd<FloatComplexMatrix>::gejsv (
+      char& joba, char& jobu, char& jobv, char& jobr, char& jobt,
+      char& jobp, F77_INT m, F77_INT n,
+      P *tmp_data, F77_INT m1, DM_P *s_vec, P *u,
+      P *v, F77_INT nrow_v1, std::vector<P>& work,
+      F77_INT& lwork, std::vector<F77_INT>& iwork,
+      F77_INT& info)
+    {
+      F77_INT lrwork = -1;          // work space size query
+      std::vector<float> rwork (1);
+      work.reserve (2);
+
+      GEJSV_COMPLEX_STEP (cgejsv, CGEJSV, F77_CMPLX_ARG);
+      
+      lwork = static_cast<F77_INT> (work[0].real ());
+      work.reserve (lwork);
+      
+      lrwork = static_cast<F77_INT> (rwork[0]);
+      rwork.reserve (lrwork);
+      
+      F77_INT liwork = static_cast<F77_INT> (iwork[0]);
+      iwork.reserve (liwork);
+
+      GEJSV_COMPLEX_STEP (cgejsv, CGEJSV, F77_CMPLX_ARG);
+    }
+
+#undef GEJSV_REAL_STEP
+#undef GEJSV_COMPLEX_STEP
+
     template<typename T>
     svd<T>::svd (const T& a, svd::Type type,
                  svd::Driver driver)
@@ -366,8 +759,9 @@
 
       std::vector<P> work (1);
 
-      F77_INT m1 = std::max (m, static_cast<F77_INT> (1));
-      F77_INT nrow_vt1 = std::max (nrow_vt, static_cast<F77_INT> (1));
+      const F77_INT F77_INT_ONE = static_cast<F77_INT> (1);
+      F77_INT m1 = std::max (m, F77_INT_ONE);
+      F77_INT nrow_vt1 = std::max (nrow_vt, F77_INT_ONE);
 
       if (m_driver == svd::Driver::GESVD)
         gesvd (jobu, jobv, m, n, tmp_data, m1, s_vec, u, vt, nrow_vt1,
@@ -382,6 +776,73 @@
           gesdd (jobz, m, n, tmp_data, m1, s_vec, u, vt, nrow_vt1,
                  work, lwork, iwork.data (), info);
         }
+      else if (m_driver == svd::Driver::GEJSV)
+        {
+          assert (jobu == jobv);  //TODO: allow request of U or V only.
+          char jobz = jobu;
+          
+          bool transposed = false;
+          if (n > m)
+            {
+              // GEJSV only accept m >= n, thus we need transpose here
+              transposed = true;
+
+              std::swap(m, n);
+              m1 = std::max (m, F77_INT_ONE);
+              nrow_vt1 = std::max (n, F77_INT_ONE); // we have m>n
+              if (m_type == svd::Type::sigma_only)
+                nrow_vt1 = 1;
+              std::swap(jobu, jobv);
+
+              atmp = atmp.hermitian ();
+              tmp_data = atmp.fortran_vec ();
+              // Swap pointers of U and V.
+              u  = right_sm.fortran_vec ();
+              vt = left_sm.fortran_vec ();
+            }
+          
+          char gsvd2jsv[256][2] = {0};  //FIXME: use map, if-else, or switch?
+          gsvd2jsv[int('A')][0] = 'F';
+          gsvd2jsv[int('A')][1] = 'J';
+          gsvd2jsv[int('S')][0] = 'U';
+          gsvd2jsv[int('S')][1] = 'V';
+          gsvd2jsv[int('O')][0] = 'W';
+          gsvd2jsv[int('O')][1] = 'W';
+          gsvd2jsv[int('N')][0] = 'N';
+          gsvd2jsv[int('N')][1] = 'N';
+          jobu = gsvd2jsv[int(jobz)][0];
+          jobv = gsvd2jsv[int(jobz)][1];
+          
+          char joba = 'F';  // 'F': most conservative
+          char jobr = 'R';  // 'R' is recommended.
+          char jobt = 'N';  //or 'T', but that require U and V appear together
+          char jobp = 'N';  // 'P' if denormal is poorly implemented.
+          
+          std::vector<F77_INT> iwork (std::max(m + 3*n, 1));
+
+          gejsv (joba, jobu, jobv, jobr, jobt, jobp, m, n, tmp_data, m1,
+                 s_vec, u, vt, nrow_vt1, work, lwork, iwork, info);
+          
+          if (iwork[2] == 1)
+            {
+              printf("GEJSV: denormal occured, possible loss of accuracy\n");
+            }
+          
+          if (info < 0)
+            {
+              printf("GEJSV: illegal argument at #%d\n", -info);
+            }
+          else if (info > 0)
+            {
+              printf("GEJSV: inaccurate result - fail to converge within max sweeps.\n");
+            }
+
+          if (transposed)
+            {
+              std::swap(m, n);
+              // other necessary reverse swap?
+            }
+        }
       else
         (*current_liboctave_error_handler) ("svd: unknown driver");
 
@@ -392,7 +853,8 @@
             sigma.dgxelem (i) = DM_P (0);
         }
 
-      if (! (jobv == 'N' || jobv == 'O'))
+      // GESVD and GESDD return VT instead of V, GEJSV return V.
+      if (! (jobv == 'N' || jobv == 'O') && (m_driver != svd::Driver::GEJSV))
         right_sm = right_sm.hermitian ();
     }
 
diff -r f85d491d2322 -r 5438e38e5cb0 liboctave/numeric/svd.h
--- a/liboctave/numeric/svd.h	Fri Feb 15 08:29:56 2019 -0800
+++ b/liboctave/numeric/svd.h	Mon Feb 18 09:08:16 2019 +0800
@@ -50,7 +50,8 @@
       enum class Driver
       {
         GESVD,
-        GESDD
+        GESDD,
+        GEJSV
       };
 
       svd (void)
@@ -110,6 +111,14 @@
                   P *vt, octave_f77_int_type nrow_vt1, std::vector<P>& work,
                   octave_f77_int_type& lwork, octave_f77_int_type *iwork,
                   octave_f77_int_type& info);
+
+      void gejsv (char& joba, char& jobu, char& jobv, char& jobr, char& jobt,
+                  char& jobp, octave_f77_int_type m, octave_f77_int_type n,
+                  P *tmp_data, octave_f77_int_type m1, DM_P *s_vec, P *u,
+                  P *v, octave_f77_int_type nrow_v1, std::vector<P>& work,
+                  octave_f77_int_type& lwork,
+                  std::vector<octave_f77_int_type>& iwork,
+                  octave_f77_int_type& info);
     };
   }
 }
