LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 transformations 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,out] 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 *> \ingroup realGEsolve
199 *
200 *> \par Contributors:
201 * ==================
202 *>
203 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
204 *> California at Berkeley, USA \n
205 *> Osni Marques, LBNL/NERSC, USA \n
206 *
207 * =====================================================================
208  SUBROUTINE sgelsd( M, N, NRHS, A, LDA, B, LDB, S, RCOND,
209  $ RANK, WORK, LWORK, IWORK, INFO )
210 *
211 * -- LAPACK driver routine --
212 * -- LAPACK is a software package provided by Univ. of Tennessee, --
213 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
214 *
215 * .. Scalar Arguments ..
216  INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
217  REAL RCOND
218 * ..
219 * .. Array Arguments ..
220  INTEGER IWORK( * )
221  REAL A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
222 * ..
223 *
224 * =====================================================================
225 *
226 * .. Parameters ..
227  REAL ZERO, ONE, TWO
228  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
229 * ..
230 * .. Local Scalars ..
231  LOGICAL LQUERY
232  INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
233  $ ldwork, liwork, maxmn, maxwrk, minmn, minwrk,
234  $ mm, mnthr, nlvl, nwork, smlsiz, wlalsd
235  REAL ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
236 * ..
237 * .. External Subroutines ..
238  EXTERNAL sgebrd, sgelqf, sgeqrf, slabad, slacpy, slalsd,
240 * ..
241 * .. External Functions ..
242  INTEGER ILAENV
243  REAL SLAMCH, SLANGE
244  EXTERNAL slamch, slange, ilaenv
245 * ..
246 * .. Intrinsic Functions ..
247  INTRINSIC int, log, max, min, real
248 * ..
249 * .. Executable Statements ..
250 *
251 * Test the input arguments.
252 *
253  info = 0
254  minmn = min( m, n )
255  maxmn = max( m, n )
256  lquery = ( lwork.EQ.-1 )
257  IF( m.LT.0 ) THEN
258  info = -1
259  ELSE IF( n.LT.0 ) THEN
260  info = -2
261  ELSE IF( nrhs.LT.0 ) THEN
262  info = -3
263  ELSE IF( lda.LT.max( 1, m ) ) THEN
264  info = -5
265  ELSE IF( ldb.LT.max( 1, maxmn ) ) THEN
266  info = -7
267  END IF
268 *
269 * Compute workspace.
270 * (Note: Comments in the code beginning "Workspace:" describe the
271 * minimal amount of workspace needed at that point in the code,
272 * as well as the preferred amount for good performance.
273 * NB refers to the optimal block size for the immediately
274 * following subroutine, as returned by ILAENV.)
275 *
276  IF( info.EQ.0 ) THEN
277  minwrk = 1
278  maxwrk = 1
279  liwork = 1
280  IF( minmn.GT.0 ) THEN
281  smlsiz = ilaenv( 9, 'SGELSD', ' ', 0, 0, 0, 0 )
282  mnthr = ilaenv( 6, 'SGELSD', ' ', m, n, nrhs, -1 )
283  nlvl = max( int( log( real( minmn ) / real( smlsiz + 1 ) ) /
284  $ log( two ) ) + 1, 0 )
285  liwork = 3*minmn*nlvl + 11*minmn
286  mm = m
287  IF( m.GE.n .AND. m.GE.mnthr ) THEN
288 *
289 * Path 1a - overdetermined, with many more rows than
290 * columns.
291 *
292  mm = n
293  maxwrk = max( maxwrk, n + n*ilaenv( 1, 'SGEQRF', ' ', m,
294  $ n, -1, -1 ) )
295  maxwrk = max( maxwrk, n + nrhs*ilaenv( 1, 'SORMQR', 'LT',
296  $ m, nrhs, n, -1 ) )
297  END IF
298  IF( m.GE.n ) THEN
299 *
300 * Path 1 - overdetermined or exactly determined.
301 *
302  maxwrk = max( maxwrk, 3*n + ( mm + n )*ilaenv( 1,
303  $ 'SGEBRD', ' ', mm, n, -1, -1 ) )
304  maxwrk = max( maxwrk, 3*n + nrhs*ilaenv( 1, 'SORMBR',
305  $ 'QLT', mm, nrhs, n, -1 ) )
306  maxwrk = max( maxwrk, 3*n + ( n - 1 )*ilaenv( 1,
307  $ 'SORMBR', 'PLN', n, nrhs, n, -1 ) )
308  wlalsd = 9*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs +
309  $ ( smlsiz + 1 )**2
310  maxwrk = max( maxwrk, 3*n + wlalsd )
311  minwrk = max( 3*n + mm, 3*n + nrhs, 3*n + wlalsd )
312  END IF
313  IF( n.GT.m ) THEN
314  wlalsd = 9*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs +
315  $ ( smlsiz + 1 )**2
316  IF( n.GE.mnthr ) THEN
317 *
318 * Path 2a - underdetermined, with many more columns
319 * than rows.
320 *
321  maxwrk = m + m*ilaenv( 1, 'SGELQF', ' ', m, n, -1,
322  $ -1 )
323  maxwrk = max( maxwrk, m*m + 4*m + 2*m*ilaenv( 1,
324  $ 'SGEBRD', ' ', m, m, -1, -1 ) )
325  maxwrk = max( maxwrk, m*m + 4*m + nrhs*ilaenv( 1,
326  $ 'SORMBR', 'QLT', m, nrhs, m, -1 ) )
327  maxwrk = max( maxwrk, m*m + 4*m + ( m - 1 )*ilaenv( 1,
328  $ 'SORMBR', '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*ilaenv( 1, 'SORMLQ',
335  $ '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, 'SGEBRD', ' ', m,
346  $ n, -1, -1 )
347  maxwrk = max( maxwrk, 3*m + nrhs*ilaenv( 1, 'SORMBR',
348  $ 'QLT', m, nrhs, n, -1 ) )
349  maxwrk = max( maxwrk, 3*m + m*ilaenv( 1, 'SORMBR',
350  $ '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  END IF
356  minwrk = min( minwrk, maxwrk )
357  work( 1 ) = maxwrk
358  iwork( 1 ) = liwork
359 *
360  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
361  info = -12
362  END IF
363  END IF
364 *
365  IF( info.NE.0 ) THEN
366  CALL xerbla( 'SGELSD', -info )
367  RETURN
368  ELSE IF( lquery ) THEN
369  RETURN
370  END IF
371 *
372 * Quick return if possible.
373 *
374  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
375  rank = 0
376  RETURN
377  END IF
378 *
379 * Get machine parameters.
380 *
381  eps = slamch( 'P' )
382  sfmin = slamch( 'S' )
383  smlnum = sfmin / eps
384  bignum = one / smlnum
385  CALL slabad( smlnum, bignum )
386 *
387 * Scale A if max entry outside range [SMLNUM,BIGNUM].
388 *
389  anrm = slange( 'M', m, n, a, lda, work )
390  iascl = 0
391  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
392 *
393 * Scale matrix norm up to SMLNUM.
394 *
395  CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
396  iascl = 1
397  ELSE IF( anrm.GT.bignum ) THEN
398 *
399 * Scale matrix norm down to BIGNUM.
400 *
401  CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
402  iascl = 2
403  ELSE IF( anrm.EQ.zero ) THEN
404 *
405 * Matrix all zero. Return zero solution.
406 *
407  CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
408  CALL slaset( 'F', minmn, 1, zero, zero, s, 1 )
409  rank = 0
410  GO TO 10
411  END IF
412 *
413 * Scale B if max entry outside range [SMLNUM,BIGNUM].
414 *
415  bnrm = slange( 'M', m, nrhs, b, ldb, work )
416  ibscl = 0
417  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
418 *
419 * Scale matrix norm up to SMLNUM.
420 *
421  CALL slascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
422  ibscl = 1
423  ELSE IF( bnrm.GT.bignum ) THEN
424 *
425 * Scale matrix norm down to BIGNUM.
426 *
427  CALL slascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
428  ibscl = 2
429  END IF
430 *
431 * If M < N make sure certain entries of B are zero.
432 *
433  IF( m.LT.n )
434  $ CALL slaset( 'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
435 *
436 * Overdetermined case.
437 *
438  IF( m.GE.n ) THEN
439 *
440 * Path 1 - overdetermined or exactly determined.
441 *
442  mm = m
443  IF( m.GE.mnthr ) THEN
444 *
445 * Path 1a - overdetermined, with many more rows than columns.
446 *
447  mm = n
448  itau = 1
449  nwork = itau + n
450 *
451 * Compute A=Q*R.
452 * (Workspace: need 2*N, prefer N+N*NB)
453 *
454  CALL sgeqrf( m, n, a, lda, work( itau ), work( nwork ),
455  $ lwork-nwork+1, info )
456 *
457 * Multiply B by transpose(Q).
458 * (Workspace: need N+NRHS, prefer N+NRHS*NB)
459 *
460  CALL sormqr( 'L', 'T', m, nrhs, n, a, lda, work( itau ), b,
461  $ ldb, work( nwork ), lwork-nwork+1, info )
462 *
463 * Zero out below R.
464 *
465  IF( n.GT.1 ) THEN
466  CALL slaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
467  END IF
468  END IF
469 *
470  ie = 1
471  itauq = ie + n
472  itaup = itauq + n
473  nwork = itaup + n
474 *
475 * Bidiagonalize R in A.
476 * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
477 *
478  CALL sgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
479  $ work( itaup ), work( nwork ), lwork-nwork+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 sormbr( 'Q', 'L', 'T', mm, nrhs, n, a, lda, work( itauq ),
486  $ b, ldb, work( nwork ), lwork-nwork+1, info )
487 *
488 * Solve the bidiagonal least squares problem.
489 *
490  CALL slalsd( 'U', smlsiz, n, nrhs, s, work( ie ), b, ldb,
491  $ rcond, rank, work( nwork ), iwork, info )
492  IF( info.NE.0 ) THEN
493  GO TO 10
494  END IF
495 *
496 * Multiply B by right bidiagonalizing vectors of R.
497 *
498  CALL sormbr( 'P', 'L', 'N', n, nrhs, n, a, lda, work( itaup ),
499  $ b, ldb, work( nwork ), lwork-nwork+1, info )
500 *
501  ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
502  $ max( m, 2*m-4, nrhs, n-3*m, wlalsd ) ) THEN
503 *
504 * Path 2a - underdetermined, with many more columns than rows
505 * and sufficient workspace for an efficient algorithm.
506 *
507  ldwork = m
508  IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
509  $ m*lda+m+m*nrhs, 4*m+m*lda+wlalsd ) )ldwork = lda
510  itau = 1
511  nwork = m + 1
512 *
513 * Compute A=L*Q.
514 * (Workspace: need 2*M, prefer M+M*NB)
515 *
516  CALL sgelqf( m, n, a, lda, work( itau ), work( nwork ),
517  $ lwork-nwork+1, info )
518  il = nwork
519 *
520 * Copy L to WORK(IL), zeroing out above its diagonal.
521 *
522  CALL slacpy( 'L', m, m, a, lda, work( il ), ldwork )
523  CALL slaset( 'U', m-1, m-1, zero, zero, work( il+ldwork ),
524  $ ldwork )
525  ie = il + ldwork*m
526  itauq = ie + m
527  itaup = itauq + m
528  nwork = itaup + m
529 *
530 * Bidiagonalize L in WORK(IL).
531 * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
532 *
533  CALL sgebrd( m, m, work( il ), ldwork, s, work( ie ),
534  $ work( itauq ), work( itaup ), work( nwork ),
535  $ lwork-nwork+1, info )
536 *
537 * Multiply B by transpose of left bidiagonalizing vectors of L.
538 * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
539 *
540  CALL sormbr( 'Q', 'L', 'T', m, nrhs, m, work( il ), ldwork,
541  $ work( itauq ), b, ldb, work( nwork ),
542  $ lwork-nwork+1, info )
543 *
544 * Solve the bidiagonal least squares problem.
545 *
546  CALL slalsd( 'U', smlsiz, m, nrhs, s, work( ie ), b, ldb,
547  $ rcond, rank, work( nwork ), iwork, info )
548  IF( info.NE.0 ) THEN
549  GO TO 10
550  END IF
551 *
552 * Multiply B by right bidiagonalizing vectors of L.
553 *
554  CALL sormbr( 'P', 'L', 'N', m, nrhs, m, work( il ), ldwork,
555  $ work( itaup ), b, ldb, work( nwork ),
556  $ lwork-nwork+1, info )
557 *
558 * Zero out below first M rows of B.
559 *
560  CALL slaset( 'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
561  nwork = itau + m
562 *
563 * Multiply transpose(Q) by B.
564 * (Workspace: need M+NRHS, prefer M+NRHS*NB)
565 *
566  CALL sormlq( 'L', 'T', n, nrhs, m, a, lda, work( itau ), b,
567  $ ldb, work( nwork ), lwork-nwork+1, info )
568 *
569  ELSE
570 *
571 * Path 2 - remaining underdetermined cases.
572 *
573  ie = 1
574  itauq = ie + m
575  itaup = itauq + m
576  nwork = itaup + m
577 *
578 * Bidiagonalize A.
579 * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
580 *
581  CALL sgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
582  $ work( itaup ), work( nwork ), lwork-nwork+1,
583  $ info )
584 *
585 * Multiply B by transpose of left bidiagonalizing vectors.
586 * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
587 *
588  CALL sormbr( 'Q', 'L', 'T', m, nrhs, n, a, lda, work( itauq ),
589  $ b, ldb, work( nwork ), lwork-nwork+1, info )
590 *
591 * Solve the bidiagonal least squares problem.
592 *
593  CALL slalsd( 'L', smlsiz, m, nrhs, s, work( ie ), b, ldb,
594  $ rcond, rank, work( nwork ), iwork, info )
595  IF( info.NE.0 ) THEN
596  GO TO 10
597  END IF
598 *
599 * Multiply B by right bidiagonalizing vectors of A.
600 *
601  CALL sormbr( 'P', 'L', 'N', n, nrhs, m, a, lda, work( itaup ),
602  $ b, ldb, work( nwork ), lwork-nwork+1, info )
603 *
604  END IF
605 *
606 * Undo scaling.
607 *
608  IF( iascl.EQ.1 ) THEN
609  CALL slascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
610  CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
611  $ info )
612  ELSE IF( iascl.EQ.2 ) THEN
613  CALL slascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
614  CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
615  $ info )
616  END IF
617  IF( ibscl.EQ.1 ) THEN
618  CALL slascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
619  ELSE IF( ibscl.EQ.2 ) THEN
620  CALL slascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
621  END IF
622 *
623  10 CONTINUE
624  work( 1 ) = maxwrk
625  iwork( 1 ) = liwork
626  RETURN
627 *
628 * End of SGELSD
629 *
630  END
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:143
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
Definition: sgeqrf.f:145
subroutine sgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
SGEBRD
Definition: sgebrd.f:205
subroutine sgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGELQF
Definition: sgelqf.f:143
subroutine sgelsd(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, IWORK, INFO)
SGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices
Definition: sgelsd.f:210
subroutine slalsd(UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, RANK, WORK, IWORK, INFO)
SLALSD uses the singular value decomposition of A to solve the least squares problem.
Definition: slalsd.f:179
subroutine sormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMBR
Definition: sormbr.f:196
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
Definition: sormqr.f:168
subroutine sormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMLQ
Definition: sormlq.f:168