LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dgelss()

subroutine dgelss ( integer  M,
integer  N,
integer  NRHS,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( * )  S,
double precision  RCOND,
integer  RANK,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DGELSS solves overdetermined or underdetermined systems for GE matrices

Download DGELSS + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DGELSS computes the minimum norm solution to a real linear least
 squares problem:

 Minimize 2-norm(| b - A*x |).

 using the singular value decomposition (SVD) of A. A is an M-by-N
 matrix which may be rank-deficient.

 Several right hand side vectors b and solution vectors x can be
 handled in a single call; they are stored as the columns of the
 M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix
 X.

 The effective rank of A is determined by treating as zero those
 singular values which are less than RCOND times the largest singular
 value.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A. M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A. N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrices B and X. NRHS >= 0.
[in,out]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, the first min(m,n) rows of A are overwritten with
          its right singular vectors, stored rowwise.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[in,out]B
          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
          On entry, the M-by-NRHS right hand side matrix B.
          On exit, B is overwritten by the N-by-NRHS solution
          matrix X.  If m >= n and RANK = n, the residual
          sum-of-squares for the solution in the i-th column is given
          by the sum of squares of elements n+1:m in that column.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1,max(M,N)).
[out]S
          S is DOUBLE PRECISION array, dimension (min(M,N))
          The singular values of A in decreasing order.
          The condition number of A in the 2-norm = S(1)/S(min(m,n)).
[in]RCOND
          RCOND is DOUBLE PRECISION
          RCOND is used to determine the effective rank of A.
          Singular values S(i) <= RCOND*S(1) are treated as zero.
          If RCOND < 0, machine precision is used instead.
[out]RANK
          RANK is INTEGER
          The effective rank of A, i.e., the number of singular values
          which are greater than RCOND*S(1).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= 1, and also:
          LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS )
          For good performance, LWORK should generally be larger.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  the algorithm for computing the SVD failed to converge;
                if INFO = i, i off-diagonal elements of an intermediate
                bidiagonal form did not converge to zero.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 170 of file dgelss.f.

