LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ sgelss()

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.

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