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


SYNOPSIS

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

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

           REAL           RCOND

           INTEGER        IWORK( * )

           REAL           RWORK( * ), S( * )

           COMPLEX        A( LDA, * ), B( LDB, * ), WORK( * )


PURPOSE

       CGELSD 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   tranformations,  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


               The number of rows of the matrix A. M >= 0.

       N       (input) INTEGER
               The number of columns of the matrix 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/output) COMPLEX 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) COMPLEX 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,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) COMPLEX 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

               or  equal  to  N  or 2 * M + M * NRHS if M is less
               than N, the code will execute correctly.  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.

       RWORK   (workspace) REAL array, dimension at least
               10*N  +  2*N*SMLSIZ  +  8*N*NLVL + 3*SMLSIZ*NRHS +
               (SMLSIZ+1)**2 if M is greater than or equal  to  N
               or  10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*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  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 )

       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