File: HLS.FT of Disk: Disks/MyPDP/m8-2-rka1-rkb1
(Source file text)
SUBROUTINE HLS(A,B,M,N,L,IER,AUX,IPIV,EPS,X)
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 HLS (A,B,M,N,L,IER,AUX,IPIV,EPS,X)
C
C DESCRIPTION OF PARAMETERS
C A - M BY N COEFFICIENT MATRIX (DESTROYED).
C B - M BY L RIGHT HAND SIDE MATRIX (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 - 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 - INPUT PARAMETER WHICH SPECIFIES AN ABSOLUTE
C TOLERANCE FOR DETERMINATION OF RANK OF MATRIX A.
C IER - A RESULTING OUTPUT PARAMETER GIVING RANK OF
C MATRIX A OR ERROR CONDITION. (SEE REMARKS)
C AUX - AUXILIARY STORAGE ARRAY OF DIMENSION MAX(N,L)+N
C ON RETURN FIRST L LOCATIONS OF AUX CONTAIN THE
C RESULTING SUM OF LEAST SQUARES.
C
C REMARKS
C (1) NO ACTION BESIDES ERROR MESSAGE IER=-1002 IN CASE
C M LESS THAN N.
C (2) NO ACTION BESIDES ERROR MESSAGE IER=-1001 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 AND ALL PARAMETERS
C IMPROVED, IER IS SET TO N.
C (5) IF ALL THE ORIGINAL PARAMETERS ARE ACCEPTABLE,
C DEPENDING ON EPS, THE PRELIMINARY TERMINATION TEST IS
C SWITCHED OFF AND A FULL PASS MADE TO PERFORM THE
C COMPLETE HH TRANSFORMATION. EXIT WITH IER=-N.
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C HHSTEP, HHUK
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
DIMENSION A(1),B(1),X(1),IPIV(1),AUX(1)
C
C ERROR TEST
IF(M-N)30,1,1
C
C CALCN OF INITIAL VECTORS S(K),T(K) IN LOCS AUX(1),AUX(K1)
1 PIV=0.
K1 = MAX0(N,L)
IST = 1 - M
K2 = K1 + 1
DO 4 K=1,N
IPIV(K)=K
IST = IST + M
CALL HHPIV(A(IST),B,1,1,M,H,G)
AUX(K)=H
AUX(K2) = G
PIVT = G*G/H
IF(PIVT-PIV)4,4,3
3 PIV = PIVT
KPIV=K
4 K2 = K2 + 1
CALL HHSMSQ(B(1),1,M,F)
C
C ERROR TEST
IF(PIV)31,31,5
C
C DEFINE TOLERANCE FOR CHECKING RANK OF A, SET ZERO IF SOLN GOOD
5 TOL = EPS*EPS
IER = 1
C
C DECOMPOSITION LOOP
9 IST = - M
JB = 0
DO 21 K=1,N
IF(EPS.EQ.0.) GO TO 12
IF(EPS.GT.0.) GO TO 11
IF(PIV.GT.TOL) GO TO 12
GO TO 34
11 IF(F.GT.TOL*FLOAT(M-K+1)) GO TO 12
34 IF(K.NE.1) GO TO 32
TOL = 0.
IER = -1
12 F = F - PIV
IST = IST + M + 1
JB = JB + 1
LV = M-K+1
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
K2 = K1 + K
K3 = K1 + KPIV
G = AUX(K2)
AUX(K2) = AUX(K3)
AUX(K3) = G
JA = IST
JD = IST + I*M
NR = M - K + 1
CALL HHSWOP(A(IST),A(JD),1,NR)
C
C COMPUTATION OF PARAMETER SIG
C GENERATION OF VECTOR UK IN K-TH COLUMN OF MATRIX A AND OF
C PARAMETER BETA
8 CALL HHUK(A(IST),1,LV,SIG,BETA)
C
J = K1 + K
AUX(J)=-SIG
C
C SAVE INTERCHANGE INFORMATION
13 IPIV(KPIV)=IPIV(K)
IPIV(K)=KPIV
IF(K-N)14,19,19
C
C TRANSFORMATION OF MATRIX A
14 NK = N - K
IA = IST + M
CALL HHSTEP(A(IST),A(IA),1,M,LV,NK,BETA)
C
C TRANSFORMATION OF RIGHT HAND SIDE MATRIX B
19 IB = K
CALL HHSTEP(A(IST),B(IB),1,M,LV,L,BETA)
IF(K-N)10,21,21
C
C UPDATING OF S(K),T(K) ELEMENTS STORED IN AUX
10 PIV = 0.
KPIV = K + 1
J1 = KPIV
ID = IST
K2 = K1 + J1
DO 18 J=J1,N
ID = ID + M
H=AUX(J) - A(ID) * A(ID)
AUX(J)=H
G = AUX(K2) - A(ID) * B(JB)
AUX(K2) = G
PIVT = G*G/H
IF(PIVT-PIV)18,18,17
17 PIV = PIVT
KPIV=J
18 K2 = K2 + 1
21 CONTINUE
GO TO 20
C END OF DECOMPOSITION LOOP
C
C
C RANK OF MATRIX LESS THAN N, ZERO X,S AND BACK SUBSTITUTE
32 KR = K - 1
KT = KR
IER = KR
JK = KR - N + 1
JL = 0
DO 15 JY=1,L
JK = JK + N
JL = JL + N
DO 15 JX=JK,JL
15 X(JX) = 0.
GO TO 16
C
C BACK SUBSTITUTION AND BACK INTERCHANGE
20 IER = N * IER
KT = N
JK = 0
IK = N - M
K = K1 + N
PIV=1./AUX(K)
DO 22 K=1,L
JK = JK + N
IK = IK + M
22 X(JK) = PIV*B(IK)
KR = N - 1
C
16 IF(KR)26,26,23
23 JST = (KR+1)*(M+1)
IEND = M*(N-1) + KR + 1
DO 25 J=1,KR
JST = JST - M - 1
IEND = IEND - 1
KST = KR - J + 1
K = K1 + KST
PIV=1./AUX(K)
ID=IPIV(KST)-KST
KSTX = KST - N
KSTB = KST - M
DO 25 K=1,L
KSTX = KSTX + N
KSTB = KSTB + M
H = B(KSTB)
II = KSTX
DO 24 I=JST,IEND,M
II = II + 1
24 H = H - A(I) * X(II)
II = KSTX + ID
X(KSTX) = X(II)
X(II)=PIV*H
25 CONTINUE
C
C COMPUTATION OF LEAST SQUARES
26 IST = KT - M + 1
N1 = KT + 1
H = 0.
DO 29 J=1,L
IST = IST + M
H=0.
JA = IST
IF(M-KT)29,29,27
27 NR = M - N1 + 1
IF(NR.EQ.1) GO TO 28
CALL HHSMSQ(B(IST),1,NR,H)
GO TO 29
28 H = B(IST)*B(IST)
C
29 AUX(J)=H
RETURN
C
C ERROR RETURN IN CASE M LESS THAN N
30 IER = -1002
RETURN
C
C ERROR RETURN IN CASE OF ZERO-MATRIX A
31 IER = -1001
RETURN
END