LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine sgelss ( integer  M,
integer  N,
integer  NRHS,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldb, * )  B,
integer  LDB,
real, dimension( * )  S,
real  RCOND,
integer  RANK,
real, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

SGELSS solves overdetermined or underdetermined systems for GE matrices

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

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

Here is the call graph for this function:

Here is the caller graph for this function: