Re: [PATCH] Add workaround for broken xGELSD workspace queries.

[John W. Eaton]

On 13-Feb-2008, David Bateman wrote:

| John W. Eaton wrote:
| > On 12-Feb-2008, Jason Riedy wrote:
| >
| > | The xGELSD workspace queries are inconsistant with the logic inside
| > | those routines in the short and fat case (m > n).  For that case,
| > | calculate a workspace size large enough to enable the efficient
| > | algorithm.  This may allocate more memory than necessary for small
| > | cases.  The calculated alternate size is at least the documented
| > | minimum workspace requirement, so this amount always will work even
| > | if the routines change.
| > | ---
| > |    Sorry it's been so long.  I don't want to describe the LAPACK
| > |    developer dynamic, but it's severely busted.  ergh.  I think
| > |    this is a minimal change, and the ILAENV declaration should
| > |    work on every Fortran compiler I know. -- Jason
| >
| > Just to avoid possible incompatible C<->Fortran calling conventions
| > for Fortran funtions, we have always wrapped them in Fortran
| > subroutines (see the files in libcruft/lapack-xtra).  I did that for
| > ilaenv and checked in your change.
| >
| > Thanks,
| >
| > jwe
| >
| >   
| This supercedes the xGELSD change I sent. However I have one
| disagreement with this patch. The xGELSD.f functions use the condition
| 
|                IF( N.GE.MNTHR ) THEN
| 
| rather that (n > m) as Jason used in his patch to check for the
| workaround. Where MNTHR is given by
| 
|             MNTHR = ILAENV( 6, 'ZGELSD', ' ', M, N, NRHS, -1 )
| 
| or 1.6 *m. Therefore I believe that the workaround should only apply
| with the condition "if (1.6 * n > m)" instead of teh "if (n > m)" that
| is currently used.

OK, does the following look OK?  It's in my public archive now, along
with another slight modification so we check n >= mnthr instead of
just m > mnthr.

Thanks,

jwe


# HG changeset patch
# User John W. Eaton <address@hidden>
# Date 1203119745 18000
# Branch release-3-0-x
# Node ID a39ca14faf2d97d420d5177605a66dca7c8defff
# Parent  096c3353e6929779d7bc966cac6d27eb57f344cf
more xGELSD workspace fixes

diff -r 096c3353e692 -r a39ca14faf2d liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc      Fri Feb 15 16:50:16 2008 -0500
+++ b/liboctave/CMatrix.cc      Fri Feb 15 18:55:45 2008 -0500
@@ -2507,6 +2507,13 @@ ComplexMatrix::lssolve (const ComplexMat
                                   F77_CHAR_ARG_LEN (6)
                                   F77_CHAR_ARG_LEN (1));
 
+      octave_idx_type mnthr;
+      F77_FUNC (xilaenv, XILAENV) (6, F77_CONST_CHAR_ARG2 ("ZGELSD", 6),
+                                  F77_CONST_CHAR_ARG2 (" ", 1),
+                                  m, n, nrhs, -1, mnthr
+                                  F77_CHAR_ARG_LEN (6)
+                                  F77_CHAR_ARG_LEN (1));
+
       // We compute the size of rwork and iwork because ZGELSD in
       // older versions of LAPACK does not return them on a query
       // call.
@@ -2539,10 +2546,10 @@ ComplexMatrix::lssolve (const ComplexMat
                                 lwork, prwork, 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
+      // through 3.1.1 when n > mnthr.  The obtuse formula below
+      // should provide sufficient workspace for ZGELSD to operate
       // efficiently.
-      if (n > m)
+      if (n > mnthr)
        {
          octave_idx_type addend = m;
 
diff -r 096c3353e692 -r a39ca14faf2d liboctave/ChangeLog
--- a/liboctave/ChangeLog       Fri Feb 15 16:50:16 2008 -0500
+++ b/liboctave/ChangeLog       Fri Feb 15 18:55:45 2008 -0500
@@ -1,3 +1,9 @@ 2008-02-12  John W. Eaton  <address@hidden
+2008-02-15  John W. Eaton  <address@hidden>
+
+       * dMatrix.cc (Matrix::lssolve): Check n > mnthr, not n > m when
+       deciding whether to calculate workspace size, with mnthr from ILAENV.
+       * CMatrix.cc (ComplexMatrix::lssolve): Likewise.
+
 2008-02-12  John W. Eaton  <address@hidden>
 
        * CMatrix.cc: Declare xilaenv instead of ilaenv.
diff -r 096c3353e692 -r a39ca14faf2d liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc      Fri Feb 15 16:50:16 2008 -0500
+++ b/liboctave/dMatrix.cc      Fri Feb 15 18:55:45 2008 -0500
@@ -2118,6 +2118,13 @@ Matrix::lssolve (const Matrix& b, octave
                                   F77_CHAR_ARG_LEN (6)
                                   F77_CHAR_ARG_LEN (1));
 
+      octave_idx_type mnthr;
+      F77_FUNC (xilaenv, XILAENV) (6, F77_CONST_CHAR_ARG2 ("DGELSD", 6),
+                                  F77_CONST_CHAR_ARG2 (" ", 1),
+                                  m, n, nrhs, -1, mnthr
+                                  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.
       double dminmn = static_cast<double> (minmn);
@@ -2142,10 +2149,10 @@ Matrix::lssolve (const Matrix& b, octave
                                 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
+      // through 3.1.1 when n > mnthr.  The obtuse formula below
       // should provide sufficient workspace for DGELSD to operate
       // efficiently.
-      if (n > m)
+      if (n > mnthr)
        {
          const octave_idx_type wlalsd
            = 9*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs + (smlsiz+1)*(smlsiz+1);