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


SYNOPSIS

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

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

           DOUBLE         PRECISION RCOND

           INTEGER        IWORK( * )

           DOUBLE         PRECISION A( LDA, * ), B( LDB, * ),  S(
                          * ), WORK( * )


PURPOSE

       DGELSD 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) DOUBLE PRECISION 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)   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 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) 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)).

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

       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)  DOUBLE  PRECISION  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

               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 exe­
               cute  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 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