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