LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
November 2011

Definition at line 174 of file dgelss.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: