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