C
C     ..................................................................
C
C        SUBROUTINE DLLSQ
C
C        PURPOSE
C           TO SOLVE LINEAR LEAST SQUARES PROBLEMS, I.E. TO MINIMIZE
C           THE EUCLIDEAN NORM OF B-A*X, WHERE A IS A M BY N MATRIX
C           WITH M NOT LESS THAN N. IN THE SPECIAL CASE M=N SYSTEMS OF
C           LINEAR EQUATIONS MAY BE SOLVED.
C
C        USAGE
C           CALL DLLSQ (A,B,M,N,L,X,IPIV,EPS,IER,AUX)
C
C        DESCRIPTION OF PARAMETERS
C           A      - DOUBLE PRECISION M BY N COEFFICIENT MATRIX
C                    (DESTROYED).
C           B      - DOUBLE PRECISION M BY L RIGHT HAND SIDE MATRIX
C                    (DESTROYED).
C           M      - ROW NUMBER OF MATRICES A AND B.
C           N      - COLUMN NUMBER OF MATRIX A, ROW NUMBER OF MATRIX X.
C           L      - COLUMN NUMBER OF MATRICES B AND X.
C           X      - DOUBLE PRECISION N BY L SOLUTION MATRIX.
C           IPIV   - INTEGER OUTPUT VECTOR OF DIMENSION N WHICH
C                    CONTAINS INFORMATIONS ON COLUMN INTERCHANGES
C                    IN MATRIX A. (SEE REMARK NO.3).
C           EPS    - SINGLE PRECISION INPUT PARAMETER WHICH SPECIFIES
C                    A RELATIVE TOLERANCE FOR DETERMINATION OF RANK OF
C                    MATRIX A.
C           IER    - A RESULTING ERROR PARAMETER.
C           AUX    - A DOUBLE PRECISION AUXILIARY STORAGE ARRAY OF
C                    DIMENSION MAX(2*N,L). ON RETURN FIRST L LOCATIONS
C                    OF AUX CONTAIN THE RESULTING LEAST SQUARES.
C
C        REMARKS
C           (1) NO ACTION BESIDES ERROR MESSAGE IER=-2 IN CASE
C               M LESS THAN N.
C           (2) NO ACTION BESIDES ERROR MESSAGE IER=-1 IN CASE
C               OF A ZERO-MATRIX A.
C           (3) IF RANK K OF MATRIX A IS FOUND TO BE LESS THAN N BUT
C               GREATER THAN 0, THE PROCEDURE RETURNS WITH ERROR CODE
C               IER=K INTO CALLING PROGRAM. THE LAST N-K ELEMENTS OF
C               VECTOR IPIV DENOTE THE USELESS COLUMNS IN MATRIX A.
C               THE REMAINING USEFUL COLUMNS FORM A BASE OF MATRIX A.
C           (4) IF THE PROCEDURE WAS SUCCESSFUL, ERROR PARAMETER IER
C               IS SET TO 0.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           HOUSEHOLDER TRANSFORMATIONS ARE USED TO TRANSFORM MATRIX A
C           TO UPPER TRIANGULAR FORM. AFTER HAVING APPLIED THE SAME
C           TRANSFORMATION TO THE RIGHT HAND SIDE MATRIX B, AN
C           APPROXIMATE SOLUTION OF THE PROBLEM IS COMPUTED BY
C           BACK SUBSTITUTION. FOR REFERENCE, SEE
C           G. GOLUB, NUMERICAL METHODS FOR SOLVING LINEAR LEAST
C           SQUARES PROBLEMS, NUMERISCHE MATHEMATIK, VOL.7,
C           ISS.3 (1965), PP.206-216.
C
C     ..................................................................
C
      SUBROUTINE DLLSQ(A,B,M,N,L,X,IPIV,EPS,IER,AUX)
C
      DIMENSION A(1),B(1),X(1),IPIV(1),AUX(1)
      DOUBLE PRECISION A,B,X,AUX,PIV,H,SIG,BETA,TOL
C
C     ERROR TEST
      IF(M-N)30,1,1
C
C     GENERATION OF INITIAL VECTOR S(K) (K=1,2,...,N) IN STORAGE
C     LOCATIONS AUX(K) (K=1,2,...,N)
    1 PIV=0.D0
      IEND=0
      DO 4 K=1,N
      IPIV(K)=K
      H=0.D0
      IST=IEND+1
      IEND=IEND+M
      DO 2 I=IST,IEND
    2 H=H+A(I)*A(I)
      AUX(K)=H
      IF(H-PIV)4,4,3
    3 PIV=H
      KPIV=K
    4 CONTINUE
C
C     ERROR TEST
      IF(PIV)31,31,5
