001:       SUBROUTINE CGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
002:      $                   WORK, LWORK, RWORK, INFO )
003: *
004: *  -- LAPACK driver routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
010:       REAL               RCOND
011: *     ..
012: *     .. Array Arguments ..
013:       REAL               RWORK( * ), S( * )
014:       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  CGELSS computes the minimum norm solution to a complex linear
021: *  least squares problem:
022: *
023: *  Minimize 2-norm(| b - A*x |).
024: *
025: *  using the singular value decomposition (SVD) of A. A is an M-by-N
026: *  matrix which may be rank-deficient.
027: *
028: *  Several right hand side vectors b and solution vectors x can be
029: *  handled in a single call; they are stored as the columns of the
030: *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix
031: *  X.
032: *
033: *  The effective rank of A is determined by treating as zero those
034: *  singular values which are less than RCOND times the largest singular
035: *  value.
036: *
037: *  Arguments
038: *  =========
039: *
040: *  M       (input) INTEGER
041: *          The number of rows of the matrix A. M >= 0.
042: *
043: *  N       (input) INTEGER
044: *          The number of columns of the matrix A. N >= 0.
045: *
046: *  NRHS    (input) INTEGER
047: *          The number of right hand sides, i.e., the number of columns
048: *          of the matrices B and X. NRHS >= 0.
049: *
050: *  A       (input/output) COMPLEX array, dimension (LDA,N)
051: *          On entry, the M-by-N matrix A.
052: *          On exit, the first min(m,n) rows of A are overwritten with
053: *          its right singular vectors, stored rowwise.
054: *
055: *  LDA     (input) INTEGER
056: *          The leading dimension of the array A. LDA >= max(1,M).
057: *
058: *  B       (input/output) COMPLEX array, dimension (LDB,NRHS)
059: *          On entry, the M-by-NRHS right hand side matrix B.
060: *          On exit, B is overwritten by the N-by-NRHS solution matrix X.
061: *          If m >= n and RANK = n, the residual sum-of-squares for
062: *          the solution in the i-th column is given by the sum of
063: *          squares of the modulus of elements n+1:m in that column.
064: *
065: *  LDB     (input) INTEGER
066: *          The leading dimension of the array B.  LDB >= max(1,M,N).
067: *
068: *  S       (output) REAL array, dimension (min(M,N))
069: *          The singular values of A in decreasing order.
070: *          The condition number of A in the 2-norm = S(1)/S(min(m,n)).
071: *
072: *  RCOND   (input) REAL
073: *          RCOND is used to determine the effective rank of A.
074: *          Singular values S(i) <= RCOND*S(1) are treated as zero.
075: *          If RCOND < 0, machine precision is used instead.
076: *
077: *  RANK    (output) INTEGER
078: *          The effective rank of A, i.e., the number of singular values
079: *          which are greater than RCOND*S(1).
080: *
081: *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
082: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
083: *
084: *  LWORK   (input) INTEGER
085: *          The dimension of the array WORK. LWORK >= 1, and also:
086: *          LWORK >=  2*min(M,N) + max(M,N,NRHS)
087: *          For good performance, LWORK should generally be larger.
088: *
089: *          If LWORK = -1, then a workspace query is assumed; the routine
090: *          only calculates the optimal size of the WORK array, returns
091: *          this value as the first entry of the WORK array, and no error
092: *          message related to LWORK is issued by XERBLA.
093: *
094: *  RWORK   (workspace) REAL array, dimension (5*min(M,N))
095: *
096: *  INFO    (output) INTEGER
097: *          = 0:  successful exit
098: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
099: *          > 0:  the algorithm for computing the SVD failed to converge;
100: *                if INFO = i, i off-diagonal elements of an intermediate
101: *                bidiagonal form did not converge to zero.
102: *
103: *  =====================================================================
104: *
105: *     .. Parameters ..
106:       REAL               ZERO, ONE
107:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
108:       COMPLEX            CZERO, CONE
109:       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
110:      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
111: *     ..
112: *     .. Local Scalars ..
113:       LOGICAL            LQUERY
114:       INTEGER            BL, CHUNK, I, IASCL, IBSCL, IE, IL, IRWORK,
115:      $                   ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN,
116:      $                   MAXWRK, MINMN, MINWRK, MM, MNTHR
117:       REAL               ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
118: *     ..
119: *     .. Local Arrays ..
120:       COMPLEX            VDUM( 1 )
121: *     ..
122: *     .. External Subroutines ..
123:       EXTERNAL           CBDSQR, CCOPY, CGEBRD, CGELQF, CGEMM, CGEMV,
124:      $                   CGEQRF, CLACPY, CLASCL, CLASET, CSRSCL, CUNGBR,
125:      $                   CUNMBR, CUNMLQ, CUNMQR, SLABAD, SLASCL, SLASET,
126:      $                   XERBLA
127: *     ..
128: *     .. External Functions ..
129:       INTEGER            ILAENV
130:       REAL               CLANGE, SLAMCH
131:       EXTERNAL           ILAENV, CLANGE, SLAMCH
132: *     ..
133: *     .. Intrinsic Functions ..
134:       INTRINSIC          MAX, MIN
135: *     ..
136: *     .. Executable Statements ..
137: *
138: *     Test the input arguments
139: *
140:       INFO = 0
141:       MINMN = MIN( M, N )
142:       MAXMN = MAX( M, N )
143:       LQUERY = ( LWORK.EQ.-1 )
144:       IF( M.LT.0 ) THEN
145:          INFO = -1
146:       ELSE IF( N.LT.0 ) THEN
147:          INFO = -2
148:       ELSE IF( NRHS.LT.0 ) THEN
149:          INFO = -3
150:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
151:          INFO = -5
152:       ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
153:          INFO = -7
154:       END IF
155: *
156: *     Compute workspace
157: *      (Note: Comments in the code beginning "Workspace:" describe the
158: *       minimal amount of workspace needed at that point in the code,
159: *       as well as the preferred amount for good performance.
160: *       CWorkspace refers to complex workspace, and RWorkspace refers
161: *       to real workspace. NB refers to the optimal block size for the
162: *       immediately following subroutine, as returned by ILAENV.)
163: *
164:       IF( INFO.EQ.0 ) THEN
165:          MINWRK = 1
166:          MAXWRK = 1
167:          IF( MINMN.GT.0 ) THEN
168:             MM = M
169:             MNTHR = ILAENV( 6, 'CGELSS', ' ', M, N, NRHS, -1 )
170:             IF( M.GE.N .AND. M.GE.MNTHR ) THEN
171: *
172: *              Path 1a - overdetermined, with many more rows than
173: *                        columns
174: *
175:                MM = N
176:                MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'CGEQRF', ' ', M,
177:      $                       N, -1, -1 ) )
178:                MAXWRK = MAX( MAXWRK, N + NRHS*ILAENV( 1, 'CUNMQR', 'LC',
179:      $                       M, NRHS, N, -1 ) )
180:             END IF
181:             IF( M.GE.N ) THEN
182: *
183: *              Path 1 - overdetermined or exactly determined
184: *
185:                MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1,
186:      $                       'CGEBRD', ' ', MM, N, -1, -1 ) )
187:                MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'CUNMBR',
188:      $                       'QLC', MM, NRHS, N, -1 ) )
189:                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
190:      $                       'CUNGBR', 'P', N, N, N, -1 ) )
191:                MAXWRK = MAX( MAXWRK, N*NRHS )
192:                MINWRK = 2*N + MAX( NRHS, M )
193:             END IF
194:             IF( N.GT.M ) THEN
195:                MINWRK = 2*M + MAX( NRHS, N )
196:                IF( N.GE.MNTHR ) THEN
197: *
198: *                 Path 2a - underdetermined, with many more columns
199: *                 than rows
200: *
201:                   MAXWRK = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1,
202:      $                     -1 )
203:                   MAXWRK = MAX( MAXWRK, 3*M + M*M + 2*M*ILAENV( 1,
204:      $                          'CGEBRD', ' ', M, M, -1, -1 ) )
205:                   MAXWRK = MAX( MAXWRK, 3*M + M*M + NRHS*ILAENV( 1,
206:      $                          'CUNMBR', 'QLC', M, NRHS, M, -1 ) )
207:                   MAXWRK = MAX( MAXWRK, 3*M + M*M + ( M - 1 )*ILAENV( 1,
208:      $                          'CUNGBR', 'P', M, M, M, -1 ) )
209:                   IF( NRHS.GT.1 ) THEN
210:                      MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
211:                   ELSE
212:                      MAXWRK = MAX( MAXWRK, M*M + 2*M )
213:                   END IF
214:                   MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'CUNMLQ',
215:      $                          'LC', N, NRHS, M, -1 ) )
216:                ELSE
217: *
218: *                 Path 2 - underdetermined
219: *
220:                   MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'CGEBRD', ' ', M,
221:      $                     N, -1, -1 )
222:                   MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'CUNMBR',
223:      $                          'QLC', M, NRHS, M, -1 ) )
224:                   MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'CUNGBR',
225:      $                          'P', M, N, M, -1 ) )
226:                   MAXWRK = MAX( MAXWRK, N*NRHS )
227:                END IF
228:             END IF
229:             MAXWRK = MAX( MINWRK, MAXWRK )
230:          END IF
231:          WORK( 1 ) = MAXWRK
232: *
233:          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
234:      $      INFO = -12
235:       END IF
236: *
237:       IF( INFO.NE.0 ) THEN
238:          CALL XERBLA( 'CGELSS', -INFO )
239:          RETURN
240:       ELSE IF( LQUERY ) THEN
241:          RETURN
242:       END IF
243: *
244: *     Quick return if possible
245: *
246:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
247:          RANK = 0
248:          RETURN
249:       END IF
250: *
251: *     Get machine parameters
252: *
253:       EPS = SLAMCH( 'P' )
254:       SFMIN = SLAMCH( 'S' )
255:       SMLNUM = SFMIN / EPS
256:       BIGNUM = ONE / SMLNUM
257:       CALL SLABAD( SMLNUM, BIGNUM )
258: *
259: *     Scale A if max element outside range [SMLNUM,BIGNUM]
260: *
261:       ANRM = CLANGE( 'M', M, N, A, LDA, RWORK )
262:       IASCL = 0
263:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
264: *
265: *        Scale matrix norm up to SMLNUM
266: *
267:          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
268:          IASCL = 1
269:       ELSE IF( ANRM.GT.BIGNUM ) THEN
270: *
271: *        Scale matrix norm down to BIGNUM
272: *
273:          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
274:          IASCL = 2
275:       ELSE IF( ANRM.EQ.ZERO ) THEN
276: *
277: *        Matrix all zero. Return zero solution.
278: *
279:          CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
280:          CALL SLASET( 'F', MINMN, 1, ZERO, ZERO, S, MINMN )
281:          RANK = 0
282:          GO TO 70
283:       END IF
284: *
285: *     Scale B if max element outside range [SMLNUM,BIGNUM]
286: *
287:       BNRM = CLANGE( 'M', M, NRHS, B, LDB, RWORK )
288:       IBSCL = 0
289:       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
290: *
291: *        Scale matrix norm up to SMLNUM
292: *
293:          CALL CLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
294:          IBSCL = 1
295:       ELSE IF( BNRM.GT.BIGNUM ) THEN
296: *
297: *        Scale matrix norm down to BIGNUM
298: *
299:          CALL CLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
300:          IBSCL = 2
301:       END IF
302: *
303: *     Overdetermined case
304: *
305:       IF( M.GE.N ) THEN
306: *
307: *        Path 1 - overdetermined or exactly determined
308: *
309:          MM = M
310:          IF( M.GE.MNTHR ) THEN
311: *
312: *           Path 1a - overdetermined, with many more rows than columns
313: *
314:             MM = N
315:             ITAU = 1
316:             IWORK = ITAU + N
317: *
318: *           Compute A=Q*R
319: *           (CWorkspace: need 2*N, prefer N+N*NB)
320: *           (RWorkspace: none)
321: *
322:             CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
323:      $                   LWORK-IWORK+1, INFO )
324: *
325: *           Multiply B by transpose(Q)
326: *           (CWorkspace: need N+NRHS, prefer N+NRHS*NB)
327: *           (RWorkspace: none)
328: *
329:             CALL CUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B,
330:      $                   LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
331: *
332: *           Zero out below R
333: *
334:             IF( N.GT.1 )
335:      $         CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
336:      $                      LDA )
337:          END IF
338: *
339:          IE = 1
340:          ITAUQ = 1
341:          ITAUP = ITAUQ + N
342:          IWORK = ITAUP + N
343: *
344: *        Bidiagonalize R in A
345: *        (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
346: *        (RWorkspace: need N)
347: *
348:          CALL CGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
349:      $                WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
350:      $                INFO )
351: *
352: *        Multiply B by transpose of left bidiagonalizing vectors of R
353: *        (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
354: *        (RWorkspace: none)
355: *
356:          CALL CUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
357:      $                B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
358: *
359: *        Generate right bidiagonalizing vectors of R in A
360: *        (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
361: *        (RWorkspace: none)
362: *
363:          CALL CUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
364:      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
365:          IRWORK = IE + N
366: *
367: *        Perform bidiagonal QR iteration
368: *          multiply B by transpose of left singular vectors
369: *          compute right singular vectors in A
370: *        (CWorkspace: none)
371: *        (RWorkspace: need BDSPAC)
372: *
373:          CALL CBDSQR( 'U', N, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM,
374:      $                1, B, LDB, RWORK( IRWORK ), INFO )
375:          IF( INFO.NE.0 )
376:      $      GO TO 70
377: *
378: *        Multiply B by reciprocals of singular values
379: *
380:          THR = MAX( RCOND*S( 1 ), SFMIN )
381:          IF( RCOND.LT.ZERO )
382:      $      THR = MAX( EPS*S( 1 ), SFMIN )
383:          RANK = 0
384:          DO 10 I = 1, N
385:             IF( S( I ).GT.THR ) THEN
386:                CALL CSRSCL( NRHS, S( I ), B( I, 1 ), LDB )
387:                RANK = RANK + 1
388:             ELSE
389:                CALL CLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
390:             END IF
391:    10    CONTINUE
392: *
393: *        Multiply B by right singular vectors
394: *        (CWorkspace: need N, prefer N*NRHS)
395: *        (RWorkspace: none)
396: *
397:          IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
398:             CALL CGEMM( 'C', 'N', N, NRHS, N, CONE, A, LDA, B, LDB,
399:      $                  CZERO, WORK, LDB )
400:             CALL CLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
401:          ELSE IF( NRHS.GT.1 ) THEN
402:             CHUNK = LWORK / N
403:             DO 20 I = 1, NRHS, CHUNK
404:                BL = MIN( NRHS-I+1, CHUNK )
405:                CALL CGEMM( 'C', 'N', N, BL, N, CONE, A, LDA, B( 1, I ),
406:      $                     LDB, CZERO, WORK, N )
407:                CALL CLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB )
408:    20       CONTINUE
409:          ELSE
410:             CALL CGEMV( 'C', N, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
411:             CALL CCOPY( N, WORK, 1, B, 1 )
412:          END IF
413: *
414:       ELSE IF( N.GE.MNTHR .AND. LWORK.GE.3*M+M*M+MAX( M, NRHS, N-2*M ) )
415:      $          THEN
416: *
417: *        Underdetermined case, M much less than N
418: *
419: *        Path 2a - underdetermined, with many more columns than rows
420: *        and sufficient workspace for an efficient algorithm
421: *
422:          LDWORK = M
423:          IF( LWORK.GE.3*M+M*LDA+MAX( M, NRHS, N-2*M ) )
424:      $      LDWORK = LDA
425:          ITAU = 1
426:          IWORK = M + 1
427: *
428: *        Compute A=L*Q
429: *        (CWorkspace: need 2*M, prefer M+M*NB)
430: *        (RWorkspace: none)
431: *
432:          CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
433:      $                LWORK-IWORK+1, INFO )
434:          IL = IWORK
435: *
436: *        Copy L to WORK(IL), zeroing out above it
437: *
438:          CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
439:          CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
440:      $                LDWORK )
441:          IE = 1
442:          ITAUQ = IL + LDWORK*M
443:          ITAUP = ITAUQ + M
444:          IWORK = ITAUP + M
445: *
446: *        Bidiagonalize L in WORK(IL)
447: *        (CWorkspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
448: *        (RWorkspace: need M)
449: *
450:          CALL CGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ),
451:      $                WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ),
452:      $                LWORK-IWORK+1, INFO )
453: *
454: *        Multiply B by transpose of left bidiagonalizing vectors of L
455: *        (CWorkspace: need M*M+3*M+NRHS, prefer M*M+3*M+NRHS*NB)
456: *        (RWorkspace: none)
457: *
458:          CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK,
459:      $                WORK( ITAUQ ), B, LDB, WORK( IWORK ),
460:      $                LWORK-IWORK+1, INFO )
461: *
462: *        Generate right bidiagonalizing vectors of R in WORK(IL)
463: *        (CWorkspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
464: *        (RWorkspace: none)
465: *
466:          CALL CUNGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ),
467:      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
468:          IRWORK = IE + M
469: *
470: *        Perform bidiagonal QR iteration, computing right singular
471: *        vectors of L in WORK(IL) and multiplying B by transpose of
472: *        left singular vectors
473: *        (CWorkspace: need M*M)
474: *        (RWorkspace: need BDSPAC)
475: *
476:          CALL CBDSQR( 'U', M, M, 0, NRHS, S, RWORK( IE ), WORK( IL ),
477:      $                LDWORK, A, LDA, B, LDB, RWORK( IRWORK ), INFO )
478:          IF( INFO.NE.0 )
479:      $      GO TO 70
480: *
481: *        Multiply B by reciprocals of singular values
482: *
483:          THR = MAX( RCOND*S( 1 ), SFMIN )
484:          IF( RCOND.LT.ZERO )
485:      $      THR = MAX( EPS*S( 1 ), SFMIN )
486:          RANK = 0
487:          DO 30 I = 1, M
488:             IF( S( I ).GT.THR ) THEN
489:                CALL CSRSCL( NRHS, S( I ), B( I, 1 ), LDB )
490:                RANK = RANK + 1
491:             ELSE
492:                CALL CLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
493:             END IF
494:    30    CONTINUE
495:          IWORK = IL + M*LDWORK
496: *
497: *        Multiply B by right singular vectors of L in WORK(IL)
498: *        (CWorkspace: need M*M+2*M, prefer M*M+M+M*NRHS)
499: *        (RWorkspace: none)
500: *
501:          IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN
502:             CALL CGEMM( 'C', 'N', M, NRHS, M, CONE, WORK( IL ), LDWORK,
503:      $                  B, LDB, CZERO, WORK( IWORK ), LDB )
504:             CALL CLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB )
505:          ELSE IF( NRHS.GT.1 ) THEN
506:             CHUNK = ( LWORK-IWORK+1 ) / M
507:             DO 40 I = 1, NRHS, CHUNK
508:                BL = MIN( NRHS-I+1, CHUNK )
509:                CALL CGEMM( 'C', 'N', M, BL, M, CONE, WORK( IL ), LDWORK,
510:      $                     B( 1, I ), LDB, CZERO, WORK( IWORK ), M )
511:                CALL CLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ),
512:      $                      LDB )
513:    40       CONTINUE
514:          ELSE
515:             CALL CGEMV( 'C', M, M, CONE, WORK( IL ), LDWORK, B( 1, 1 ),
516:      $                  1, CZERO, WORK( IWORK ), 1 )
517:             CALL CCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 )
518:          END IF
519: *
520: *        Zero out below first M rows of B
521: *
522:          CALL CLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
523:          IWORK = ITAU + M
524: *
525: *        Multiply transpose(Q) by B
526: *        (CWorkspace: need M+NRHS, prefer M+NHRS*NB)
527: *        (RWorkspace: none)
528: *
529:          CALL CUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B,
530:      $                LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
531: *
532:       ELSE
533: *
534: *        Path 2 - remaining underdetermined cases
535: *
536:          IE = 1
537:          ITAUQ = 1
538:          ITAUP = ITAUQ + M
539:          IWORK = ITAUP + M
540: *
541: *        Bidiagonalize A
542: *        (CWorkspace: need 3*M, prefer 2*M+(M+N)*NB)
543: *        (RWorkspace: need N)
544: *
545:          CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
546:      $                WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
547:      $                INFO )
548: *
549: *        Multiply B by transpose of left bidiagonalizing vectors
550: *        (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
551: *        (RWorkspace: none)
552: *
553:          CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ),
554:      $                B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
555: *
556: *        Generate right bidiagonalizing vectors in A
557: *        (CWorkspace: need 3*M, prefer 2*M+M*NB)
558: *        (RWorkspace: none)
559: *
560:          CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
561:      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
562:          IRWORK = IE + M
563: *
564: *        Perform bidiagonal QR iteration,
565: *           computing right singular vectors of A in A and
566: *           multiplying B by transpose of left singular vectors
567: *        (CWorkspace: none)
568: *        (RWorkspace: need BDSPAC)
569: *
570:          CALL CBDSQR( 'L', M, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM,
571:      $                1, B, LDB, RWORK( IRWORK ), INFO )
572:          IF( INFO.NE.0 )
573:      $      GO TO 70
574: *
575: *        Multiply B by reciprocals of singular values
576: *
577:          THR = MAX( RCOND*S( 1 ), SFMIN )
578:          IF( RCOND.LT.ZERO )
579:      $      THR = MAX( EPS*S( 1 ), SFMIN )
580:          RANK = 0
581:          DO 50 I = 1, M
582:             IF( S( I ).GT.THR ) THEN
583:                CALL CSRSCL( NRHS, S( I ), B( I, 1 ), LDB )
584:                RANK = RANK + 1
585:             ELSE
586:                CALL CLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
587:             END IF
588:    50    CONTINUE
589: *
590: *        Multiply B by right singular vectors of A
591: *        (CWorkspace: need N, prefer N*NRHS)
592: *        (RWorkspace: none)
593: *
594:          IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
595:             CALL CGEMM( 'C', 'N', N, NRHS, M, CONE, A, LDA, B, LDB,
596:      $                  CZERO, WORK, LDB )
597:             CALL CLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
598:          ELSE IF( NRHS.GT.1 ) THEN
599:             CHUNK = LWORK / N
600:             DO 60 I = 1, NRHS, CHUNK
601:                BL = MIN( NRHS-I+1, CHUNK )
602:                CALL CGEMM( 'C', 'N', N, BL, M, CONE, A, LDA, B( 1, I ),
603:      $                     LDB, CZERO, WORK, N )
604:                CALL CLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB )
605:    60       CONTINUE
606:          ELSE
607:             CALL CGEMV( 'C', M, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
608:             CALL CCOPY( N, WORK, 1, B, 1 )
609:          END IF
610:       END IF
611: *
612: *     Undo scaling
613: *
614:       IF( IASCL.EQ.1 ) THEN
615:          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
616:          CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
617:      $                INFO )
618:       ELSE IF( IASCL.EQ.2 ) THEN
619:          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
620:          CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
621:      $                INFO )
622:       END IF
623:       IF( IBSCL.EQ.1 ) THEN
624:          CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
625:       ELSE IF( IBSCL.EQ.2 ) THEN
626:          CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
627:       END IF
628:    70 CONTINUE
629:       WORK( 1 ) = MAXWRK
630:       RETURN
631: *
632: *     End of CGELSS
633: *
634:       END
635: