SGELSD  - compute the minimum-norm solution to a real lin­
       ear least squares problem


SYNOPSIS

       SUBROUTINE SGELSD( M, N, NRHS, A, LDA, B, LDB,  S,  RCOND,
                          RANK, WORK, LWORK, IWORK, INFO )

           INTEGER        INFO, LDA, LDB, LWORK, M, N, NRHS, RANK

           REAL           RCOND

           INTEGER        IWORK( * )

           REAL           A( LDA, * ), B( LDB, * ), S( * ), WORK(
                          * )


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 tranformations 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.

       The  divide  and conquer algorithm makes very mild assump­
       tions about floating point arithmetic.  It  will  work  on
       machines  with  a guard digit in add/subtract, or on those
       binary machines without guard digits which  subtract  like
       the  Cray  X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could
       conceivably fail on hexadecimal or decimal machines  with­
       out guard digits, but we know of none.


ARGUMENTS

       M       (input) INTEGER
               The number of rows of A. M >= 0.

               The number of columns of A. N >= 0.

       NRHS    (input) INTEGER
               The  number  of right hand sides, i.e., the number
               of columns of the matrices B and X. NRHS >= 0.

       A       (input) REAL array, dimension (LDA,N)
               On entry, the M-by-N matrix A.   On  exit,  A  has
               been destroyed.

       LDA     (input) INTEGER
               The  leading  dimension  of  the  array A.  LDA >=
               max(1,M).

       B       (input/output) 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 solu­
               tion matrix X.  If m >= n and RANK = n, the resid­
               ual  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.

       LDB     (input) INTEGER
               The  leading  dimension  of  the  array  B. LDB >=
               max(1,max(M,N)).

       S       (output) 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)).

       RCOND   (input) 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.

       RANK    (output) INTEGER
               The  effective rank of A, i.e., the number of sin­
               gular values which are greater than RCOND*S(1).

       WORK    (workspace/output) REAL array, dimension (LWORK)
               On exit, if INFO = 0, WORK(1) returns the  optimal
               LWORK.

       LWORK   (input) 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 + (SML­
               SIZ+1)**2,  if  M  is  less  than N, the code will

               and  is  equal to the maximum size of the subprob­
               lems at the bottom of the computation  tree  (usu­
               ally  about  25),  and  NLVL = MAX( 0, INT( LOG_2(
               MIN( M,N )/(SMLSIZ+1) ) ) + 1 ) For  good  perfor­
               mance, 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.

       IWORK   (workspace) INTEGER array, dimension (LIWORK)
               LIWORK  >=  3  *  MINMN * NLVL + 11 * MINMN, where
               MINMN = MIN( M,N ).

       INFO    (output) INTEGER
               = 0:  successful exit
               < 0:  if INFO = -i, the i-th argument had an ille­
               gal 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 con­
               verge to zero.


FURTHER DETAILS

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


Man(1) output converted with man2html