sgelsd()

subroutine sgelsd	(	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, dimension( * )	iwork,
		integer	info
	)

SGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices

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

Purpose:

 SGELSD 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 problem is solved in three steps:
 (1) Reduce the coefficient matrix A to bidiagonal form with
     Householder transformations, reducing the original problem
     into a "bidiagonal least squares problem" (BLS)
 (2) Solve the BLS using a divide and conquer approach.
 (3) Apply back all the Householder transformations to solve
     the original least squares problem.

 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 A. M >= 0.
[in]	N	N is INTEGER The number of columns of 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, A has been destroyed.
[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 must be at least 1. The exact minimum amount of workspace needed depends on M, N and NRHS. As long as LWORK is at least 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2, if M is greater than or equal to N or 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2, if M is less than N, the code will execute correctly. SMLSIZ is returned by ILAENV and is equal to the maximum size of the subproblems at the bottom of the computation tree (usually about 25), and NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 ) 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 array WORK and the minimum size of the array IWORK, and returns these values as the first entries of the WORK and IWORK arrays, and no error message related to LWORK is issued by XERBLA.
[out]	IWORK	IWORK is INTEGER array, dimension (MAX(1,LIWORK)) LIWORK >= max(1, 3*MINMN*NLVL + 11*MINMN), where MINMN = MIN( M,N ). On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
[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.

Contributors:

Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

Definition at line 202 of file sgelsd.f.

*
*  -- LAPACK driver routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
      REAL               RCOND
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      REAL               A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE, TWO
      parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY
      INTEGER            IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
     $                   LDWORK, LIWORK, MAXMN, MAXWRK, MINMN, MINWRK,
     $                   MM, MNTHR, NLVL, NWORK, SMLSIZ, WLALSD
      REAL               ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
*     ..
*     .. External Subroutines ..
      EXTERNAL           sgebrd, sgelqf, sgeqrf, slacpy, slalsd, slascl,
     $                   slaset, sormbr, sormlq, sormqr, xerbla
*     ..
*     .. External Functions ..
      INTEGER            ILAENV
      REAL               SLAMCH, SLANGE, SROUNDUP_LWORK
      EXTERNAL           slamch, slange, ilaenv, sroundup_lwork
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          int, log, max, min, real
*     ..
*     .. Executable Statements ..
*
*     Test the input arguments.
*
      info = 0
      minmn = min( m, n )
      maxmn = max( m, n )
      lquery = ( lwork.EQ.-1 )
      IF( m.LT.0 ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE IF( nrhs.LT.0 ) THEN
         info = -3
      ELSE IF( lda.LT.max( 1, m ) ) THEN
         info = -5
      ELSE IF( ldb.LT.max( 1, maxmn ) ) THEN
         info = -7
      END IF
*
*     Compute workspace.
*     (Note: Comments in the code beginning "Workspace:" describe the
*     minimal amount of workspace needed at that point in the code,
*     as well as the preferred amount for good performance.
*     NB refers to the optimal block size for the immediately
*     following subroutine, as returned by ILAENV.)
*
      IF( info.EQ.0 ) THEN
         minwrk = 1
         maxwrk = 1
         liwork = 1
         IF( minmn.GT.0 ) THEN
            smlsiz = ilaenv( 9, 'SGELSD', ' ', 0, 0, 0, 0 )
            mnthr = ilaenv( 6, 'SGELSD', ' ', m, n, nrhs, -1 )
            nlvl = max( int( log( real( minmn ) / real( smlsiz + 1 ) ) /
     $                  log( two ) ) + 1, 0 )
            liwork = 3*minmn*nlvl + 11*minmn
            mm = m
            IF( m.GE.n .AND. m.GE.mnthr ) THEN
*
*              Path 1a - overdetermined, with many more rows than
*                        columns.
*
               mm = n
               maxwrk = max( maxwrk, n + n*ilaenv( 1, 'SGEQRF', ' ', m,
     $                       n, -1, -1 ) )
               maxwrk = max( maxwrk, n + nrhs*ilaenv( 1, 'SORMQR', 'LT',
     $                       m, nrhs, n, -1 ) )
            END IF
            IF( m.GE.n ) THEN
*
*              Path 1 - overdetermined or exactly determined.
*
               maxwrk = max( maxwrk, 3*n + ( mm + n )*ilaenv( 1,
     $                       'SGEBRD', ' ', mm, n, -1, -1 ) )
               maxwrk = max( maxwrk, 3*n + nrhs*ilaenv( 1, 'SORMBR',
     $                       'QLT', mm, nrhs, n, -1 ) )
               maxwrk = max( maxwrk, 3*n + ( n - 1 )*ilaenv( 1,
     $                       'SORMBR', 'PLN', n, nrhs, n, -1 ) )
               wlalsd = 9*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs +
     $                  ( smlsiz + 1 )**2
               maxwrk = max( maxwrk, 3*n + wlalsd )
               minwrk = max( 3*n + mm, 3*n + nrhs, 3*n + wlalsd )
            END IF
            IF( n.GT.m ) THEN
               wlalsd = 9*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs +
     $                  ( smlsiz + 1 )**2
               IF( n.GE.mnthr ) THEN
*
*                 Path 2a - underdetermined, with many more columns
*                           than rows.
*
                  maxwrk = m + m*ilaenv( 1, 'SGELQF', ' ', m, n, -1,
     $                                  -1 )
                  maxwrk = max( maxwrk, m*m + 4*m + 2*m*ilaenv( 1,
     $                          'SGEBRD', ' ', m, m, -1, -1 ) )
                  maxwrk = max( maxwrk, m*m + 4*m + nrhs*ilaenv( 1,
     $                          'SORMBR', 'QLT', m, nrhs, m, -1 ) )
                  maxwrk = max( maxwrk, m*m + 4*m + ( m - 1 )*ilaenv( 1,
     $                          'SORMBR', 'PLN', m, nrhs, m, -1 ) )
                  IF( nrhs.GT.1 ) THEN
                     maxwrk = max( maxwrk, m*m + m + m*nrhs )
                  ELSE
                     maxwrk = max( maxwrk, m*m + 2*m )
                  END IF
                  maxwrk = max( maxwrk, m + nrhs*ilaenv( 1, 'SORMLQ',
     $                          'LT', n, nrhs, m, -1 ) )
                  maxwrk = max( maxwrk, m*m + 4*m + wlalsd )
!     XXX: Ensure the Path 2a case below is triggered.  The workspace
!     calculation should use queries for all routines eventually.
                  maxwrk = max( maxwrk,
     $                 4*m+m*m+max( m, 2*m-4, nrhs, n-3*m ) )
               ELSE
*
*                 Path 2 - remaining underdetermined cases.
*
                  maxwrk = 3*m + ( n + m )*ilaenv( 1, 'SGEBRD', ' ', m,
     $                     n, -1, -1 )
                  maxwrk = max( maxwrk, 3*m + nrhs*ilaenv( 1, 'SORMBR',
     $                          'QLT', m, nrhs, n, -1 ) )
                  maxwrk = max( maxwrk, 3*m + m*ilaenv( 1, 'SORMBR',
     $                          'PLN', n, nrhs, m, -1 ) )
                  maxwrk = max( maxwrk, 3*m + wlalsd )
               END IF
               minwrk = max( 3*m + nrhs, 3*m + m, 3*m + wlalsd )
            END IF
         END IF
         minwrk = min( minwrk, maxwrk )
         work( 1 ) = sroundup_lwork(maxwrk)
         iwork( 1 ) = liwork
*
         IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
            info = -12
         END IF
      END IF
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'SGELSD', -info )
         RETURN
      ELSE IF( lquery ) THEN
         RETURN
      END IF
*
*     Quick return if possible.
*
      IF( m.EQ.0 .OR. n.EQ.0 ) THEN
         rank = 0
         RETURN
      END IF
*
*     Get machine parameters.
*
      eps = slamch( 'P' )
      sfmin = slamch( 'S' )
      smlnum = sfmin / eps
      bignum = one / smlnum
*
*     Scale A if max entry outside range [SMLNUM,BIGNUM].
*
      anrm = slange( 'M', m, n, a, lda, work )
      iascl = 0
      IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
*
*        Scale matrix norm up to SMLNUM.
*
         CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
         iascl = 1
      ELSE IF( anrm.GT.bignum ) THEN
*
*        Scale matrix norm down to BIGNUM.
*
         CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
         iascl = 2
      ELSE IF( anrm.EQ.zero ) THEN
*
*        Matrix all zero. Return zero solution.
*
         CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
         CALL slaset( 'F', minmn, 1, zero, zero, s, 1 )
         rank = 0
         GO TO 10
      END IF
*
*     Scale B if max entry outside range [SMLNUM,BIGNUM].
*
      bnrm = slange( 'M', m, nrhs, b, ldb, work )
      ibscl = 0
      IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
*
*        Scale matrix norm up to SMLNUM.
*
         CALL slascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
         ibscl = 1
      ELSE IF( bnrm.GT.bignum ) THEN
*
*        Scale matrix norm down to BIGNUM.
*
         CALL slascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
         ibscl = 2
      END IF
*
*     If M < N make sure certain entries of B are zero.
*
      IF( m.LT.n )
     $   CALL slaset( 'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
*
*     Overdetermined case.
*
      IF( m.GE.n ) THEN
*
*        Path 1 - overdetermined or exactly determined.
*
         mm = m
         IF( m.GE.mnthr ) THEN
*
*           Path 1a - overdetermined, with many more rows than columns.
*
            mm = n
            itau = 1
            nwork = itau + n
*
*           Compute A=Q*R.
*           (Workspace: need 2*N, prefer N+N*NB)
*
            CALL sgeqrf( m, n, a, lda, work( itau ), work( nwork ),
     $                   lwork-nwork+1, info )
*
*           Multiply B by transpose(Q).
*           (Workspace: need N+NRHS, prefer N+NRHS*NB)
*
            CALL sormqr( 'L', 'T', m, nrhs, n, a, lda, work( itau ), b,
     $                   ldb, work( nwork ), lwork-nwork+1, info )
*
*           Zero out below R.
*
            IF( n.GT.1 ) THEN
               CALL slaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
            END IF
         END IF
*
         ie = 1
         itauq = ie + n
         itaup = itauq + n
         nwork = itaup + n
*
*        Bidiagonalize R in A.
*        (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
*
         CALL sgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
     $                work( itaup ), work( nwork ), lwork-nwork+1,
     $                info )
*
*        Multiply B by transpose of left bidiagonalizing vectors of R.
*        (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
*
         CALL sormbr( 'Q', 'L', 'T', mm, nrhs, n, a, lda, work( itauq ),
     $                b, ldb, work( nwork ), lwork-nwork+1, info )
*
*        Solve the bidiagonal least squares problem.
*
         CALL slalsd( 'U', smlsiz, n, nrhs, s, work( ie ), b, ldb,
     $                rcond, rank, work( nwork ), iwork, info )
         IF( info.NE.0 ) THEN
            GO TO 10
         END IF
*
*        Multiply B by right bidiagonalizing vectors of R.
*
         CALL sormbr( 'P', 'L', 'N', n, nrhs, n, a, lda, work( itaup ),
     $                b, ldb, work( nwork ), lwork-nwork+1, info )
*
      ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
     $         max( m, 2*m-4, nrhs, n-3*m, wlalsd ) ) THEN
*
*        Path 2a - underdetermined, with many more columns than rows
*        and sufficient workspace for an efficient algorithm.
*
         ldwork = m
         IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
     $       m*lda+m+m*nrhs, 4*m+m*lda+wlalsd ) )ldwork = lda
         itau = 1
         nwork = m + 1
