Contents


NAME

     dgelsy - compute the minimum-norm solution to a real  linear
     least squares problem

SYNOPSIS

     SUBROUTINE DGELSY(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
           WORK, LWORK, INFO)

     INTEGER M, N, NRHS, LDA, LDB, RANK, LWORK, INFO
     INTEGER JPVT(*)
     DOUBLE PRECISION RCOND
     DOUBLE PRECISION A(LDA,*), B(LDB,*), WORK(*)

     SUBROUTINE DGELSY_64(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
           WORK, LWORK, INFO)

     INTEGER*8 M, N, NRHS, LDA, LDB, RANK, LWORK, INFO
     INTEGER*8 JPVT(*)
     DOUBLE PRECISION RCOND
     DOUBLE PRECISION A(LDA,*), B(LDB,*), WORK(*)

  F95 INTERFACE
     SUBROUTINE GELSY([M], [N], [NRHS], A, [LDA], B, [LDB], JPVT, RCOND,
            RANK, [WORK], [LWORK], [INFO])

     INTEGER :: M, N, NRHS, LDA, LDB, RANK, LWORK, INFO
     INTEGER, DIMENSION(:) :: JPVT
     REAL(8) :: RCOND
     REAL(8), DIMENSION(:) :: WORK
     REAL(8), DIMENSION(:,:) :: A, B

     SUBROUTINE GELSY_64([M], [N], [NRHS], A, [LDA], B, [LDB], JPVT,
            RCOND, RANK, [WORK], [LWORK], [INFO])

     INTEGER(8) :: M, N, NRHS, LDA, LDB, RANK, LWORK, INFO
     INTEGER(8), DIMENSION(:) :: JPVT
     REAL(8) :: RCOND
     REAL(8), DIMENSION(:) :: WORK
     REAL(8), DIMENSION(:,:) :: A, B

  C INTERFACE
     #include <sunperf.h>

     void dgelsy(int m, int n, int nrhs, double *a, int lda, dou-
               ble  *b,  int  ldb,  int  *jpvt, double rcond, int
               *rank, int *info);
     void dgelsy_64(long m, long n, long nrhs,  double  *a,  long
               lda,  double  *b,  long  ldb,  long  *jpvt, double
               rcond, long *rank, long *info);

PURPOSE

     dgelsy computes the minimum-norm solution to a  real  linear
     least squares problem:
         minimize || A * X - B ||
     using a complete orthogonal factorization 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 routine first computes a QR  factorization  with  column
     pivoting:
         A * P = Q * [ R11 R12 ]
                     [  0  R22 ]
     with R11 defined as  the  largest  leading  submatrix  whose
     estimated  condition number is less than 1/RCOND.  The order
     of R11, RANK, is the effective rank of A.

     Then, R22 is considered to be negligible, and R12 is annihi-
     lated by orthogonal transformations from the right, arriving
     at the complete orthogonal factorization:
        A * P = Q * [ T11 0 ] * Z
                    [  0  0 ]
     The minimum-norm solution is then
        X = P * Z' [ inv(T11)*Q1'*B ]
                   [        0       ]
     where Q1 consists of the first RANK columns of Q.

     This routine is basically identical to the  original  xGELSX
     except three differences:
       o The call to the subroutine xGEQPF has  been  substituted
     by the
         the call to the subroutine xGEQP3. This subroutine is  a
     Blas-3
         version of the QR factorization with column pivoting.
       o Matrix B (the right hand side) is updated with Blas-3.
       o The permutation of matrix B (the  right  hand  side)  is
     faster and
         more simple.

ARGUMENTS

     M (input) The number of rows of the matrix A.  M >= 0.
     N (input) The number of columns of the matrix A.  N >= 0.

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

     A (input/output)
               On entry, the M-by-N matrix A.   On  exit,  A  has
               been   overwritten  by  details  of  its  complete
               orthogonal factorization.

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

     B (input/output)
               On entry, the M-by-NRHS right hand side matrix  B.
               On exit, the N-by-NRHS solution matrix X.

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

     JPVT (input/output)
               On entry, if JPVT(i) .ne. 0, the i-th column of  A
               is permuted to the front of AP, otherwise column i
               is a free column.  On exit, if JPVT(i) =  k,  then
               the i-th column of AP was the k-th column of A.

     RCOND (input)
               RCOND is used to determine the effective  rank  of
               A,  which  is  defined as the order of the largest
               leading triangular submatrix R11 in the QR factor-
               ization with pivoting of A, whose estimated condi-
               tion number < 1/RCOND.

     RANK (output)
               The effective rank of A, i.e., the  order  of  the
               submatrix  R11.   This is the same as the order of
               the submatrix T11 in the complete orthogonal  fac-
               torization of A.

     WORK (workspace)
               On exit, if INFO = 0, WORK(1) returns the  optimal
               LWORK.

     LWORK (input)
               The dimension of the array  WORK.   The  unblocked
               strategy  requires  that:  LWORK >= MAX( MN+3*N+1,
               2*MN+NRHS ), where MN = min( M, N  ).   The  block
               algorithm    requires   that:    LWORK   >=   MAX(
               MN+2*N+NB*(N+1), 2*MN+NB*NRHS ), where  NB  is  an
               upper  bound  on  the blocksize returned by ILAENV
               for the routines SGEQP3, STZRZF,  STZRQF,  SORMQR,
               and SORMRZ.

               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.

     INFO (output)
               = 0: successful exit
               < 0: If INFO = -i, the i-th argument had an  ille-
               gal value.

FURTHER DETAILS

     Based on contributions by
       A. Petitet, Computer Science Dept., Univ. of Tenn.,  Knox-
     ville, USA
       E. Quintana-Orti, Depto. de Informatica, Universidad Jaime
     I, Spain
       G. Quintana-Orti, Depto. de Informatica, Universidad Jaime
     I, Spain