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