xGELSD changes July 2024

[frjohnst]

This PR for xGELSD has two purposes:

1. fix work space checks in "Path 2a" in CGELSD and ZGELSD

2. add quick return for NRHS=0 to all four xGELSD, since lower
   level xLALSD does not accept NRHS=0.

These are decsribed below in more detail:

---

1. fix work space checks in "Path 2a" in CGELSD and ZGELSD

We believe the workspace check for "Path 2a" in
cgelsd.f and zgelsd.f is incorrect. Here is the
check currently in zgelsd.f (around line 531):

  +531        ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
  +532       $         MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN
  +533  *
  +534  *        Path 2a - underdetermined, with many more columns than rows
  +535  *        and sufficient workspace for an efficient algorithm.
  +536  *
"lapack/SRC/zgelsd.f"

We believe
$ MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN
should be changed to
$ MAX( M, 2*M-4, NRHS, N-3*M, M*NRHS ) ) THEN

This part of the workspace will be passed to ZLALSD
which requires workspace with M*NRHS elements:

*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX*16 array, dimension (N * NRHS)
*> \endverbatim
*>
"lapack/SRC/zlalsd.f"

Note that within ZLALSD the name "N" is used for the third argument:

      SUBROUTINE ZLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
     $                   RANK, WORK, RWORK, IWORK, INFO )
*
"./lapack/SRC/zlalsd.f"

while ZGELSD passes "M" to ZLALSD in "Path 2a":

  +578           CALL ZLALSD( 'U', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB,
  +579       $                RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
  +580       $                IWORK, INFO )
"lapack/SRC/zgelsd.f"

To be consistent with the preceding, the following test in Path 2a

  +534  *        Path 2a - underdetermined, with many more columns than rows
  +535  *        and sufficient workspace for an efficient algorithm.
  +536  *
  +537           LDWORK = M
  +538           IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
  +539       $       M*LDA+M+M*NRHS ) )LDWORK = LDA
"lapack/SRC/zgelsd.f"

should be changed to (add M*NRHS to inner MAX list)

*        Path 2a - underdetermined, with many more columns than rows
*        and sufficient workspace for an efficient algorithm.
*
         LDWORK = M
         IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M,
     $       M*NRHS ), M*LDA+M+M*NRHS ) )LDWORK = LDA

Both cgelsd.f and zgelsd.f have these issues.

---

2. add quick return for NRHS=0 to all four xGELSD

As currently coded, CGELSD, DGELSD, SGELSD, ZGELSD
all accept NRHS=0, but lower level routines return
non-zero INFO when passed NRHS=0.

For example, DGELSD accepts NRHS=0

  +256        ELSE IF( NRHS.LT.0 ) THEN
  +257           INFO = -3
"lapack/SRC/dgelsd.f"

but NRHS=0 will be passed to DLALSD, for example,

+490 CALL DLALSD( 'U', SMLSIZ, N, NRHS, S, WORK( IE ), B, LDB,
+491 $ RCOND, RANK, WORK( NWORK ), IWORK, INFO )
"lapack/SRC/dgelsd.f"

and DLALSD does not accept NRHS=0

  +222        ELSE IF( NRHS.LT.1 ) THEN
  +223           INFO = -4
"lapack/SRC/dlalsd.f"

We believe this should be fixed by adding a quick return for the NRHS=0
case to all four xGELSD routines.

For example, change (around line 370 in DGELSD)

*
*     Quick return if possible.
*
      IF( M.EQ.0 .OR. N.EQ.0 ) THEN
         RANK = 0
         RETURN
      END IF

to (add NRHS.EQ.0 to the IF list)

*
*     Quick return if possible.
*
      IF( M.EQ.0 .OR. N.EQ.0 .OR. NRHS.EQ.0) THEN
         RANK = 0
         RETURN
      END IF

To improve the end-user experience, this should be fixed
in CGELSD, DGELSD, SGELSD, ZGELSD.