LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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,
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=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=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=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=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=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=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=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=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=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=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=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=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=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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
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
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DBDSQR
Definition: dbdsqr.f:241
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:156
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
subroutine dorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGBR
Definition: dorgbr.f:157
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
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:145
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 drscl(N, SA, SX, INCX)
DRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: drscl.f:84
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 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: