--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -59,6 +59,13 @@
 
 extern "C"
 {
+  octave_idx_type
+  F77_FUNC (ilaenv, ILAENV) (const octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
+			     F77_CONST_CHAR_ARG_DECL,
+			     const octave_idx_type&, const octave_idx_type&,
+			     const octave_idx_type&, const octave_idx_type&
+			     F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL);
+
   F77_RET_T
   F77_FUNC (dgebal, DGEBAL) (F77_CONST_CHAR_ARG_DECL,
 			     const octave_idx_type&, double*, const octave_idx_type&, octave_idx_type&,
@@ -2103,8 +2110,12 @@
 
       Array<double> work (1);
 
-      // FIXME: Can SMLSIZ be other than 25?
-      octave_idx_type smlsiz = 25;
+      const octave_idx_type smlsiz
+	= F77_FUNC (ilaenv, ILAENV) (9, F77_CONST_CHAR_ARG2 ("DGELSD", 6),
+				     F77_CONST_CHAR_ARG2 (" ", 1),
+				     0, 0, 0, 0
+				     F77_CHAR_ARG_LEN (6)
+				     F77_CHAR_ARG_LEN (1));
 
       // We compute the size of iwork because DGELSD in older versions
       // of LAPACK does not return it on a query call.
@@ -2129,6 +2140,35 @@
 				 ps, rcond, rank, work.fortran_vec (),
 				 lwork, piwork, info));
 
+      // The workspace query is broken in at least LAPACK 3.0.0
+      // through 3.1.1 when n > m.  The obtuse formula below
+      // should provide sufficient workspace for DGELSD to operate
+      // efficiently.
+      if (n > m)
+	{
+	  const octave_idx_type wlalsd
+	    = 9*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs + (smlsiz+1)*(smlsiz+1);
+
+	  octave_idx_type addend = m;
+
+	  if (2*m-4 > addend)
+	    addend = 2*m-4;
+
+	  if (nrhs > addend)
+	    addend = nrhs;
+
+	  if (n-3*m > addend)
+	    addend = n-3*m;
+
+	  if (wlalsd > addend)
+	    addend = wlalsd;
+
+	  const octave_idx_type lworkaround = 4*m + m*m + addend;
+
+	  if (work(0) < lworkaround)
+	    work(0) = lworkaround;
+	}
+
       if (f77_exception_encountered)
 	(*current_liboctave_error_handler) 
 	  ("unrecoverable error in dgelsd");