C
C     DEFINE TOLERANCE FOR CHECKING RANK OF A
    5 SIG=DSQRT(PIV)
      TOL=SIG*ABS(EPS)
C
C
C     DECOMPOSITION LOOP
      LM=L*M
      IST=-M
      DO 21 K=1,N
      IST=IST+M+1
      IEND=IST+M-K
      I=KPIV-K
      IF(I)8,8,6
C
C     INTERCHANGE K-TH COLUMN OF A WITH KPIV-TH IN CASE KPIV.GT.K
    6 H=AUX(K)
      AUX(K)=AUX(KPIV)
      AUX(KPIV)=H
      ID=I*M
      DO 7 I=IST,IEND
      J=I+ID
      H=A(I)
      A(I)=A(J)
    7 A(J)=H
C
C     COMPUTATION OF PARAMETER SIG
    8 IF(K-1)11,11,9
    9 SIG=0.D0
      DO 10 I=IST,IEND
   10 SIG=SIG+A(I)*A(I)
      SIG=DSQRT(SIG)
C
C     TEST ON SINGULARITY
      IF(SIG-TOL)32,32,11
C
C     GENERATE CORRECT SIGN OF PARAMETER SIG
   11 H=A(IST)
      IF(H)12,13,13
   12 SIG=-SIG
C
C     SAVE INTERCHANGE INFORMATION
   13 IPIV(KPIV)=IPIV(K)
      IPIV(K)=KPIV
C
C     GENERATION OF VECTOR UK IN K-TH COLUMN OF MATRIX A AND OF
C     PARAMETER BETA
      BETA=H+SIG
      A(IST)=BETA
      BETA=1.D0/(SIG*BETA)
      J=N+K
      AUX(J)=-SIG
      IF(K-N)14,19,19
C
C     TRANSFORMATION OF MATRIX A
   14 PIV=0.D0
      ID=0
      JST=K+1
      KPIV=JST
      DO 18 J=JST,N
      ID=ID+M
      H=0.D0
      DO 15 I=IST,IEND
      II=I+ID
   15 H=H+A(I)*A(II)
      H=BETA*H
      DO 16 I=IST,IEND
      II=I+ID
   16 A(II)=A(II)-A(I)*H
C
C     UPDATING OF ELEMENT S(J) STORED IN LOCATION AUX(J)
      II=IST+ID
      H=AUX(J)-A(II)*A(II)
      AUX(J)=H
      IF(H-PIV)18,18,17
   17 PIV=H
      KPIV=J
   18 CONTINUE
C
C     TRANSFORMATION OF RIGHT HAND SIDE MATRIX B
   19 DO 21 J=K,LM,M
      H=0.D0
      IEND=J+M-K
      II=IST
      DO 20 I=J,IEND
      H=H+A(II)*B(I)
   20 II=II+1
      H=BETA*H
      II=IST
      DO 21 I=J,IEND
      B(I)=B(I)-A(II)*H
   21 II=II+1
C     END OF DECOMPOSITION LOOP
C
C
C     BACK SUBSTITUTION AND BACK INTERCHANGE
      IER=0
      I=N
      LN=L*N
      PIV=1.D0/AUX(2*N)
      DO 22 K=N,LN,N
      X(K)=PIV*B(I)
   22 I=I+M
      IF(N-1)26,26,23
   23 JST=(N-1)*M+N
      DO 25 J=2,N
      JST=JST-M-1
      K=N+N+1-J
      PIV=1.D0/AUX(K)
      KST=K-N
      ID=IPIV(KST)-KST
      IST=2-J
      DO 25 K=1,L
      H=B(KST)
      IST=IST+N
      IEND=IST+J-2
      II=JST
      DO 24 I=IST,IEND
      II=II+M
   24 H=H-A(II)*X(I)
      I=IST-1
      II=I+ID
      X(I)=X(II)
      X(II)=PIV*H
   25 KST=KST+M
C
C
C     COMPUTATION OF LEAST SQUARES
   26 IST=N+1
      IEND=0
      DO 29 J=1,L
      IEND=IEND+M
      H=0.D0
      IF(M-N)29,29,27
   27 DO 28 I=IST,IEND
   28 H=H+B(I)*B(I)
      IST=IST+M
   29 AUX(J)=H
      RETURN
C
C     ERROR RETURN IN CASE M LESS THAN N
   30 IER=-2
      RETURN
C
C     ERROR RETURN IN CASE OF ZERO-MATRIX A
   31 IER=-1
      RETURN
C
C     ERROR RETURN IN CASE OF RANK OF MATRIX A LESS THAN N
   32 IER=K-1
      RETURN
      END
