LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dgelsd.f
Go to the documentation of this file.
1 *> \brief <b> DGELSD computes the minimum-norm solution to a linear least squares problem 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 DGELSD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelsd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelsd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelsd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
22 * WORK, LWORK, IWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
26 * DOUBLE PRECISION RCOND
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IWORK( * )
30 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> DGELSD computes the minimum-norm solution to a real linear least
40 *> squares problem:
41 *> minimize 2-norm(| b - A*x |)
42 *> using the singular value decomposition (SVD) of A. A is an M-by-N
43 *> matrix which may be rank-deficient.
44 *>
45 *> Several right hand side vectors b and solution vectors x can be
46 *> handled in a single call; they are stored as the columns of the
47 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
48 *> matrix X.
49 *>
50 *> The problem is solved in three steps:
51 *> (1) Reduce the coefficient matrix A to bidiagonal form with
52 *> Householder transformations, reducing the original problem
53 *> into a "bidiagonal least squares problem" (BLS)
54 *> (2) Solve the BLS using a divide and conquer approach.
55 *> (3) Apply back all the Householder tranformations to solve
56 *> the original least squares problem.
57 *>
58 *> The effective rank of A is determined by treating as zero those
59 *> singular values which are less than RCOND times the largest singular
60 *> value.
61 *>
62 *> The divide and conquer algorithm makes very mild assumptions about
63 *> floating point arithmetic. It will work on machines with a guard
64 *> digit in add/subtract, or on those binary machines without guard
65 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
66 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
67 *> without guard digits, but we know of none.
68 *> \endverbatim
69 *
70 * Arguments:
71 * ==========
72 *
73 *> \param[in] M
74 *> \verbatim
75 *> M is INTEGER
76 *> The number of rows of A. M >= 0.
77 *> \endverbatim
78 *>
79 *> \param[in] N
80 *> \verbatim
81 *> N is INTEGER
82 *> The number of columns of A. N >= 0.
83 *> \endverbatim
84 *>
85 *> \param[in] NRHS
86 *> \verbatim
87 *> NRHS is INTEGER
88 *> The number of right hand sides, i.e., the number of columns
89 *> of the matrices B and X. NRHS >= 0.
90 *> \endverbatim
91 *>
92 *> \param[in] A
93 *> \verbatim
94 *> A is DOUBLE PRECISION array, dimension (LDA,N)
95 *> On entry, the M-by-N matrix A.
96 *> On exit, A has been destroyed.
97 *> \endverbatim
98 *>
99 *> \param[in] LDA
100 *> \verbatim
101 *> LDA is INTEGER
102 *> The leading dimension of the array A. LDA >= max(1,M).
103 *> \endverbatim
104 *>
105 *> \param[in,out] B
106 *> \verbatim
107 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
108 *> On entry, the M-by-NRHS right hand side matrix B.
109 *> On exit, B is overwritten by the N-by-NRHS solution
110 *> matrix X. If m >= n and RANK = n, the residual
111 *> sum-of-squares for the solution in the i-th column is given
112 *> by the sum of squares of elements n+1:m in that column.
113 *> \endverbatim
114 *>
115 *> \param[in] LDB
116 *> \verbatim
117 *> LDB is INTEGER
118 *> The leading dimension of the array B. LDB >= max(1,max(M,N)).
119 *> \endverbatim
120 *>
121 *> \param[out] S
122 *> \verbatim
123 *> S is DOUBLE PRECISION array, dimension (min(M,N))
124 *> The singular values of A in decreasing order.
125 *> The condition number of A in the 2-norm = S(1)/S(min(m,n)).
126 *> \endverbatim
127 *>
128 *> \param[in] RCOND
129 *> \verbatim
130 *> RCOND is DOUBLE PRECISION
131 *> RCOND is used to determine the effective rank of A.
132 *> Singular values S(i) <= RCOND*S(1) are treated as zero.
133 *> If RCOND < 0, machine precision is used instead.
134 *> \endverbatim
135 *>
136 *> \param[out] RANK
137 *> \verbatim
138 *> RANK is INTEGER
139 *> The effective rank of A, i.e., the number of singular values
140 *> which are greater than RCOND*S(1).
141 *> \endverbatim
142 *>
143 *> \param[out] WORK
144 *> \verbatim
145 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
146 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
147 *> \endverbatim
148 *>
149 *> \param[in] LWORK
150 *> \verbatim
151 *> LWORK is INTEGER
152 *> The dimension of the array WORK. LWORK must be at least 1.
153 *> The exact minimum amount of workspace needed depends on M,
154 *> N and NRHS. As long as LWORK is at least
155 *> 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2,
156 *> if M is greater than or equal to N or
157 *> 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2,
158 *> if M is less than N, the code will execute correctly.
159 *> SMLSIZ is returned by ILAENV and is equal to the maximum
160 *> size of the subproblems at the bottom of the computation
161 *> tree (usually about 25), and
162 *> NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
163 *> For good performance, LWORK should generally be larger.
164 *>
165 *> If LWORK = -1, then a workspace query is assumed; the routine
166 *> only calculates the optimal size of the WORK array, returns
167 *> this value as the first entry of the WORK array, and no error
168 *> message related to LWORK is issued by XERBLA.
169 *> \endverbatim
170 *>
171 *> \param[out] IWORK
172 *> \verbatim
173 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
174 *> LIWORK >= max(1, 3 * MINMN * NLVL + 11 * MINMN),
175 *> where MINMN = MIN( M,N ).
176 *> On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
177 *> \endverbatim
178 *>
179 *> \param[out] INFO
180 *> \verbatim
181 *> INFO is INTEGER
182 *> = 0: successful exit
183 *> < 0: if INFO = -i, the i-th argument had an illegal value.
184 *> > 0: the algorithm for computing the SVD failed to converge;
185 *> if INFO = i, i off-diagonal elements of an intermediate
186 *> bidiagonal form did not converge to zero.
187 *> \endverbatim
188 *
189 * Authors:
190 * ========
191 *
192 *> \author Univ. of Tennessee
193 *> \author Univ. of California Berkeley
194 *> \author Univ. of Colorado Denver
195 *> \author NAG Ltd.
196 *
197 *> \date November 2011
198 *
199 *> \ingroup doubleGEsolve
200 *
201 *> \par Contributors:
202 * ==================
203 *>
204 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
205 *> California at Berkeley, USA \n
206 *> Osni Marques, LBNL/NERSC, USA \n
207 *
208 * =====================================================================
209  SUBROUTINE dgelsd( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
210  $ work, lwork, iwork, info )
211 *
212 * -- LAPACK driver routine (version 3.4.0) --
213 * -- LAPACK is a software package provided by Univ. of Tennessee, --
214 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
215 * November 2011
216 *
217 * .. Scalar Arguments ..
218  INTEGER info, lda, ldb, lwork, m, n, nrhs, rank
219  DOUBLE PRECISION rcond
220 * ..
221 * .. Array Arguments ..
222  INTEGER iwork( * )
223  DOUBLE PRECISION a( lda, * ), b( ldb, * ), s( * ), work( * )
224 * ..
225 *
226 * =====================================================================
227 *
228 * .. Parameters ..
229  DOUBLE PRECISION zero, one, two
230  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
231 * ..
232 * .. Local Scalars ..
233  LOGICAL lquery
234  INTEGER iascl, ibscl, ie, il, itau, itaup, itauq,
235  $ ldwork, liwork, maxmn, maxwrk, minmn, minwrk,
236  $ mm, mnthr, nlvl, nwork, smlsiz, wlalsd
237  DOUBLE PRECISION anrm, bignum, bnrm, eps, sfmin, smlnum
238 * ..
239 * .. External Subroutines ..
240  EXTERNAL dgebrd, dgelqf, dgeqrf, dlabad, dlacpy, dlalsd,
242 * ..
243 * .. External Functions ..
244  INTEGER ilaenv
245  DOUBLE PRECISION dlamch, dlange
246  EXTERNAL ilaenv, dlamch, dlange
247 * ..
248 * .. Intrinsic Functions ..
249  INTRINSIC dble, int, log, max, min
250 * ..
251 * .. Executable Statements ..
252 *
253 * Test the input arguments.
254 *
255  info = 0
256  minmn = min( m, n )
257  maxmn = max( m, n )
258  mnthr = ilaenv( 6, 'DGELSD', ' ', m, n, nrhs, -1 )
259  lquery = ( lwork.EQ.-1 )
260  IF( m.LT.0 ) THEN
261  info = -1
262  ELSE IF( n.LT.0 ) THEN
263  info = -2
264  ELSE IF( nrhs.LT.0 ) THEN
265  info = -3
266  ELSE IF( lda.LT.max( 1, m ) ) THEN
267  info = -5
268  ELSE IF( ldb.LT.max( 1, maxmn ) ) THEN
269  info = -7
270  END IF
271 *
272  smlsiz = ilaenv( 9, 'DGELSD', ' ', 0, 0, 0, 0 )
273 *
274 * Compute workspace.
275 * (Note: Comments in the code beginning "Workspace:" describe the
276 * minimal amount of workspace needed at that point in the code,
277 * as well as the preferred amount for good performance.
278 * NB refers to the optimal block size for the immediately
279 * following subroutine, as returned by ILAENV.)
280 *
281  minwrk = 1
282  liwork = 1
283  minmn = max( 1, minmn )
284  nlvl = max( int( log( dble( minmn ) / dble( smlsiz+1 ) ) /
285  $ log( two ) ) + 1, 0 )
286 *
287  IF( info.EQ.0 ) THEN
288  maxwrk = 0
289  liwork = 3*minmn*nlvl + 11*minmn
290  mm = m
291  IF( m.GE.n .AND. m.GE.mnthr ) THEN
292 *
293 * Path 1a - overdetermined, with many more rows than columns.
294 *
295  mm = n
296  maxwrk = max( maxwrk, n+n*ilaenv( 1, 'DGEQRF', ' ', m, n,
297  $ -1, -1 ) )
298  maxwrk = max( maxwrk, n+nrhs*
299  $ ilaenv( 1, 'DORMQR', 'LT', m, nrhs, n, -1 ) )
300  END IF
301  IF( m.GE.n ) THEN
302 *
303 * Path 1 - overdetermined or exactly determined.
304 *
305  maxwrk = max( maxwrk, 3*n+( mm+n )*
306  $ ilaenv( 1, 'DGEBRD', ' ', mm, n, -1, -1 ) )
307  maxwrk = max( maxwrk, 3*n+nrhs*
308  $ ilaenv( 1, 'DORMBR', 'QLT', mm, nrhs, n, -1 ) )
309  maxwrk = max( maxwrk, 3*n+( n-1 )*
310  $ ilaenv( 1, 'DORMBR', 'PLN', n, nrhs, n, -1 ) )
311  wlalsd = 9*n+2*n*smlsiz+8*n*nlvl+n*nrhs+(smlsiz+1)**2
312  maxwrk = max( maxwrk, 3*n+wlalsd )
313  minwrk = max( 3*n+mm, 3*n+nrhs, 3*n+wlalsd )
314  END IF
315  IF( n.GT.m ) THEN
316  wlalsd = 9*m+2*m*smlsiz+8*m*nlvl+m*nrhs+(smlsiz+1)**2
317  IF( n.GE.mnthr ) THEN
318 *
319 * Path 2a - underdetermined, with many more columns
320 * than rows.
321 *
322  maxwrk = m + m*ilaenv( 1, 'DGELQF', ' ', m, n, -1, -1 )
323  maxwrk = max( maxwrk, m*m+4*m+2*m*
324  $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
325  maxwrk = max( maxwrk, m*m+4*m+nrhs*
326  $ ilaenv( 1, 'DORMBR', 'QLT', m, nrhs, m, -1 ) )
327  maxwrk = max( maxwrk, m*m+4*m+( m-1 )*
328  $ ilaenv( 1, 'DORMBR', 'PLN', m, nrhs, m, -1 ) )
329  IF( nrhs.GT.1 ) THEN
330  maxwrk = max( maxwrk, m*m+m+m*nrhs )
331  ELSE
332  maxwrk = max( maxwrk, m*m+2*m )
333  END IF
334  maxwrk = max( maxwrk, m+nrhs*
335  $ ilaenv( 1, 'DORMLQ', 'LT', n, nrhs, m, -1 ) )
336  maxwrk = max( maxwrk, m*m+4*m+wlalsd )
337 ! XXX: Ensure the Path 2a case below is triggered. The workspace
338 ! calculation should use queries for all routines eventually.
339  maxwrk = max( maxwrk,
340  $ 4*m+m*m+max( m, 2*m-4, nrhs, n-3*m ) )
341  ELSE
342 *
343 * Path 2 - remaining underdetermined cases.
344 *
345  maxwrk = 3*m + ( n+m )*ilaenv( 1, 'DGEBRD', ' ', m, n,
346  $ -1, -1 )
347  maxwrk = max( maxwrk, 3*m+nrhs*
348  $ ilaenv( 1, 'DORMBR', 'QLT', m, nrhs, n, -1 ) )
349  maxwrk = max( maxwrk, 3*m+m*
350  $ ilaenv( 1, 'DORMBR', 'PLN', n, nrhs, m, -1 ) )
351  maxwrk = max( maxwrk, 3*m+wlalsd )
352  END IF
353  minwrk = max( 3*m+nrhs, 3*m+m, 3*m+wlalsd )
354  END IF
355  minwrk = min( minwrk, maxwrk )
356  work( 1 ) = maxwrk
357  iwork( 1 ) = liwork
358 
359  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
360  info = -12
361  END IF
362  END IF
363 *
364  IF( info.NE.0 ) THEN
365  CALL xerbla( 'DGELSD', -info )
366  return
367  ELSE IF( lquery ) THEN
368  go to 10
369  END IF
370 *
371 * Quick return if possible.
372 *
373  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
374  rank = 0
375  return
376  END IF
377 *
378 * Get machine parameters.
379 *
380  eps = dlamch( 'P' )
381  sfmin = dlamch( 'S' )
382  smlnum = sfmin / eps
383  bignum = one / smlnum
384  CALL dlabad( smlnum, bignum )
385 *
386 * Scale A if max entry outside range [SMLNUM,BIGNUM].
387 *
388  anrm = dlange( 'M', m, n, a, lda, work )
389  iascl = 0
390  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
391 *
392 * Scale matrix norm up to SMLNUM.
393 *
394  CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
395  iascl = 1
396  ELSE IF( anrm.GT.bignum ) THEN
397 *
398 * Scale matrix norm down to BIGNUM.
399 *
400  CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
401  iascl = 2
402  ELSE IF( anrm.EQ.zero ) THEN
403 *
404 * Matrix all zero. Return zero solution.
405 *
406  CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
407  CALL dlaset( 'F', minmn, 1, zero, zero, s, 1 )
408  rank = 0
409  go to 10
410  END IF
411 *
412 * Scale B if max entry outside range [SMLNUM,BIGNUM].
413 *
414  bnrm = dlange( 'M', m, nrhs, b, ldb, work )
415  ibscl = 0
416  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
417 *
418 * Scale matrix norm up to SMLNUM.
419 *
420  CALL dlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
421  ibscl = 1
422  ELSE IF( bnrm.GT.bignum ) THEN
423 *
424 * Scale matrix norm down to BIGNUM.
425 *
426  CALL dlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
427  ibscl = 2
428  END IF
429 *
430 * If M < N make sure certain entries of B are zero.
431 *
432  IF( m.LT.n )
433  $ CALL dlaset( 'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
434 *
435 * Overdetermined case.
436 *
437  IF( m.GE.n ) THEN
438 *
439 * Path 1 - overdetermined or exactly determined.
440 *
441  mm = m
442  IF( m.GE.mnthr ) THEN
443 *
444 * Path 1a - overdetermined, with many more rows than columns.
445 *
446  mm = n
447  itau = 1
448  nwork = itau + n
449 *
450 * Compute A=Q*R.
451 * (Workspace: need 2*N, prefer N+N*NB)
452 *
453  CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
454  $ lwork-nwork+1, info )
455 *
456 * Multiply B by transpose(Q).
457 * (Workspace: need N+NRHS, prefer N+NRHS*NB)
458 *
459  CALL dormqr( 'L', 'T', m, nrhs, n, a, lda, work( itau ), b,
460  $ ldb, work( nwork ), lwork-nwork+1, info )
461 *
462 * Zero out below R.
463 *
464  IF( n.GT.1 ) THEN
465  CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
466  END IF
467  END IF
468 *
469  ie = 1
470  itauq = ie + n
471  itaup = itauq + n
472  nwork = itaup + n
473 *
474 * Bidiagonalize R in A.
475 * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
476 *
477  CALL dgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
478  $ work( itaup ), work( nwork ), lwork-nwork+1,
479  $ info )
480 *
481 * Multiply B by transpose of left bidiagonalizing vectors of R.
482 * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
483 *
484  CALL dormbr( 'Q', 'L', 'T', mm, nrhs, n, a, lda, work( itauq ),
485  $ b, ldb, work( nwork ), lwork-nwork+1, info )
486 *
487 * Solve the bidiagonal least squares problem.
488 *
489  CALL dlalsd( 'U', smlsiz, n, nrhs, s, work( ie ), b, ldb,
490  $ rcond, rank, work( nwork ), iwork, info )
491  IF( info.NE.0 ) THEN
492  go to 10
493  END IF
494 *
495 * Multiply B by right bidiagonalizing vectors of R.
496 *
497  CALL dormbr( 'P', 'L', 'N', n, nrhs, n, a, lda, work( itaup ),
498  $ b, ldb, work( nwork ), lwork-nwork+1, info )
499 *
500  ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
501  $ max( m, 2*m-4, nrhs, n-3*m, wlalsd ) ) THEN
502 *
503 * Path 2a - underdetermined, with many more columns than rows
504 * and sufficient workspace for an efficient algorithm.
505 *
506  ldwork = m
507  IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
508  $ m*lda+m+m*nrhs, 4*m+m*lda+wlalsd ) )ldwork = lda
509  itau = 1
510  nwork = m + 1
511 *
512 * Compute A=L*Q.
513 * (Workspace: need 2*M, prefer M+M*NB)
514 *
515  CALL dgelqf( m, n, a, lda, work( itau ), work( nwork ),
516  $ lwork-nwork+1, info )
517  il = nwork
518 *
519 * Copy L to WORK(IL), zeroing out above its diagonal.
520 *
521  CALL dlacpy( 'L', m, m, a, lda, work( il ), ldwork )
522  CALL dlaset( 'U', m-1, m-1, zero, zero, work( il+ldwork ),
523  $ ldwork )
524  ie = il + ldwork*m
525  itauq = ie + m
526  itaup = itauq + m
527  nwork = itaup + m
528 *
529 * Bidiagonalize L in WORK(IL).
530 * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
531 *
532  CALL dgebrd( m, m, work( il ), ldwork, s, work( ie ),
533  $ work( itauq ), work( itaup ), work( nwork ),
534  $ lwork-nwork+1, info )
535 *
536 * Multiply B by transpose of left bidiagonalizing vectors of L.
537 * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
538 *
539  CALL dormbr( 'Q', 'L', 'T', m, nrhs, m, work( il ), ldwork,
540  $ work( itauq ), b, ldb, work( nwork ),
541  $ lwork-nwork+1, info )
542 *
543 * Solve the bidiagonal least squares problem.
544 *
545  CALL dlalsd( 'U', smlsiz, m, nrhs, s, work( ie ), b, ldb,
546  $ rcond, rank, work( nwork ), iwork, info )
547  IF( info.NE.0 ) THEN
548  go to 10
549  END IF
550 *
551 * Multiply B by right bidiagonalizing vectors of L.
552 *
553  CALL dormbr( 'P', 'L', 'N', m, nrhs, m, work( il ), ldwork,
554  $ work( itaup ), b, ldb, work( nwork ),
555  $ lwork-nwork+1, info )
556 *
557 * Zero out below first M rows of B.
558 *
559  CALL dlaset( 'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
560  nwork = itau + m
561 *
562 * Multiply transpose(Q) by B.
563 * (Workspace: need M+NRHS, prefer M+NRHS*NB)
564 *
565  CALL dormlq( 'L', 'T', n, nrhs, m, a, lda, work( itau ), b,
566  $ ldb, work( nwork ), lwork-nwork+1, info )
567 *
568  ELSE
569 *
570 * Path 2 - remaining underdetermined cases.
571 *
572  ie = 1
573  itauq = ie + m
574  itaup = itauq + m
575  nwork = itaup + m
576 *
577 * Bidiagonalize A.
578 * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
579 *
580  CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
581  $ work( itaup ), work( nwork ), lwork-nwork+1,
582  $ info )
583 *
584 * Multiply B by transpose of left bidiagonalizing vectors.
585 * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
586 *
587  CALL dormbr( 'Q', 'L', 'T', m, nrhs, n, a, lda, work( itauq ),
588  $ b, ldb, work( nwork ), lwork-nwork+1, info )
589 *
590 * Solve the bidiagonal least squares problem.
591 *
592  CALL dlalsd( 'L', smlsiz, m, nrhs, s, work( ie ), b, ldb,
593  $ rcond, rank, work( nwork ), iwork, info )
594  IF( info.NE.0 ) THEN
595  go to 10
596  END IF
597 *
598 * Multiply B by right bidiagonalizing vectors of A.
599 *
600  CALL dormbr( 'P', 'L', 'N', n, nrhs, m, a, lda, work( itaup ),
601  $ b, ldb, work( nwork ), lwork-nwork+1, info )
602 *
603  END IF
604 *
605 * Undo scaling.
606 *
607  IF( iascl.EQ.1 ) THEN
608  CALL dlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
609  CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
610  $ info )
611  ELSE IF( iascl.EQ.2 ) THEN
612  CALL dlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
613  CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
614  $ info )
615  END IF
616  IF( ibscl.EQ.1 ) THEN
617  CALL dlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
618  ELSE IF( ibscl.EQ.2 ) THEN
619  CALL dlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
620  END IF
621 *
622  10 continue
623  work( 1 ) = maxwrk
624  iwork( 1 ) = liwork
625  return
626 *
627 * End of DGELSD
628 *
629  END