LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
sgelsd.f
Go to the documentation of this file.
1 *> \brief <b> SGELSD 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 SGELSD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgelsd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgelsd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgelsd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND,
22 * RANK, WORK, LWORK, IWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
26 * REAL RCOND
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IWORK( * )
30 * REAL A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> SGELSD 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 REAL 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 REAL 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 REAL 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 REAL
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 REAL 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 array WORK and the
167 *> minimum size of the array IWORK, and returns these values as
168 *> the first entries of the WORK and IWORK arrays, and no error
169 *> message related to LWORK is issued by XERBLA.
170 *> \endverbatim
171 *>
172 *> \param[out] IWORK
173 *> \verbatim
174 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
175 *> LIWORK >= max(1, 3*MINMN*NLVL + 11*MINMN),
176 *> where MINMN = MIN( M,N ).
177 *> On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
178 *> \endverbatim
179 *>
180 *> \param[out] INFO
181 *> \verbatim
182 *> INFO is INTEGER
183 *> = 0: successful exit
184 *> < 0: if INFO = -i, the i-th argument had an illegal value.
185 *> > 0: the algorithm for computing the SVD failed to converge;
186 *> if INFO = i, i off-diagonal elements of an intermediate
187 *> bidiagonal form did not converge to zero.
188 *> \endverbatim
189 *
190 * Authors:
191 * ========
192 *
193 *> \author Univ. of Tennessee
194 *> \author Univ. of California Berkeley
195 *> \author Univ. of Colorado Denver
196 *> \author NAG Ltd.
197 *
198 *> \date November 2011
199 *
200 *> \ingroup realGEsolve
201 *
202 *> \par Contributors:
203 * ==================
204 *>
205 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
206 *> California at Berkeley, USA \n
207 *> Osni Marques, LBNL/NERSC, USA \n
208 *
209 * =====================================================================
210  SUBROUTINE sgelsd( M, N, NRHS, A, LDA, B, LDB, S, RCOND,
211  $ rank, work, lwork, iwork, info )
212 *
213 * -- LAPACK driver routine (version 3.4.0) --
214 * -- LAPACK is a software package provided by Univ. of Tennessee, --
215 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
216 * November 2011
217 *
218 * .. Scalar Arguments ..
219  INTEGER info, lda, ldb, lwork, m, n, nrhs, rank
220  REAL rcond
221 * ..
222 * .. Array Arguments ..
223  INTEGER iwork( * )
224  REAL a( lda, * ), b( ldb, * ), s( * ), work( * )
225 * ..
226 *
227 * =====================================================================
228 *
229 * .. Parameters ..
230  REAL zero, one, two
231  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
232 * ..
233 * .. Local Scalars ..
234  LOGICAL lquery
235  INTEGER iascl, ibscl, ie, il, itau, itaup, itauq,
236  $ ldwork, liwork, maxmn, maxwrk, minmn, minwrk,
237  $ mm, mnthr, nlvl, nwork, smlsiz, wlalsd
238  REAL anrm, bignum, bnrm, eps, sfmin, smlnum
239 * ..
240 * .. External Subroutines ..
241  EXTERNAL sgebrd, sgelqf, sgeqrf, slabad, slacpy, slalsd,
243 * ..
244 * .. External Functions ..
245  INTEGER ilaenv
246  REAL slamch, slange
247  EXTERNAL slamch, slange, ilaenv
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC int, log, max, min, real
251 * ..
252 * .. Executable Statements ..
253 *
254 * Test the input arguments.
255 *
256  info = 0
257  minmn = min( m, n )
258  maxmn = max( m, n )
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 * Compute workspace.
273 * (Note: Comments in the code beginning "Workspace:" describe the
274 * minimal amount of workspace needed at that point in the code,
275 * as well as the preferred amount for good performance.
276 * NB refers to the optimal block size for the immediately
277 * following subroutine, as returned by ILAENV.)
278 *
279  IF( info.EQ.0 ) THEN
280  minwrk = 1
281  maxwrk = 1
282  liwork = 1
283  IF( minmn.GT.0 ) THEN
284  smlsiz = ilaenv( 9, 'SGELSD', ' ', 0, 0, 0, 0 )
285  mnthr = ilaenv( 6, 'SGELSD', ' ', m, n, nrhs, -1 )
286  nlvl = max( int( log( REAL( MINMN ) / REAL( SMLSIZ + 1 ) ) /
287  $ log( two ) ) + 1, 0 )
288  liwork = 3*minmn*nlvl + 11*minmn
289  mm = m
290  IF( m.GE.n .AND. m.GE.mnthr ) THEN
291 *
292 * Path 1a - overdetermined, with many more rows than
293 * columns.
294 *
295  mm = n
296  maxwrk = max( maxwrk, n + n*ilaenv( 1, 'SGEQRF', ' ', m,
297  $ n, -1, -1 ) )
298  maxwrk = max( maxwrk, n + nrhs*ilaenv( 1, 'SORMQR', 'LT',
299  $ 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 )*ilaenv( 1,
306  $ 'SGEBRD', ' ', mm, n, -1, -1 ) )
307  maxwrk = max( maxwrk, 3*n + nrhs*ilaenv( 1, 'SORMBR',
308  $ 'QLT', mm, nrhs, n, -1 ) )
309  maxwrk = max( maxwrk, 3*n + ( n - 1 )*ilaenv( 1,
310  $ 'SORMBR', 'PLN', n, nrhs, n, -1 ) )
311  wlalsd = 9*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs +
312  $ ( smlsiz + 1 )**2
313  maxwrk = max( maxwrk, 3*n + wlalsd )
314  minwrk = max( 3*n + mm, 3*n + nrhs, 3*n + wlalsd )
315  END IF
316  IF( n.GT.m ) THEN
317  wlalsd = 9*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs +
318  $ ( smlsiz + 1 )**2
319  IF( n.GE.mnthr ) THEN
320 *
321 * Path 2a - underdetermined, with many more columns
322 * than rows.
323 *
324  maxwrk = m + m*ilaenv( 1, 'SGELQF', ' ', m, n, -1,
325  $ -1 )
326  maxwrk = max( maxwrk, m*m + 4*m + 2*m*ilaenv( 1,
327  $ 'SGEBRD', ' ', m, m, -1, -1 ) )
328  maxwrk = max( maxwrk, m*m + 4*m + nrhs*ilaenv( 1,
329  $ 'SORMBR', 'QLT', m, nrhs, m, -1 ) )
330  maxwrk = max( maxwrk, m*m + 4*m + ( m - 1 )*ilaenv( 1,
331  $ 'SORMBR', 'PLN', m, nrhs, m, -1 ) )
332  IF( nrhs.GT.1 ) THEN
333  maxwrk = max( maxwrk, m*m + m + m*nrhs )
334  ELSE
335  maxwrk = max( maxwrk, m*m + 2*m )
336  END IF
337  maxwrk = max( maxwrk, m + nrhs*ilaenv( 1, 'SORMLQ',
338  $ 'LT', n, nrhs, m, -1 ) )
339  maxwrk = max( maxwrk, m*m + 4*m + wlalsd )
340 ! XXX: Ensure the Path 2a case below is triggered. The workspace
341 ! calculation should use queries for all routines eventually.
342  maxwrk = max( maxwrk,
343  $ 4*m+m*m+max( m, 2*m-4, nrhs, n-3*m ) )
344  ELSE
345 *
346 * Path 2 - remaining underdetermined cases.
347 *
348  maxwrk = 3*m + ( n + m )*ilaenv( 1, 'SGEBRD', ' ', m,
349  $ n, -1, -1 )
350  maxwrk = max( maxwrk, 3*m + nrhs*ilaenv( 1, 'SORMBR',
351  $ 'QLT', m, nrhs, n, -1 ) )
352  maxwrk = max( maxwrk, 3*m + m*ilaenv( 1, 'SORMBR',
353  $ 'PLN', n, nrhs, m, -1 ) )
354  maxwrk = max( maxwrk, 3*m + wlalsd )
355  END IF
356  minwrk = max( 3*m + nrhs, 3*m + m, 3*m + wlalsd )
357  END IF
358  END IF
359  minwrk = min( minwrk, maxwrk )
360  work( 1 ) = maxwrk
361  iwork( 1 ) = liwork
362 *
363  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
364  info = -12
365  END IF
366  END IF
367 *
368  IF( info.NE.0 ) THEN
369  CALL xerbla( 'SGELSD', -info )
370  RETURN
371  ELSE IF( lquery ) THEN
372  RETURN
373  END IF
374 *
375 * Quick return if possible.
376 *
377  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
378  rank = 0
379  RETURN
380  END IF
381 *
382 * Get machine parameters.
383 *
384  eps = slamch( 'P' )
385  sfmin = slamch( 'S' )
386  smlnum = sfmin / eps
387  bignum = one / smlnum
388  CALL slabad( smlnum, bignum )
389 *
390 * Scale A if max entry outside range [SMLNUM,BIGNUM].
391 *
392  anrm = slange( 'M', m, n, a, lda, work )
393  iascl = 0
394  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
395 *
396 * Scale matrix norm up to SMLNUM.
397 *
398  CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
399  iascl = 1
400  ELSE IF( anrm.GT.bignum ) THEN
401 *
402 * Scale matrix norm down to BIGNUM.
403 *
404  CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
405  iascl = 2
406  ELSE IF( anrm.EQ.zero ) THEN
407 *
408 * Matrix all zero. Return zero solution.
409 *
410  CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
411  CALL slaset( 'F', minmn, 1, zero, zero, s, 1 )
412  rank = 0
413  go to 10
414  END IF
415 *
416 * Scale B if max entry outside range [SMLNUM,BIGNUM].
417 *
418  bnrm = slange( 'M', m, nrhs, b, ldb, work )
419  ibscl = 0
420  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
421 *
422 * Scale matrix norm up to SMLNUM.
423 *
424  CALL slascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
425  ibscl = 1
426  ELSE IF( bnrm.GT.bignum ) THEN
427 *
428 * Scale matrix norm down to BIGNUM.
429 *
430  CALL slascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
431  ibscl = 2
432  END IF
433 *
434 * If M < N make sure certain entries of B are zero.
435 *
436  IF( m.LT.n )
437  $ CALL slaset( 'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
438 *
439 * Overdetermined case.
440 *
441  IF( m.GE.n ) THEN
442 *
443 * Path 1 - overdetermined or exactly determined.
444 *
445  mm = m
446  IF( m.GE.mnthr ) THEN
447 *
448 * Path 1a - overdetermined, with many more rows than columns.
449 *
450  mm = n
451  itau = 1
452  nwork = itau + n
453 *
454 * Compute A=Q*R.
455 * (Workspace: need 2*N, prefer N+N*NB)
456 *
457  CALL sgeqrf( m, n, a, lda, work( itau ), work( nwork ),
458  $ lwork-nwork+1, info )
459 *
460 * Multiply B by transpose(Q).
461 * (Workspace: need N+NRHS, prefer N+NRHS*NB)
462 *
463  CALL sormqr( 'L', 'T', m, nrhs, n, a, lda, work( itau ), b,
464  $ ldb, work( nwork ), lwork-nwork+1, info )
465 *
466 * Zero out below R.
467 *
468  IF( n.GT.1 ) THEN
469  CALL slaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
470  END IF
471  END IF
472 *
473  ie = 1
474  itauq = ie + n
475  itaup = itauq + n
476  nwork = itaup + n
477 *
478 * Bidiagonalize R in A.
479 * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
480 *
481  CALL sgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
482  $ work( itaup ), work( nwork ), lwork-nwork+1,
483  $ info )
484 *
485 * Multiply B by transpose of left bidiagonalizing vectors of R.
486 * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
487 *
488  CALL sormbr( 'Q', 'L', 'T', mm, nrhs, n, a, lda, work( itauq ),
489  $ b, ldb, work( nwork ), lwork-nwork+1, info )
490 *
491 * Solve the bidiagonal least squares problem.
492 *
493  CALL slalsd( 'U', smlsiz, n, nrhs, s, work( ie ), b, ldb,
494  $ rcond, rank, work( nwork ), iwork, info )
495  IF( info.NE.0 ) THEN
496  go to 10
497  END IF
498 *
499 * Multiply B by right bidiagonalizing vectors of R.
500 *
501  CALL sormbr( 'P', 'L', 'N', n, nrhs, n, a, lda, work( itaup ),
502  $ b, ldb, work( nwork ), lwork-nwork+1, info )
503 *
504  ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
505  $ max( m, 2*m-4, nrhs, n-3*m, wlalsd ) ) THEN
506 *
507 * Path 2a - underdetermined, with many more columns than rows
508 * and sufficient workspace for an efficient algorithm.
509 *
510  ldwork = m
511  IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
512  $ m*lda+m+m*nrhs, 4*m+m*lda+wlalsd ) )ldwork = lda
513  itau = 1
514  nwork = m + 1
515 *
516 * Compute A=L*Q.
517 * (Workspace: need 2*M, prefer M+M*NB)
518 *
519  CALL sgelqf( m, n, a, lda, work( itau ), work( nwork ),
520  $ lwork-nwork+1, info )
521  il = nwork
522 *
523 * Copy L to WORK(IL), zeroing out above its diagonal.
524 *
525  CALL slacpy( 'L', m, m, a, lda, work( il ), ldwork )
526  CALL slaset( 'U', m-1, m-1, zero, zero, work( il+ldwork ),
527  $ ldwork )
528  ie = il + ldwork*m
529  itauq = ie + m
530  itaup = itauq + m
531  nwork = itaup + m
532 *
533 * Bidiagonalize L in WORK(IL).
534 * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
535 *
536  CALL sgebrd( m, m, work( il ), ldwork, s, work( ie ),
537  $ work( itauq ), work( itaup ), work( nwork ),
538  $ lwork-nwork+1, info )
539 *
540 * Multiply B by transpose of left bidiagonalizing vectors of L.
541 * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
542 *
543  CALL sormbr( 'Q', 'L', 'T', m, nrhs, m, work( il ), ldwork,
544  $ work( itauq ), b, ldb, work( nwork ),
545  $ lwork-nwork+1, info )
546 *
547 * Solve the bidiagonal least squares problem.
548 *
549  CALL slalsd( 'U', smlsiz, m, nrhs, s, work( ie ), b, ldb,
550  $ rcond, rank, work( nwork ), iwork, info )
551  IF( info.NE.0 ) THEN
552  go to 10
553  END IF
554 *
555 * Multiply B by right bidiagonalizing vectors of L.
556 *
557  CALL sormbr( 'P', 'L', 'N', m, nrhs, m, work( il ), ldwork,
558  $ work( itaup ), b, ldb, work( nwork ),
559  $ lwork-nwork+1, info )
560 *
561 * Zero out below first M rows of B.
562 *
563  CALL slaset( 'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
564  nwork = itau + m
565 *
566 * Multiply transpose(Q) by B.
567 * (Workspace: need M+NRHS, prefer M+NRHS*NB)
568 *
569  CALL sormlq( 'L', 'T', n, nrhs, m, a, lda, work( itau ), b,
570  $ ldb, work( nwork ), lwork-nwork+1, info )
571 *
572  ELSE
573 *
574 * Path 2 - remaining underdetermined cases.
575 *
576  ie = 1
577  itauq = ie + m
578  itaup = itauq + m
579  nwork = itaup + m
580 *
581 * Bidiagonalize A.
582 * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
583 *
584  CALL sgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
585  $ work( itaup ), work( nwork ), lwork-nwork+1,
586  $ info )
587 *
588 * Multiply B by transpose of left bidiagonalizing vectors.
589 * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
590 *
591  CALL sormbr( 'Q', 'L', 'T', m, nrhs, n, a, lda, work( itauq ),
592  $ b, ldb, work( nwork ), lwork-nwork+1, info )
593 *
594 * Solve the bidiagonal least squares problem.
595 *
596  CALL slalsd( 'L', smlsiz, m, nrhs, s, work( ie ), b, ldb,
597  $ rcond, rank, work( nwork ), iwork, info )
598  IF( info.NE.0 ) THEN
599  go to 10
600  END IF
601 *
602 * Multiply B by right bidiagonalizing vectors of A.
603 *
604  CALL sormbr( 'P', 'L', 'N', n, nrhs, m, a, lda, work( itaup ),
605  $ b, ldb, work( nwork ), lwork-nwork+1, info )
606 *
607  END IF
608 *
609 * Undo scaling.
610 *
611  IF( iascl.EQ.1 ) THEN
612  CALL slascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
613  CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
614  $ info )
615  ELSE IF( iascl.EQ.2 ) THEN
616  CALL slascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
617  CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
618  $ info )
619  END IF
620  IF( ibscl.EQ.1 ) THEN
621  CALL slascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
622  ELSE IF( ibscl.EQ.2 ) THEN
623  CALL slascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
624  END IF
625 *
626  10 CONTINUE
627  work( 1 ) = maxwrk
628  iwork( 1 ) = liwork
629  RETURN
630 *
631 * End of SGELSD
632 *
633  END