Directory of image this file is from
This file as a plain text file

      SUBROUTINE LLSQ(A,B,M,N,L,X,IPIV,EPS,IER,AUX)
      DIMENSION A(1),B(1),X(1),IPIV(1),AUX(1)
      IF(M-N)30,1,1
    1 PIV=0.
      IEND=0
      DO 4 K=1,N
      IPIV(K)=K
      H=0.
      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
      IF(PIV)31,31,5
    5 SIG=SQRT(PIV)
      TOL=SIG*ABS(EPS)
      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
    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
    8 IF(K-1)11,11,9
    9 SIG=0.
      DO 10 I=IST,IEND
   10 SIG=SIG+A(I)*A(I)
      SIG=SQRT(SIG)
      IF(SIG-TOL)32,32,11
   11 H=A(IST)
      IF(H)12,13,13
   12 SIG=-SIG
   13 IPIV(KPIV)=IPIV(K)
      IPIV(K)=KPIV
      BETA=H+SIG
      A(IST)=BETA
      BETA=1./(SIG*BETA)
      J=N+K
      AUX(J)=-SIG
      IF(K-N)14,19,19
   14 PIV=0.
      ID=0
      JST=K+1
      KPIV=JST
      DO 18 J=JST,N
      ID=ID+M
      H=0.
      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
      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
   19 DO 21 J=K,LM,M
      H=0.
      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
      IER=0
      I=N
      LN=L*N
      PIV=1./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./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
   26 IST=N+1
      IEND=0
      DO 29 J=1,L
      IEND=IEND+M
      H=0.
      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
   30 IER=-2
      RETURN
   31 IER=-1
      RETURN
   32 IER=K-1
      RETURN
      END