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