*
*        Compute A=L*Q.
*        (Workspace: need 2*M, prefer M+M*NB)
*
         CALL sgelqf( m, n, a, lda, work( itau ), work( nwork ),
     $                lwork-nwork+1, info )
         il = nwork
*
*        Copy L to WORK(IL), zeroing out above its diagonal.
*
         CALL slacpy( 'L', m, m, a, lda, work( il ), ldwork )
         CALL slaset( 'U', m-1, m-1, zero, zero, work( il+ldwork ),
     $                ldwork )
         ie = il + ldwork*m
         itauq = ie + m
         itaup = itauq + m
         nwork = itaup + m
*
*        Bidiagonalize L in WORK(IL).
*        (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
*
         CALL sgebrd( m, m, work( il ), ldwork, s, work( ie ),
     $                work( itauq ), work( itaup ), work( nwork ),
     $                lwork-nwork+1, info )
*
*        Multiply B by transpose of left bidiagonalizing vectors of L.
*        (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
*
         CALL sormbr( 'Q', 'L', 'T', m, nrhs, m, work( il ), ldwork,
     $                work( itauq ), b, ldb, work( nwork ),
     $                lwork-nwork+1, info )
*
*        Solve the bidiagonal least squares problem.
*
         CALL slalsd( 'U', smlsiz, m, nrhs, s, work( ie ), b, ldb,
     $                rcond, rank, work( nwork ), iwork, info )
         IF( info.NE.0 ) THEN
            GO TO 10
         END IF
*
*        Multiply B by right bidiagonalizing vectors of L.
*
         CALL sormbr( 'P', 'L', 'N', m, nrhs, m, work( il ), ldwork,
     $                work( itaup ), b, ldb, work( nwork ),
     $                lwork-nwork+1, info )
*
*        Zero out below first M rows of B.
*
         CALL slaset( 'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
         nwork = itau + m
*
*        Multiply transpose(Q) by B.
*        (Workspace: need M+NRHS, prefer M+NRHS*NB)
*
         CALL sormlq( 'L', 'T', n, nrhs, m, a, lda, work( itau ), b,
     $                ldb, work( nwork ), lwork-nwork+1, info )
*
      ELSE
*
*        Path 2 - remaining underdetermined cases.
*
         ie = 1
         itauq = ie + m
         itaup = itauq + m
         nwork = itaup + m
*
*        Bidiagonalize A.
*        (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
*
         CALL sgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
     $                work( itaup ), work( nwork ), lwork-nwork+1,
     $                info )
*
*        Multiply B by transpose of left bidiagonalizing vectors.
*        (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
*
         CALL sormbr( 'Q', 'L', 'T', m, nrhs, n, a, lda, work( itauq ),
     $                b, ldb, work( nwork ), lwork-nwork+1, info )
*
*        Solve the bidiagonal least squares problem.
*
         CALL slalsd( 'L', smlsiz, m, nrhs, s, work( ie ), b, ldb,
     $                rcond, rank, work( nwork ), iwork, info )
         IF( info.NE.0 ) THEN
            GO TO 10
         END IF
*
*        Multiply B by right bidiagonalizing vectors of A.
*
         CALL sormbr( 'P', 'L', 'N', n, nrhs, m, a, lda, work( itaup ),
     $                b, ldb, work( nwork ), lwork-nwork+1, info )
*
      END IF
*
*     Undo scaling.
*
      IF( iascl.EQ.1 ) THEN
         CALL slascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
         CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
     $                info )
      ELSE IF( iascl.EQ.2 ) THEN
         CALL slascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
         CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
     $                info )
      END IF
      IF( ibscl.EQ.1 ) THEN
         CALL slascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
      ELSE IF( ibscl.EQ.2 ) THEN
         CALL slascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
      END IF
*
   10 CONTINUE
      work( 1 ) = sroundup_lwork(maxwrk)
      iwork( 1 ) = liwork
      RETURN
*
*     End of SGELSD
*

Here is the call graph for this function:

Here is the caller graph for this function: