SRC/cgelss.f

   1 *> \brief <b> CGELSS solves overdetermined or underdetermined systems for GE matrices</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CGELSS + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgelss.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgelss.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgelss.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
  22 *                          WORK, LWORK, RWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
  26 *       REAL               RCOND
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       REAL               RWORK( * ), S( * )
  30 *       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> CGELSS computes the minimum norm solution to a complex linear
  40 *> least squares problem:
  41 *>
  42 *> Minimize 2-norm(| b - A*x |).
  43 *>
  44 *> using the singular value decomposition (SVD) of A. A is an M-by-N
  45 *> matrix which may be rank-deficient.
  46 *>
  47 *> Several right hand side vectors b and solution vectors x can be
  48 *> handled in a single call; they are stored as the columns of the
  49 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix
  50 *> X.
  51 *>
  52 *> The effective rank of A is determined by treating as zero those
  53 *> singular values which are less than RCOND times the largest singular
  54 *> value.
  55 *> \endverbatim
  56 *
  57 *  Arguments:
  58 *  ==========
  59 *
  60 *> \param[in] M
  61 *> \verbatim
  62 *>          M is INTEGER
  63 *>          The number of rows of the matrix A. M >= 0.
  64 *> \endverbatim
  65 *>
  66 *> \param[in] N
  67 *> \verbatim
  68 *>          N is INTEGER
  69 *>          The number of columns of the matrix A. N >= 0.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] NRHS
  73 *> \verbatim
  74 *>          NRHS is INTEGER
  75 *>          The number of right hand sides, i.e., the number of columns
  76 *>          of the matrices B and X. NRHS >= 0.
  77 *> \endverbatim
  78 *>
  79 *> \param[in,out] A
  80 *> \verbatim
  81 *>          A is COMPLEX array, dimension (LDA,N)
  82 *>          On entry, the M-by-N matrix A.
  83 *>          On exit, the first min(m,n) rows of A are overwritten with
  84 *>          its right singular vectors, stored rowwise.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] LDA
  88 *> \verbatim
  89 *>          LDA is INTEGER
  90 *>          The leading dimension of the array A. LDA >= max(1,M).
  91 *> \endverbatim
  92 *>
  93 *> \param[in,out] B
  94 *> \verbatim
  95 *>          B is COMPLEX array, dimension (LDB,NRHS)
  96 *>          On entry, the M-by-NRHS right hand side matrix B.
  97 *>          On exit, B is overwritten by the N-by-NRHS solution matrix X.
  98 *>          If m >= n and RANK = n, the residual sum-of-squares for
  99 *>          the solution in the i-th column is given by the sum of
 100 *>          squares of the modulus of elements n+1:m in that column.
 101 *> \endverbatim
 102 *>
 103 *> \param[in] LDB
 104 *> \verbatim
 105 *>          LDB is INTEGER
 106 *>          The leading dimension of the array B.  LDB >= max(1,M,N).
 107 *> \endverbatim
 108 *>
 109 *> \param[out] S
 110 *> \verbatim
 111 *>          S is REAL array, dimension (min(M,N))
 112 *>          The singular values of A in decreasing order.
 113 *>          The condition number of A in the 2-norm = S(1)/S(min(m,n)).
 114 *> \endverbatim
 115 *>
 116 *> \param[in] RCOND
 117 *> \verbatim
 118 *>          RCOND is REAL
 119 *>          RCOND is used to determine the effective rank of A.
 120 *>          Singular values S(i) <= RCOND*S(1) are treated as zero.
 121 *>          If RCOND < 0, machine precision is used instead.
 122 *> \endverbatim
 123 *>
 124 *> \param[out] RANK
 125 *> \verbatim
 126 *>          RANK is INTEGER
 127 *>          The effective rank of A, i.e., the number of singular values
 128 *>          which are greater than RCOND*S(1).
 129 *> \endverbatim
 130 *>
 131 *> \param[out] WORK
 132 *> \verbatim
 133 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
 134 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 135 *> \endverbatim
 136 *>
 137 *> \param[in] LWORK
 138 *> \verbatim
 139 *>          LWORK is INTEGER
 140 *>          The dimension of the array WORK. LWORK >= 1, and also:
 141 *>          LWORK >=  2*min(M,N) + max(M,N,NRHS)
 142 *>          For good performance, LWORK should generally be larger.
 143 *>
 144 *>          If LWORK = -1, then a workspace query is assumed; the routine
 145 *>          only calculates the optimal size of the WORK array, returns
 146 *>          this value as the first entry of the WORK array, and no error
 147 *>          message related to LWORK is issued by XERBLA.
 148 *> \endverbatim
 149 *>
 150 *> \param[out] RWORK
 151 *> \verbatim
 152 *>          RWORK is REAL array, dimension (5*min(M,N))
 153 *> \endverbatim
 154 *>
 155 *> \param[out] INFO
 156 *> \verbatim
 157 *>          INFO is INTEGER
 158 *>          = 0:  successful exit
 159 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 160 *>          > 0:  the algorithm for computing the SVD failed to converge;
 161 *>                if INFO = i, i off-diagonal elements of an intermediate
 162 *>                bidiagonal form did not converge to zero.
 163 *> \endverbatim
 164 *
 165 *  Authors:
 166 *  ========
 167 *
 168 *> \author Univ. of Tennessee
 169 *> \author Univ. of California Berkeley
 170 *> \author Univ. of Colorado Denver
 171 *> \author NAG Ltd.
 172 *
 173 *> \date June 2016
 174 *
 175 *> \ingroup complexGEsolve
 176 *
 177 *  =====================================================================
 178       SUBROUTINE CGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
 179      $                   WORK, LWORK, RWORK, INFO )
 180 *
 181 *  -- LAPACK driver routine (version 3.6.1) --
 182 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 183 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 184 *     June 2016
 185 *
 186 *     .. Scalar Arguments ..
 187       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
 188       REAL               RCOND
 189 *     ..
 190 *     .. Array Arguments ..
 191       REAL               RWORK( * ), S( * )
 192       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
 193 *     ..
 194 *
 195 *  =====================================================================
 196 *
 197 *     .. Parameters ..
 198       REAL               ZERO, ONE
 199       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 200       COMPLEX            CZERO, CONE
 201       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
 202      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
 203 *     ..
 204 *     .. Local Scalars ..
 205       LOGICAL            LQUERY
 206       INTEGER            BL, CHUNK, I, IASCL, IBSCL, IE, IL, IRWORK,
 207      $                   ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN,
 208      $                   MAXWRK, MINMN, MINWRK, MM, MNTHR
 209       INTEGER            LWORK_CGEQRF, LWORK_CUNMQR, LWORK_CGEBRD,
 210      $                   LWORK_CUNMBR, LWORK_CUNGBR, LWORK_CUNMLQ,
 211      $                   LWORK_CGELQF
 212       REAL               ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
 213 *     ..
 214 *     .. Local Arrays ..
 215       COMPLEX            DUM( 1 )
 216 *     ..
 217 *     .. External Subroutines ..
 218       EXTERNAL           CBDSQR, CCOPY, CGEBRD, CGELQF, CGEMM, CGEMV,
 219      $                   CGEQRF, CLACPY, CLASCL, CLASET, CSRSCL, CUNGBR,
 220      $                   CUNMBR, CUNMLQ, CUNMQR, SLABAD, SLASCL, SLASET,
 221      $                   XERBLA
 222 *     ..
 223 *     .. External Functions ..
 224       INTEGER            ILAENV
 225       REAL               CLANGE, SLAMCH
 226       EXTERNAL           ILAENV, CLANGE, SLAMCH
 227 *     ..
 228 *     .. Intrinsic Functions ..
 229       INTRINSIC          MAX, MIN
 230 *     ..
 231 *     .. Executable Statements ..
 232 *
 233 *     Test the input arguments
 234 *
 235       INFO = 0
 236       MINMN = MIN( M, N )
 237       MAXMN = MAX( M, N )
 238       LQUERY = ( LWORK.EQ.-1 )
 239       IF( M.LT.0 ) THEN
 240          INFO = -1
 241       ELSE IF( N.LT.0 ) THEN
 242          INFO = -2
 243       ELSE IF( NRHS.LT.0 ) THEN
 244          INFO = -3
 245       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 246          INFO = -5
 247       ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
 248          INFO = -7
 249       END IF
 250 *
 251 *     Compute workspace
 252 *      (Note: Comments in the code beginning "Workspace:" describe the
 253 *       minimal amount of workspace needed at that point in the code,
 254 *       as well as the preferred amount for good performance.
 255 *       CWorkspace refers to complex workspace, and RWorkspace refers
 256 *       to real workspace. NB refers to the optimal block size for the
 257 *       immediately following subroutine, as returned by ILAENV.)
 258 *
 259       IF( INFO.EQ.0 ) THEN
 260          MINWRK = 1
 261          MAXWRK = 1
 262          IF( MINMN.GT.0 ) THEN
 263             MM = M
 264             MNTHR = ILAENV( 6, 'CGELSS', ' ', M, N, NRHS, -1 )
 265             IF( M.GE.N .AND. M.GE.MNTHR ) THEN
 266 *
 267 *              Path 1a - overdetermined, with many more rows than
 268 *                        columns
 269 *
 270 *              Compute space needed for CGEQRF
 271                CALL CGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO )
 272                LWORK_CGEQRF=DUM(1)
 273 *              Compute space needed for CUNMQR
 274                CALL CUNMQR( 'L', 'C', M, NRHS, N, A, LDA, DUM(1), B,
 275      $                   LDB, DUM(1), -1, INFO )
 276                LWORK_CUNMQR=DUM(1)
 277                MM = N
 278                MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'CGEQRF', ' ', M,
 279      $                       N, -1, -1 ) )
 280                MAXWRK = MAX( MAXWRK, N + NRHS*ILAENV( 1, 'CUNMQR', 'LC',
 281      $                       M, NRHS, N, -1 ) )
 282             END IF
 283             IF( M.GE.N ) THEN
 284 *
 285 *              Path 1 - overdetermined or exactly determined
 286 *
 287 *              Compute space needed for CGEBRD
 288                CALL CGEBRD( MM, N, A, LDA, S, S, DUM(1), DUM(1), DUM(1),
 289      $                      -1, INFO )
 290                LWORK_CGEBRD=DUM(1)
 291 *              Compute space needed for CUNMBR
 292                CALL CUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, DUM(1),
 293      $                B, LDB, DUM(1), -1, INFO )
 294                LWORK_CUNMBR=DUM(1)
 295 *              Compute space needed for CUNGBR
 296                CALL CUNGBR( 'P', N, N, N, A, LDA, DUM(1),
 297      $                   DUM(1), -1, INFO )
 298                LWORK_CUNGBR=DUM(1)
 299 *              Compute total workspace needed
 300                MAXWRK = MAX( MAXWRK, 2*N + LWORK_CGEBRD )
 301                MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNMBR )
 302                MAXWRK = MAX( MAXWRK, 2*N + LWORK_CUNGBR )
 303                MAXWRK = MAX( MAXWRK, N*NRHS )
 304                MINWRK = 2*N + MAX( NRHS, M )
 305             END IF
 306             IF( N.GT.M ) THEN
 307                MINWRK = 2*M + MAX( NRHS, N )
 308                IF( N.GE.MNTHR ) THEN
 309 *
 310 *                 Path 2a - underdetermined, with many more columns
 311 *                 than rows
 312 *
 313 *                 Compute space needed for CGELQF
 314                   CALL CGELQF( M, N, A, LDA, DUM(1), DUM(1),
 315      $                -1, INFO )
 316                   LWORK_CGELQF=DUM(1)
 317 *                 Compute space needed for CGEBRD
 318                   CALL CGEBRD( M, M, A, LDA, S, S, DUM(1), DUM(1),
 319      $                         DUM(1), -1, INFO )
 320                   LWORK_CGEBRD=DUM(1)
 321 *                 Compute space needed for CUNMBR
 322                   CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA,
 323      $                DUM(1), B, LDB, DUM(1), -1, INFO )
 324                   LWORK_CUNMBR=DUM(1)
 325 *                 Compute space needed for CUNGBR
 326                   CALL CUNGBR( 'P', M, M, M, A, LDA, DUM(1),
 327      $                   DUM(1), -1, INFO )
 328                   LWORK_CUNGBR=DUM(1)
 329 *                 Compute space needed for CUNMLQ
 330                   CALL CUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, DUM(1),
 331      $                 B, LDB, DUM(1), -1, INFO )
 332                   LWORK_CUNMLQ=DUM(1)
 333 *                 Compute total workspace needed
 334                   MAXWRK = M + LWORK_CGELQF
 335                   MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_CGEBRD )
 336                   MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_CUNMBR )
 337                   MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_CUNGBR )
 338                   IF( NRHS.GT.1 ) THEN
 339                      MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
 340                   ELSE
 341                      MAXWRK = MAX( MAXWRK, M*M + 2*M )
 342                   END IF
 343                   MAXWRK = MAX( MAXWRK, M + LWORK_CUNMLQ )
 344                ELSE
 345 *
 346 *                 Path 2 - underdetermined
 347 *
 348 *                 Compute space needed for CGEBRD
 349                   CALL CGEBRD( M, N, A, LDA, S, S, DUM(1), DUM(1),
 350      $                         DUM(1), -1, INFO )
 351                   LWORK_CGEBRD=DUM(1)
 352 *                 Compute space needed for CUNMBR
 353                   CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, M, A, LDA,
 354      $                DUM(1), B, LDB, DUM(1), -1, INFO )
 355                   LWORK_CUNMBR=DUM(1)
 356 *                 Compute space needed for CUNGBR
 357                   CALL CUNGBR( 'P', M, N, M, A, LDA, DUM(1),
 358      $                   DUM(1), -1, INFO )
 359                   LWORK_CUNGBR=DUM(1)
 360                   MAXWRK = 2*M + LWORK_CGEBRD
 361                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNMBR )
 362                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_CUNGBR )
 363                   MAXWRK = MAX( MAXWRK, N*NRHS )
 364                END IF
 365             END IF
 366             MAXWRK = MAX( MINWRK, MAXWRK )
 367          END IF
 368          WORK( 1 ) = MAXWRK
 369 *
 370          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
 371      $      INFO = -12
 372       END IF
 373 *
 374       IF( INFO.NE.0 ) THEN
 375          CALL XERBLA( 'CGELSS', -INFO )
 376          RETURN
 377       ELSE IF( LQUERY ) THEN
 378          RETURN
 379       END IF
 380 *
 381 *     Quick return if possible
 382 *
 383       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
 384          RANK = 0
 385          RETURN
 386       END IF
 387 *
 388 *     Get machine parameters
 389 *
 390       EPS = SLAMCH( 'P' )
 391       SFMIN = SLAMCH( 'S' )
 392       SMLNUM = SFMIN / EPS
 393       BIGNUM = ONE / SMLNUM
 394       CALL SLABAD( SMLNUM, BIGNUM )
 395 *
 396 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 397 *
 398       ANRM = CLANGE( 'M', M, N, A, LDA, RWORK )
 399       IASCL = 0
 400       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 401 *
 402 *        Scale matrix norm up to SMLNUM
 403 *
 404          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
 405          IASCL = 1
 406       ELSE IF( ANRM.GT.BIGNUM ) THEN
 407 *
 408 *        Scale matrix norm down to BIGNUM
 409 *
 410          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
 411          IASCL = 2
 412       ELSE IF( ANRM.EQ.ZERO ) THEN
 413 *
 414 *        Matrix all zero. Return zero solution.
 415 *
 416          CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
 417          CALL SLASET( 'F', MINMN, 1, ZERO, ZERO, S, MINMN )
 418          RANK = 0
 419          GO TO 70
 420       END IF
 421 *
 422 *     Scale B if max element outside range [SMLNUM,BIGNUM]
 423 *
 424       BNRM = CLANGE( 'M', M, NRHS, B, LDB, RWORK )
 425       IBSCL = 0
 426       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
 427 *
 428 *        Scale matrix norm up to SMLNUM
 429 *
 430          CALL CLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
 431          IBSCL = 1
 432       ELSE IF( BNRM.GT.BIGNUM ) THEN
 433 *
 434 *        Scale matrix norm down to BIGNUM
 435 *
 436          CALL CLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
 437          IBSCL = 2
 438       END IF
 439 *
 440 *     Overdetermined case
 441 *
 442       IF( M.GE.N ) THEN
 443 *
 444 *        Path 1 - overdetermined or exactly determined
 445 *
 446          MM = M
 447          IF( M.GE.MNTHR ) THEN
 448 *
 449 *           Path 1a - overdetermined, with many more rows than columns
 450 *
 451             MM = N
 452             ITAU = 1
 453             IWORK = ITAU + N
 454 *
 455 *           Compute A=Q*R
 456 *           (CWorkspace: need 2*N, prefer N+N*NB)
 457 *           (RWorkspace: none)
 458 *
 459             CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
 460      $                   LWORK-IWORK+1, INFO )
 461 *
 462 *           Multiply B by transpose(Q)
 463 *           (CWorkspace: need N+NRHS, prefer N+NRHS*NB)
 464 *           (RWorkspace: none)
 465 *
 466             CALL CUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B,
 467      $                   LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
 468 *
 469 *           Zero out below R
 470 *
 471             IF( N.GT.1 )
 472      $         CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
 473      $                      LDA )
 474          END IF
 475 *
 476          IE = 1
 477          ITAUQ = 1
 478          ITAUP = ITAUQ + N
 479          IWORK = ITAUP + N
 480 *
 481 *        Bidiagonalize R in A
 482 *        (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
 483 *        (RWorkspace: need N)
 484 *
 485          CALL CGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
 486      $                WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
 487      $                INFO )
 488 *
 489 *        Multiply B by transpose of left bidiagonalizing vectors of R
 490 *        (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
 491 *        (RWorkspace: none)
 492 *
 493          CALL CUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
 494      $                B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
 495 *
 496 *        Generate right bidiagonalizing vectors of R in A
 497 *        (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
 498 *        (RWorkspace: none)
 499 *
 500          CALL CUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
 501      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
 502          IRWORK = IE + N
 503 *
 504 *        Perform bidiagonal QR iteration
 505 *          multiply B by transpose of left singular vectors
 506 *          compute right singular vectors in A
 507 *        (CWorkspace: none)
 508 *        (RWorkspace: need BDSPAC)
 509 *
 510          CALL CBDSQR( 'U', N, N, 0, NRHS, S, RWORK( IE ), A, LDA, DUM,
 511      $                1, B, LDB, RWORK( IRWORK ), INFO )
 512          IF( INFO.NE.0 )
 513      $      GO TO 70
 514 *
 515 *        Multiply B by reciprocals of singular values
 516 *
 517          THR = MAX( RCOND*S( 1 ), SFMIN )
 518          IF( RCOND.LT.ZERO )
 519      $      THR = MAX( EPS*S( 1 ), SFMIN )
 520          RANK = 0
 521          DO 10 I = 1, N
 522             IF( S( I ).GT.THR ) THEN
 523                CALL CSRSCL( NRHS, S( I ), B( I, 1 ), LDB )
 524                RANK = RANK + 1
 525             ELSE
 526                CALL CLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
 527             END IF
 528    10    CONTINUE
 529 *
 530 *        Multiply B by right singular vectors
 531 *        (CWorkspace: need N, prefer N*NRHS)
 532 *        (RWorkspace: none)
 533 *
 534          IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
 535             CALL CGEMM( 'C', 'N', N, NRHS, N, CONE, A, LDA, B, LDB,
 536      $                  CZERO, WORK, LDB )
 537             CALL CLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
 538          ELSE IF( NRHS.GT.1 ) THEN
 539             CHUNK = LWORK / N
 540             DO 20 I = 1, NRHS, CHUNK
 541                BL = MIN( NRHS-I+1, CHUNK )
 542                CALL CGEMM( 'C', 'N', N, BL, N, CONE, A, LDA, B( 1, I ),
 543      $                     LDB, CZERO, WORK, N )
 544                CALL CLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB )
 545    20       CONTINUE
 546          ELSE
 547             CALL CGEMV( 'C', N, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
 548             CALL CCOPY( N, WORK, 1, B, 1 )
 549          END IF
 550 *
 551       ELSE IF( N.GE.MNTHR .AND. LWORK.GE.3*M+M*M+MAX( M, NRHS, N-2*M ) )
 552      $          THEN
 553 *
 554 *        Underdetermined case, M much less than N
 555 *
 556 *        Path 2a - underdetermined, with many more columns than rows
 557 *        and sufficient workspace for an efficient algorithm
 558 *
 559          LDWORK = M
 560          IF( LWORK.GE.3*M+M*LDA+MAX( M, NRHS, N-2*M ) )
 561      $      LDWORK = LDA
 562          ITAU = 1
 563          IWORK = M + 1
 564 *
 565 *        Compute A=L*Q
 566 *        (CWorkspace: need 2*M, prefer M+M*NB)
 567 *        (RWorkspace: none)
 568 *
 569          CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
 570      $                LWORK-IWORK+1, INFO )
 571          IL = IWORK
 572 *
 573 *        Copy L to WORK(IL), zeroing out above it
 574 *
 575          CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
 576          CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
 577      $                LDWORK )
 578          IE = 1
 579          ITAUQ = IL + LDWORK*M
 580          ITAUP = ITAUQ + M
 581          IWORK = ITAUP + M
 582 *
 583 *        Bidiagonalize L in WORK(IL)
 584 *        (CWorkspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
 585 *        (RWorkspace: need M)
 586 *
 587          CALL CGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ),
 588      $                WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ),
 589      $                LWORK-IWORK+1, INFO )
 590 *
 591 *        Multiply B by transpose of left bidiagonalizing vectors of L
 592 *        (CWorkspace: need M*M+3*M+NRHS, prefer M*M+3*M+NRHS*NB)
 593 *        (RWorkspace: none)
 594 *
 595          CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK,
 596      $                WORK( ITAUQ ), B, LDB, WORK( IWORK ),
 597      $                LWORK-IWORK+1, INFO )
 598 *
 599 *        Generate right bidiagonalizing vectors of R in WORK(IL)
 600 *        (CWorkspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
 601 *        (RWorkspace: none)
 602 *
 603          CALL CUNGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ),
 604      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
 605          IRWORK = IE + M
 606 *
 607 *        Perform bidiagonal QR iteration, computing right singular
 608 *        vectors of L in WORK(IL) and multiplying B by transpose of
 609 *        left singular vectors
 610 *        (CWorkspace: need M*M)
 611 *        (RWorkspace: need BDSPAC)
 612 *
 613          CALL CBDSQR( 'U', M, M, 0, NRHS, S, RWORK( IE ), WORK( IL ),
 614      $                LDWORK, A, LDA, B, LDB, RWORK( IRWORK ), INFO )
 615          IF( INFO.NE.0 )
 616      $      GO TO 70
 617 *
 618 *        Multiply B by reciprocals of singular values
 619 *
 620          THR = MAX( RCOND*S( 1 ), SFMIN )
 621          IF( RCOND.LT.ZERO )
 622      $      THR = MAX( EPS*S( 1 ), SFMIN )
 623          RANK = 0
 624          DO 30 I = 1, M
 625             IF( S( I ).GT.THR ) THEN
 626                CALL CSRSCL( NRHS, S( I ), B( I, 1 ), LDB )
 627                RANK = RANK + 1
 628             ELSE
 629                CALL CLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
 630             END IF
 631    30    CONTINUE
 632          IWORK = IL + M*LDWORK
 633 *
 634 *        Multiply B by right singular vectors of L in WORK(IL)
 635 *        (CWorkspace: need M*M+2*M, prefer M*M+M+M*NRHS)
 636 *        (RWorkspace: none)
 637 *
 638          IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN
 639             CALL CGEMM( 'C', 'N', M, NRHS, M, CONE, WORK( IL ), LDWORK,
 640      $                  B, LDB, CZERO, WORK( IWORK ), LDB )
 641             CALL CLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB )
 642          ELSE IF( NRHS.GT.1 ) THEN
 643             CHUNK = ( LWORK-IWORK+1 ) / M
 644             DO 40 I = 1, NRHS, CHUNK
 645                BL = MIN( NRHS-I+1, CHUNK )
 646                CALL CGEMM( 'C', 'N', M, BL, M, CONE, WORK( IL ), LDWORK,
 647      $                     B( 1, I ), LDB, CZERO, WORK( IWORK ), M )
 648                CALL CLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ),
 649      $                      LDB )
 650    40       CONTINUE
 651          ELSE
 652             CALL CGEMV( 'C', M, M, CONE, WORK( IL ), LDWORK, B( 1, 1 ),
 653      $                  1, CZERO, WORK( IWORK ), 1 )
 654             CALL CCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 )
 655          END IF
 656 *
 657 *        Zero out below first M rows of B
 658 *
 659          CALL CLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
 660          IWORK = ITAU + M
 661 *
 662 *        Multiply transpose(Q) by B
 663 *        (CWorkspace: need M+NRHS, prefer M+NHRS*NB)
 664 *        (RWorkspace: none)
 665 *
 666          CALL CUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B,
 667      $                LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
 668 *
 669       ELSE
 670 *
 671 *        Path 2 - remaining underdetermined cases
 672 *
 673          IE = 1
 674          ITAUQ = 1
 675          ITAUP = ITAUQ + M
 676          IWORK = ITAUP + M
 677 *
 678 *        Bidiagonalize A
 679 *        (CWorkspace: need 3*M, prefer 2*M+(M+N)*NB)
 680 *        (RWorkspace: need N)
 681 *
 682          CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
 683      $                WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
 684      $                INFO )
 685 *
 686 *        Multiply B by transpose of left bidiagonalizing vectors
 687 *        (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
 688 *        (RWorkspace: none)
 689 *
 690          CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ),
 691      $                B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
 692 *
 693 *        Generate right bidiagonalizing vectors in A
 694 *        (CWorkspace: need 3*M, prefer 2*M+M*NB)
 695 *        (RWorkspace: none)
 696 *
 697          CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
 698      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
 699          IRWORK = IE + M
 700 *
 701 *        Perform bidiagonal QR iteration,
 702 *           computing right singular vectors of A in A and
 703 *           multiplying B by transpose of left singular vectors
 704 *        (CWorkspace: none)
 705 *        (RWorkspace: need BDSPAC)
 706 *
 707          CALL CBDSQR( 'L', M, N, 0, NRHS, S, RWORK( IE ), A, LDA, DUM,
 708      $                1, B, LDB, RWORK( IRWORK ), INFO )
 709          IF( INFO.NE.0 )
 710      $      GO TO 70
 711 *
 712 *        Multiply B by reciprocals of singular values
 713 *
 714          THR = MAX( RCOND*S( 1 ), SFMIN )
 715          IF( RCOND.LT.ZERO )
 716      $      THR = MAX( EPS*S( 1 ), SFMIN )
 717          RANK = 0
 718          DO 50 I = 1, M
 719             IF( S( I ).GT.THR ) THEN
 720                CALL CSRSCL( NRHS, S( I ), B( I, 1 ), LDB )
 721                RANK = RANK + 1
 722             ELSE
 723                CALL CLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
 724             END IF
 725    50    CONTINUE
 726 *
 727 *        Multiply B by right singular vectors of A
 728 *        (CWorkspace: need N, prefer N*NRHS)
 729 *        (RWorkspace: none)
 730 *
 731          IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
 732             CALL CGEMM( 'C', 'N', N, NRHS, M, CONE, A, LDA, B, LDB,
 733      $                  CZERO, WORK, LDB )
 734             CALL CLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
 735          ELSE IF( NRHS.GT.1 ) THEN
 736             CHUNK = LWORK / N
 737             DO 60 I = 1, NRHS, CHUNK
 738                BL = MIN( NRHS-I+1, CHUNK )
 739                CALL CGEMM( 'C', 'N', N, BL, M, CONE, A, LDA, B( 1, I ),
 740      $                     LDB, CZERO, WORK, N )
 741                CALL CLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB )
 742    60       CONTINUE
 743          ELSE
 744             CALL CGEMV( 'C', M, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
 745             CALL CCOPY( N, WORK, 1, B, 1 )
 746          END IF
 747       END IF
 748 *
 749 *     Undo scaling
 750 *
 751       IF( IASCL.EQ.1 ) THEN
 752          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
 753          CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
 754      $                INFO )
 755       ELSE IF( IASCL.EQ.2 ) THEN
 756          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
 757          CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
 758      $                INFO )
 759       END IF
 760       IF( IBSCL.EQ.1 ) THEN
 761          CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
 762       ELSE IF( IBSCL.EQ.2 ) THEN
 763          CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
 764       END IF
 765    70 CONTINUE
 766       WORK( 1 ) = MAXWRK
 767       RETURN
 768 *
 769 *     End of CGELSS
 770 *
 771       END