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