SUBROUTINE DECOMP(N,K,ETA,TOL,IRANK,INSING) C PURPOSE--THIS SUBROUTINE DECOMPOSES THE WEIGHTED DATA C MATRIX Q WHICH ORIGINALLY = THE N BY K DATA MATRIX X C TIMES THE SQUARE ROOT OF THE WEIGHTS (IN W). C THE ORIGINAL Q IS DECOMPOSED INTO A NEW Q TIMES THE C INVERSE OF A DIAGONAL MATRIX D TIMES THE DIAGONAL MATRIX D C TIMES AN UPPER TRIANGULAR MATRIX R. C THE NEW N BY K Q HAS ORTHOGONAL COLUMNS. C A SECOND OUTPUT FROM THIS SUBROUTINE IS THE RANK AND C STATUS (NON-SINGULAR OR SINGULAR) OF THE DATA MATRIX X. C A THIRD OURPUT FROM THIS SUBROUTINE IS THE NUMERICALLY C OPTIMAL PIVOT POINTS FOR THE DECOMPOSITION. C X--NOT USED C Q--USED AND CHANGED C R--DEFINED C D--PERMANENTLY DEFINED C IPIVOT--PERMANENTLY DEFINED C UPDATED --NOVEMBER 1975. C UPDATED --FEBRUARY 1976. C C--------------------------------------------------------------------- C LOGICAL FSUM DIMENSION Q(10000),R(2500),D(50),IPIVOT(50) COMMON /BLOCK2/ WS(15000) COMMON /BLOCK3/ DUM1(3000),DUM2(3000) EQUIVALENCE (Q(1),WS(1)) EQUIVALENCE (R(1),WS(10001)) EQUIVALENCE (D(1),WS(12501)) EQUIVALENCE (IPIVOT(1),WS(12551)) C C-----START POINT----------------------------------------------------- C C ZERO-OUT SOME VARIABLES, VECTORS, AND ARRAYS C INSING=0 IRANK=0 DO 5 J=1,K D(J) = 0.0 DO 6 I=1,K IRARG=(I-1)*K+J R(IRARG)=0.0 6 CONTINUE 5 CONTINUE C TOL2=TOL*TOL DO 10 J=1,K 10 IPIVOT(J)=J DO 200 IS=1,K C C BEGIN STEP NUMBER IS IN THE DECOMPOSITION C IF (IS.EQ.1) GO TO 20 GO TO 30 20 FSUM=.TRUE. 30 DIS=0.0 IP=IS C C BEGIN THE PIVOT SEARCH C DO 80 J=IS,K M=IPIVOT(J) IF (FSUM) GO TO 40 GO TO 60 40 DO 50 L=1,N IQARG=(L-1)*K+M DUM1(L)=Q(IQARG) 50 DUM2(L)=Q(IQARG) C CALL DOT(DUM1,DUM2,1,N,0.0,D(J)) C 60 IF (DIS.LT.D(J)) GO TO 70 GO TO 80 70 DIS=D(J) IP=J 80 CONTINUE C C END THE PIVOT SEARCH C M=IPIVOT(IP) IF (FSUM) DN=DIS IF (DIS.LT.ETA*DN) GO TO 90 FSUM=.FALSE. GO TO 100 90 FSUM=.TRUE. 100 IF (FSUM) GO TO 30 IF (IP.NE.IS) GO TO 110 GO TO 130 C C BEGIN COLUMN INTERCHANGES C 110 D(IP)=D(IS) IPIVOT(IP)=IPIVOT(IS) IPIVOT(IS)=M IF (IS.EQ.1) GO TO 130 ISM1=IS-1 DO 120 I=1,ISM1 IRARG1=(I-1)*K+IP IRARG2=(I-1)*K+IS HOLD=R(IRARG1) R(IRARG1)=R(IRARG2) 120 R(IRARG2)=HOLD C C END COLUMN INTERCHANGES C 130 DO 140 L=1,N IQARG=(L-1)*K+M DUM1(L)=Q(IQARG) 140 DUM2(L)=Q(IQARG) C CALL DOT(DUM1,DUM2,1,N,0.0,D(IS)) C DIS=D(IS) IF (DIS.LE.TOL2*D(1)) RETURN IF(DIS.NE.0.0)GOTO150 INSING=0 RETURN 150 ISP1=IS+1 IF (ISP1.GT.K) GO TO 190 C C BEGIN ORTHOGONALIZATION C DO 180 J=ISP1,K IP=IPIVOT(J) DO 160 L=1,N IQARG1=(L-1)*K+M IQARG2=(L-1)*K+IP DUM1(L)=Q(IQARG1) 160 DUM2(L)=Q(IQARG2) C IRARG=(IS-1)*K+J CALL DOT(DUM1,DUM2,1,N,0.0,R(IRARG)) R(IRARG)=R(IRARG)/DIS C RISJ=R(IRARG) DO 170 I=1,N IQARG1=(I-1)*K+IP IQARG2=(I-1)*K+M 170 Q(IQARG1)=Q(IQARG1)-RISJ*Q(IQARG2) 180 D(J)=D(J)-DIS*RISJ*RISJ C C END ORTHOGONALIZATION C 190 IRANK=IS 200 CONTINUE C C END STEP NUMBER IS INTHE DECOMPOSITION C INSING=1 RETURN END