172*
173* -- LAPACK driver routine --
174* -- LAPACK is a software package provided by Univ. of Tennessee, --
175* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
176*
177* .. Scalar Arguments ..
178 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
179 DOUBLE PRECISION RCOND
180* ..
181* .. Array Arguments ..
182 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
183* ..
184*
185* =====================================================================
186*
187* .. Parameters ..
188 DOUBLE PRECISION ZERO, ONE
189 parameter( zero = 0.0d+0, one = 1.0d+0 )
190* ..
191* .. Local Scalars ..
192 LOGICAL LQUERY
193 INTEGER BDSPAC, BL, CHUNK, I, IASCL, IBSCL, IE, IL,
194 $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN,
195 $ MAXWRK, MINMN, MINWRK, MM, MNTHR
196 INTEGER LWORK_DGEQRF, LWORK_DORMQR, LWORK_DGEBRD,
197 $ LWORK_DORMBR, LWORK_DORGBR, LWORK_DORMLQ,
198 $ LWORK_DGELQF
199 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
200* ..
201* .. Local Arrays ..
202 DOUBLE PRECISION DUM( 1 )
203* ..
204* .. External Subroutines ..
205 EXTERNAL dbdsqr, dcopy, dgebrd, dgelqf, dgemm, dgemv,
206 $ dgeqrf, dlabad, dlacpy, dlascl, dlaset, dorgbr,
207 $ dormbr, dormlq, dormqr, drscl, xerbla
208* ..
209* .. External Functions ..
210 INTEGER ILAENV
211 DOUBLE PRECISION DLAMCH, DLANGE
212 EXTERNAL ilaenv, dlamch, dlange
213* ..
214* .. Intrinsic Functions ..
215 INTRINSIC max, min
216* ..
217* .. Executable Statements ..
218*
219* Test the input arguments
220*
221 info = 0
222 minmn = min( m, n )
223 maxmn = max( m, n )
224 lquery = ( lwork.EQ.-1 )
225 IF( m.LT.0 ) THEN
226 info = -1
227 ELSE IF( n.LT.0 ) THEN
228 info = -2
229 ELSE IF( nrhs.LT.0 ) THEN
230 info = -3
231 ELSE IF( lda.LT.max( 1, m ) ) THEN
232 info = -5
233 ELSE IF( ldb.LT.max( 1, maxmn ) ) THEN
234 info = -7
235 END IF
236*
237* Compute workspace
238* (Note: Comments in the code beginning "Workspace:" describe the
239* minimal amount of workspace needed at that point in the code,
240* as well as the preferred amount for good performance.
241* NB refers to the optimal block size for the immediately
242* following subroutine, as returned by ILAENV.)
243*
244 IF( info.EQ.0 ) THEN
245 minwrk = 1
246 maxwrk = 1
247 IF( minmn.GT.0 ) THEN
248 mm = m
249 mnthr = ilaenv( 6, 'DGELSS', ' ', m, n, nrhs, -1 )
250 IF( m.GE.n .AND. m.GE.mnthr ) THEN
251*
252* Path 1a - overdetermined, with many more rows than
253* columns
254*
255* Compute space needed for DGEQRF
256 CALL dgeqrf( m, n, a, lda, dum(1), dum(1), -1, info )
257 lwork_dgeqrf = int( dum(1) )
258* Compute space needed for DORMQR
259 CALL dormqr( 'L', 'T', m, nrhs, n, a, lda, dum(1), b,
260 $ ldb, dum(1), -1, info )
261 lwork_dormqr = int( dum(1) )
262 mm = n
263 maxwrk = max( maxwrk, n + lwork_dgeqrf )
264 maxwrk = max( maxwrk, n + lwork_dormqr )
265 END IF
266 IF( m.GE.n ) THEN
267*
268* Path 1 - overdetermined or exactly determined
269*
270* Compute workspace needed for DBDSQR
271*
272 bdspac = max( 1, 5*n )
273* Compute space needed for DGEBRD
274 CALL dgebrd( mm, n, a, lda, s, dum(1), dum(1),
275 $ dum(1), dum(1), -1, info )
276 lwork_dgebrd = int( dum(1) )
277* Compute space needed for DORMBR
278 CALL dormbr( 'Q', 'L', 'T', mm, nrhs, n, a, lda, dum(1),
279 $ b, ldb, dum(1), -1, info )
280 lwork_dormbr = int( dum(1) )
281* Compute space needed for DORGBR
282 CALL dorgbr( 'P', n, n, n, a, lda, dum(1),
283 $ dum(1), -1, info )
284 lwork_dorgbr = int( dum(1) )
285* Compute total workspace needed
286 maxwrk = max( maxwrk, 3*n + lwork_dgebrd )
287 maxwrk = max( maxwrk, 3*n + lwork_dormbr )
288 maxwrk = max( maxwrk, 3*n + lwork_dorgbr )
289 maxwrk = max( maxwrk, bdspac )
290 maxwrk = max( maxwrk, n*nrhs )
291 minwrk = max( 3*n + mm, 3*n + nrhs, bdspac )
292 maxwrk = max( minwrk, maxwrk )
293 END IF
294 IF( n.GT.m ) THEN
295*
296* Compute workspace needed for DBDSQR
297*
298 bdspac = max( 1, 5*m )
299 minwrk = max( 3*m+nrhs, 3*m+n, bdspac )
300 IF( n.GE.mnthr ) THEN
301*
302* Path 2a - underdetermined, with many more columns
303* than rows
304*
305* Compute space needed for DGELQF
306 CALL dgelqf( m, n, a, lda, dum(1), dum(1),
307 $ -1, info )
308 lwork_dgelqf = int( dum(1) )
309* Compute space needed for DGEBRD
310 CALL dgebrd( m, m, a, lda, s, dum(1), dum(1),
311 $ dum(1), dum(1), -1, info )
312 lwork_dgebrd = int( dum(1) )
313* Compute space needed for DORMBR
314 CALL dormbr( 'Q', 'L', 'T', m, nrhs, n, a, lda,
315 $ dum(1), b, ldb, dum(1), -1, info )
316 lwork_dormbr = int( dum(1) )
317* Compute space needed for DORGBR
318 CALL dorgbr( 'P', m, m, m, a, lda, dum(1),
319 $ dum(1), -1, info )
320 lwork_dorgbr = int( dum(1) )
321* Compute space needed for DORMLQ
322 CALL dormlq( 'L', 'T', n, nrhs, m, a, lda, dum(1),
323 $ b, ldb, dum(1), -1, info )
324 lwork_dormlq = int( dum(1) )
325* Compute total workspace needed
326 maxwrk = m + lwork_dgelqf
327 maxwrk = max( maxwrk, m*m + 4*m + lwork_dgebrd )
328 maxwrk = max( maxwrk, m*m + 4*m + lwork_dormbr )
329 maxwrk = max( maxwrk, m*m + 4*m + lwork_dorgbr )
330 maxwrk = max( maxwrk, m*m + m + bdspac )
331 IF( nrhs.GT.1 ) THEN
332 maxwrk = max( maxwrk, m*m + m + m*nrhs )
333 ELSE
334 maxwrk = max( maxwrk, m*m + 2*m )
335 END IF
336 maxwrk = max( maxwrk, m + lwork_dormlq )
337 ELSE
338*
339* Path 2 - underdetermined
340*
341* Compute space needed for DGEBRD
342 CALL dgebrd( m, n, a, lda, s, dum(1), dum(1),
343 $ dum(1), dum(1), -1, info )
344 lwork_dgebrd = int( dum(1) )
345* Compute space needed for DORMBR
346 CALL dormbr( 'Q', 'L', 'T', m, nrhs, m, a, lda,
347 $ dum(1), b, ldb, dum(1), -1, info )
348 lwork_dormbr = int( dum(1) )
349* Compute space needed for DORGBR
350 CALL dorgbr( 'P', m, n, m, a, lda, dum(1),
351 $ dum(1), -1, info )
352 lwork_dorgbr = int( dum(1) )
353 maxwrk = 3*m + lwork_dgebrd
354 maxwrk = max( maxwrk, 3*m + lwork_dormbr )
355 maxwrk = max( maxwrk, 3*m + lwork_dorgbr )
356 maxwrk = max( maxwrk, bdspac )
357 maxwrk = max( maxwrk, n*nrhs )
358 END IF
359 END IF
360 maxwrk = max( minwrk, maxwrk )
361 END IF
362 work( 1 ) = maxwrk
363*
364 IF( lwork.LT.minwrk .AND. .NOT.lquery )
365 $ info = -12
366 END IF
367*
368 IF( info.NE.0 ) THEN
369 CALL xerbla( 'DGELSS', -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 = dlamch( 'P' )
385 sfmin = dlamch( 'S' )
386 smlnum = sfmin / eps
387 bignum = one / smlnum
388 CALL dlabad( smlnum, bignum )
389*
390* Scale A if max element outside range [SMLNUM,BIGNUM]
391*
392 anrm = dlange( '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 dlascl( '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 dlascl( '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 dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
411 CALL dlaset( 'F', minmn, 1, zero, zero, s, minmn )
412 rank = 0
413 GO TO 70
414 END IF
415*
416* Scale B if max element outside range [SMLNUM,BIGNUM]
417*
418 bnrm = dlange( '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 dlascl( '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 dlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
431 ibscl = 2
432 END IF
433*
434* Overdetermined case
435*
436 IF( m.GE.n ) THEN
437*
438* Path 1 - overdetermined or exactly determined
439*
440 mm = m
441 IF( m.GE.mnthr ) THEN
442*
443* Path 1a - overdetermined, with many more rows than columns
444*
445 mm = n
446 itau = 1
447 iwork = itau + n
448*
449* Compute A=Q*R
450* (Workspace: need 2*N, prefer N+N*NB)
451*
452 CALL dgeqrf( m, n, a, lda, work( itau ), work( iwork ),
453 $ lwork-iwork+1, info )
454*
455* Multiply B by transpose(Q)
456* (Workspace: need N+NRHS, prefer N+NRHS*NB)
457*
458 CALL dormqr( 'L', 'T', m, nrhs, n, a, lda, work( itau ), b,
459 $ ldb, work( iwork ), lwork-iwork+1, info )
460*
461* Zero out below R
462*
463 IF( n.GT.1 )
464 $ CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
465 END IF
466*
467 ie = 1
468 itauq = ie + n
469 itaup = itauq + n
470 iwork = itaup + n
471*
472* Bidiagonalize R in A
473* (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
474*
475 CALL dgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
476 $ work( itaup ), work( iwork ), lwork-iwork+1,
477 $ info )
478*
479* Multiply B by transpose of left bidiagonalizing vectors of R
480* (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
481*
482 CALL dormbr( 'Q', 'L', 'T', mm, nrhs, n, a, lda, work( itauq ),
483 $ b, ldb, work( iwork ), lwork-iwork+1, info )
484*
485* Generate right bidiagonalizing vectors of R in A
486* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
487*
488 CALL dorgbr( 'P', n, n, n, a, lda, work( itaup ),
489 $ work( iwork ), lwork-iwork+1, info )
490 iwork = ie + n
491*
492* Perform bidiagonal QR iteration
493* multiply B by transpose of left singular vectors
494* compute right singular vectors in A
495* (Workspace: need BDSPAC)
496*
497 CALL dbdsqr( 'U', n, n, 0, nrhs, s, work( ie ), a, lda, dum,
498 $ 1, b, ldb, work( iwork ), info )
499 IF( info.NE.0 )
500 $ GO TO 70
501*
502* Multiply B by reciprocals of singular values
503*
504 thr = max( rcond*s( 1 ), sfmin )
505 IF( rcond.LT.zero )
506 $ thr = max( eps*s( 1 ), sfmin )
507 rank = 0
508 DO 10 i = 1, n
509 IF( s( i ).GT.thr ) THEN
510 CALL drscl( nrhs, s( i ), b( i, 1 ), ldb )
511 rank = rank + 1
512 ELSE
513 CALL dlaset( 'F', 1, nrhs, zero, zero, b( i, 1 ), ldb )
514 END IF
515 10 CONTINUE
516*
517* Multiply B by right singular vectors
518* (Workspace: need N, prefer N*NRHS)
519*
520 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 ) THEN
521 CALL dgemm( 'T', 'N', n, nrhs, n, one, a, lda, b, ldb, zero,
522 $ work, ldb )
523 CALL dlacpy( 'G', n, nrhs, work, ldb, b, ldb )
524 ELSE IF( nrhs.GT.1 ) THEN
525 chunk = lwork / n
526 DO 20 i = 1, nrhs, chunk
527 bl = min( nrhs-i+1, chunk )
528 CALL dgemm( 'T', 'N', n, bl, n, one, a, lda, b( 1, i ),
529 $ ldb, zero, work, n )
530 CALL dlacpy( 'G', n, bl, work, n, b( 1, i ), ldb )
531 20 CONTINUE
532 ELSE
533 CALL dgemv( 'T', n, n, one, a, lda, b, 1, zero, work, 1 )
534 CALL dcopy( n, work, 1, b, 1 )
535 END IF
536*
537 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
538 $ max( m, 2*m-4, nrhs, n-3*m ) ) THEN
539*
540* Path 2a - underdetermined, with many more columns than rows
541* and sufficient workspace for an efficient algorithm
542*
543 ldwork = m
544 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
545 $ m*lda+m+m*nrhs ) )ldwork = lda
546 itau = 1
547 iwork = m + 1
548*
549* Compute A=L*Q
550* (Workspace: need 2*M, prefer M+M*NB)
551*
552 CALL dgelqf( m, n, a, lda, work( itau ), work( iwork ),
553 $ lwork-iwork+1, info )
554 il = iwork
555*
556* Copy L to WORK(IL), zeroing out above it
557*
558 CALL dlacpy( 'L', m, m, a, lda, work( il ), ldwork )
559 CALL dlaset( 'U', m-1, m-1, zero, zero, work( il+ldwork ),
560 $ ldwork )
561 ie = il + ldwork*m
562 itauq = ie + m
563 itaup = itauq + m
564 iwork = itaup + m
565*
566* Bidiagonalize L in WORK(IL)
567* (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
568*
569 CALL dgebrd( m, m, work( il ), ldwork, s, work( ie ),
570 $ work( itauq ), work( itaup ), work( iwork ),
571 $ lwork-iwork+1, info )
572*
573* Multiply B by transpose of left bidiagonalizing vectors of L
574* (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
575*
576 CALL dormbr( 'Q', 'L', 'T', m, nrhs, m, work( il ), ldwork,
577 $ work( itauq ), b, ldb, work( iwork ),
578 $ lwork-iwork+1, info )
579*
580* Generate right bidiagonalizing vectors of R in WORK(IL)
581* (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB)
582*
583 CALL dorgbr( 'P', m, m, m, work( il ), ldwork, work( itaup ),
584 $ work( iwork ), lwork-iwork+1, info )
585 iwork = ie + m
586*
587* Perform bidiagonal QR iteration,
588* computing right singular vectors of L in WORK(IL) and
589* multiplying B by transpose of left singular vectors
590* (Workspace: need M*M+M+BDSPAC)
591*
592 CALL dbdsqr( 'U', m, m, 0, nrhs, s, work( ie ), work( il ),
593 $ ldwork, a, lda, b, ldb, work( iwork ), info )
594 IF( info.NE.0 )
595 $ GO TO 70
596*
597* Multiply B by reciprocals of singular values
598*
599 thr = max( rcond*s( 1 ), sfmin )
600 IF( rcond.LT.zero )
601 $ thr = max( eps*s( 1 ), sfmin )
602 rank = 0
603 DO 30 i = 1, m
604 IF( s( i ).GT.thr ) THEN
605 CALL drscl( nrhs, s( i ), b( i, 1 ), ldb )
606 rank = rank + 1
607 ELSE
608 CALL dlaset( 'F', 1, nrhs, zero, zero, b( i, 1 ), ldb )
609 END IF
610 30 CONTINUE
611 iwork = ie
612*
613* Multiply B by right singular vectors of L in WORK(IL)
614* (Workspace: need M*M+2*M, prefer M*M+M+M*NRHS)
615*
616 IF( lwork.GE.ldb*nrhs+iwork-1 .AND. nrhs.GT.1 ) THEN
617 CALL dgemm( 'T', 'N', m, nrhs, m, one, work( il ), ldwork,
618 $ b, ldb, zero, work( iwork ), ldb )
619 CALL dlacpy( 'G', m, nrhs, work( iwork ), ldb, b, ldb )
620 ELSE IF( nrhs.GT.1 ) THEN
621 chunk = ( lwork-iwork+1 ) / m
622 DO 40 i = 1, nrhs, chunk
623 bl = min( nrhs-i+1, chunk )
624 CALL dgemm( 'T', 'N', m, bl, m, one, work( il ), ldwork,
625 $ b( 1, i ), ldb, zero, work( iwork ), m )
626 CALL dlacpy( 'G', m, bl, work( iwork ), m, b( 1, i ),
627 $ ldb )
628 40 CONTINUE
629 ELSE
630 CALL dgemv( 'T', m, m, one, work( il ), ldwork, b( 1, 1 ),
631 $ 1, zero, work( iwork ), 1 )
632 CALL dcopy( m, work( iwork ), 1, b( 1, 1 ), 1 )
633 END IF
634*
635* Zero out below first M rows of B
636*
637 CALL dlaset( 'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
638 iwork = itau + m
639*
640* Multiply transpose(Q) by B
641* (Workspace: need M+NRHS, prefer M+NRHS*NB)
642*
643 CALL dormlq( 'L', 'T', n, nrhs, m, a, lda, work( itau ), b,
644 $ ldb, work( iwork ), lwork-iwork+1, info )
645*
646 ELSE
647*
648* Path 2 - remaining underdetermined cases
649*
650 ie = 1
651 itauq = ie + m
652 itaup = itauq + m
653 iwork = itaup + m
654*
655* Bidiagonalize A
656* (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
657*
658 CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
659 $ work( itaup ), work( iwork ), lwork-iwork+1,
660 $ info )
661*
662* Multiply B by transpose of left bidiagonalizing vectors
663* (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
664*
665 CALL dormbr( 'Q', 'L', 'T', m, nrhs, n, a, lda, work( itauq ),
666 $ b, ldb, work( iwork ), lwork-iwork+1, info )
667*
668* Generate right bidiagonalizing vectors in A
669* (Workspace: need 4*M, prefer 3*M+M*NB)
670*
671 CALL dorgbr( 'P', m, n, m, a, lda, work( itaup ),
672 $ work( iwork ), lwork-iwork+1, info )
673 iwork = ie + m
674*
675* Perform bidiagonal QR iteration,
676* computing right singular vectors of A in A and
677* multiplying B by transpose of left singular vectors
678* (Workspace: need BDSPAC)
679*
680 CALL dbdsqr( 'L', m, n, 0, nrhs, s, work( ie ), a, lda, dum,
681 $ 1, b, ldb, work( iwork ), info )
682 IF( info.NE.0 )
683 $ GO TO 70
684*
685* Multiply B by reciprocals of singular values
686*
687 thr = max( rcond*s( 1 ), sfmin )
688 IF( rcond.LT.zero )
689 $ thr = max( eps*s( 1 ), sfmin )
690 rank = 0
691 DO 50 i = 1, m
692 IF( s( i ).GT.thr ) THEN
693 CALL drscl( nrhs, s( i ), b( i, 1 ), ldb )
694 rank = rank + 1
695 ELSE
696 CALL dlaset( 'F', 1, nrhs, zero, zero, b( i, 1 ), ldb )
697 END IF
698 50 CONTINUE
699*
700* Multiply B by right singular vectors of A
701* (Workspace: need N, prefer N*NRHS)
702*
703 IF( lwork.GE.ldb*nrhs .AND. nrhs.GT.1 ) THEN
704 CALL dgemm( 'T', 'N', n, nrhs, m, one, a, lda, b, ldb, zero,
705 $ work, ldb )
706 CALL dlacpy( 'F', n, nrhs, work, ldb, b, ldb )
707 ELSE IF( nrhs.GT.1 ) THEN
708 chunk = lwork / n
709 DO 60 i = 1, nrhs, chunk
710 bl = min( nrhs-i+1, chunk )
711 CALL dgemm( 'T', 'N', n, bl, m, one, a, lda, b( 1, i ),
712 $ ldb, zero, work, n )
713 CALL dlacpy( 'F', n, bl, work, n, b( 1, i ), ldb )
714 60 CONTINUE
715 ELSE
716 CALL dgemv( 'T', m, n, one, a, lda, b, 1, zero, work, 1 )
717 CALL dcopy( n, work, 1, b, 1 )
718 END IF
719 END IF
720*
721* Undo scaling
722*
723 IF( iascl.EQ.1 ) THEN
724 CALL dlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
725 CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
726 $ info )
727 ELSE IF( iascl.EQ.2 ) THEN
728 CALL dlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
729 CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
730 $ info )
731 END IF
732 IF( ibscl.EQ.1 ) THEN
733 CALL dlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
734 ELSE IF( ibscl.EQ.2 ) THEN
735 CALL dlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
736 END IF
737*
738 70 CONTINUE
739 work( 1 ) = maxwrk
740 RETURN
741*
742* End of DGELSS
743*
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
dlascl
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
dlacpy
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
dlaset
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
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
dbdsqr
subroutine dbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DBDSQR
Definition: dbdsqr.f:241
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
dgemv
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:156
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
dorgbr
subroutine dorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGBR
Definition: dorgbr.f:157
dlange
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:114
dgeqrf
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:146
dgelqf
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
Definition: dgelqf.f:143
dgebrd
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
Definition: dgebrd.f:205
drscl
subroutine drscl(N, SA, SX, INCX)
DRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: drscl.f:84
dormqr
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:167
dormlq
subroutine dormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMLQ
Definition: dormlq.f:167
dormbr
subroutine dormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMBR
Definition: dormbr.f:195
Here is the call graph for this function:
Here is the caller graph for this function: