*DECK DGLSS
SUBROUTINE DGLSS (A, MDA, M, N, B, MDB, NB, RNORM, WORK, LW,
+ IWORK, LIW, INFO)
C***BEGIN PROLOGUE DGLSS
C***PURPOSE Solve a linear least squares problems by performing a QR
C factorization of the input matrix using Householder
C transformations. Emphasis is put on detecting possible
C rank deficiency.
C***LIBRARY SLATEC
C***CATEGORY D9, D5
C***TYPE DOUBLE PRECISION (SGLSS-S, DGLSS-D)
C***KEYWORDS LINEAR LEAST SQUARES, LQ FACTORIZATION, QR FACTORIZATION,
C UNDERDETERMINED LINEAR SYSTEMS
C***AUTHOR Manteuffel, T. A., (LANL)
C***DESCRIPTION
C
C DGLSS solves both underdetermined and overdetermined
C LINEAR systems AX = B, where A is an M by N matrix
C and B is an M by NB matrix of right hand sides. If
C M.GE.N, the least squares solution is computed by
C decomposing the matrix A into the product of an
C orthogonal matrix Q and an upper triangular matrix
C R (QR factorization). If M.LT.N, the minimal
C length solution is computed by factoring the
C matrix A into the product of a lower triangular
C matrix L and an orthogonal matrix Q (LQ factor-
C ization). If the matrix A is determined to be rank
C deficient, that is the rank of A is less than
C MIN(M,N), then the minimal length least squares
C solution is computed.
C
C DGLSS assumes full machine precision in the data.
C If more control over the uncertainty in the data
C is desired, the codes DLLSIA and DULSIA are
C recommended.
C
C DGLSS requires MDA*N + (MDB + 1)*NB + 5*MIN(M,N) dimensioned
C real space and M+N dimensioned integer space.
C
C
C ******************************************************************
C * *
C * WARNING - All input arrays are changed on exit. *
C * *
C ******************************************************************
C SUBROUTINE DGLSS(A,MDA,M,N,B,MDB,NB,RNORM,WORK,LW,IWORK,LIW,INFO)
C
C Input..All TYPE REAL variables are DOUBLE PRECISION
C
C A(,) Linear coefficient matrix of AX=B, with MDA the
C MDA,M,N actual first dimension of A in the calling program.
C M is the row dimension (no. of EQUATIONS of the
C problem) and N the col dimension (no. of UNKNOWNS).
C
C B(,) Right hand side(s), with MDB the actual first
C MDB,NB dimension of B in the calling program. NB is the
C number of M by 1 right hand sides. Must have
C MDB.GE.MAX(M,N). If NB = 0, B is never accessed.
C
C
C RNORM() Vector of length at least NB. On input the contents
C of RNORM are unused.
C
C WORK() A real work array dimensioned 5*MIN(M,N).
C
C LW Actual dimension of WORK.
C
C IWORK() Integer work array dimensioned at least N+M.
C
C LIW Actual dimension of IWORK.
C
C
C INFO A flag which provides for the efficient
C solution of subsequent problems involving the
C same A but different B.
C If INFO = 0 original call
C INFO = 1 subsequent calls
C On subsequent calls, the user must supply A, INFO,
C LW, IWORK, LIW, and the first 2*MIN(M,N) locations
C of WORK as output by the original call to DGLSS.
C
C
C Output..All TYPE REAL variables are DOUBLE PRECISION
C
C A(,) Contains the triangular part of the reduced matrix
C and the transformation information. It together with
C the first 2*MIN(M,N) elements of WORK (see below)
C completely specify the factorization of A.
C
C B(,) Contains the N by NB solution matrix X.
C
C
C RNORM() Contains the Euclidean length of the NB residual
C vectors B(I)-AX(I), I=1,NB.
C
C WORK() The first 2*MIN(M,N) locations of WORK contain value
C necessary to reproduce the factorization of A.
C
C IWORK() The first M+N locations contain the order in
C which the rows and columns of A were used.
C If M.GE.N columns then rows. If M.LT.N rows
C then columns.
C
C INFO Flag to indicate status of computation on completion
C -1 Parameter error(s)
C 0 - Full rank
C N.GT.0 - Reduced rank rank=MIN(M,N)-INFO
C
C***REFERENCES T. Manteuffel, An interval analysis approach to rank
C determination in linear least squares problems,
C Report SAND80-0655, Sandia Laboratories, June 1980.
C***ROUTINES CALLED DLLSIA, DULSIA
C***REVISION HISTORY (YYMMDD)
C 810801 DATE WRITTEN
C 890831 Modified array declarations. (WRB)
C 891006 Cosmetic changes to prologue. (WRB)
C 891006 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE DGLSS
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION A(MDA,*),B(MDB,*),RNORM(*),WORK(*)
INTEGER IWORK(*)
C
C***FIRST EXECUTABLE STATEMENT DGLSS
RE=0.D0
AE=0.D0
KEY=0
MODE=2
NP=0
C
C IF M.GE.N CALL DLLSIA
C IF M.LT.N CALL DULSIA
C
IF(M.LT.N) GO TO 10
CALL DLLSIA(A,MDA,M,N,B,MDB,NB,RE,AE,KEY,MODE,NP,
1 KRANK,KSURE,RNORM,WORK,LW,IWORK,LIW,INFO)
IF(INFO.EQ.-1) RETURN
INFO=N-KRANK
RETURN
10 CALL DULSIA(A,MDA,M,N,B,MDB,NB,RE,AE,KEY,MODE,NP,
1 KRANK,KSURE,RNORM,WORK,LW,IWORK,LIW,INFO)
IF(INFO.EQ.-1) RETURN
INFO=M-KRANK
RETURN
END