*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