SUBROUTINE INVXWX(N,K) C PURPOSE--THIS SUBROUTINE COMPUTES THE INVERSE OF X'WX C WHICH IS DONE BY COMPUTING THE INVERSE OF R'R (WHERE C R HAS JUST RECENTLY BEEN MODIFIED BEFORE CALLING THIS C SUBROUTINE. THE INPUT R = THE SQUARE ROOT OF C THE DIAGONAL MATRIX D TIMES THE OLD MATRIX R. C THE INVERSE OF X'WX WILL BE IDENTICAL C (EXCEPT FOR THE ABSENCE OF S**2 = THE RESIDUAL C VARIANCE) TO THE COVARIANCE MATRIX OF THE COEFFICIENTS. C THE ONLY REASON THIS SUBROUTINE EXISTS IS FOR THE C CALCULATION OF SUCH COVARIANCES. C UNPIVOTING HAS ALSO BEEN DONE HEREIN SO AS TO UNDO C THE PIVOTING DONE IN THE DECOMPOSITION SUBROUTINE (DECOMP). C THE MATRIX C USED HEREIN IS AN INTERMEDIATE RESULT MATRIX. C X--NOT USED C Q--NOT USED C R--USED AND CHANGED C D--NOT USED C IPIVOT--USED C INVERSION ALGORITHM USED--CHOLESKI DECOMPOSITION C UPDATED --NOVEMBER 1975. C UPDATED --FEBRUARY 1976. C C--------------------------------------------------------------------- C 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)) DIMENSION DUM3(200) C C-----START POINT----------------------------------------------------- C DO 10 I=1,K IM1=I-1 IF(IM1.LT.1)GOTO10 DO15J=1,IM1 IRARG=(I-1)*K+J R(IRARG)=0.0 15 CONTINUE 10 CONTINUE DO30JJ=1,K J=K+1-JJ DO 30 II=1,J I=J+1-II IP1=I+1 IF(IP1.GT.K)GOTO25 DO20L=IP1,K IRARG1=(I-1)*K+L IRARG2=(J-1)*K+L IRARG3=(L-1)*K+J DUM1(L)=R(IRARG1) IF(L.LT.J)DUM2(L)=R(IRARG2) IF(L.EQ.J)DUM2(L)=DUM3(L) IF(L.GT.J)DUM2(L)=R(IRARG3) 20 CONTINUE 25 RI=0.0 IRARG=(I-1)*K+I IF (I.EQ.J) RI=1.0/R(IRARG) ANEGRI=-RI C CALL DOT(DUM1,DUM2,IP1,K,ANEGRI,DOTPRO) C IRARG=(I-1)*K+I DOTPRO=-DOTPRO/R(IRARG) IF(I.EQ.J)DUM3(I)=DOTPRO IRARG=(J-1)*K+I IF(I.LT.J)R(IRARG)=DOTPRO 30 CONTINUE DO35I=1,K IRARG=(I-1)*K+I R(IRARG)=DUM3(I) 35 CONTINUE C C MATRIX C NOW EQUALS THE INVERSE OF R'R. C NOW 'UNPIVOT' ON C AND PUT THE RESULTS BACK INTO R. C DO40I=1,K II=IPIVOT(I) DO40J=1,I JJ=IPIVOT(J) IRARG1=(II-1)*K+JJ IRARG2=(I-1)*K+J IRARG3=(JJ-1)*K+II IF(II.LT.JJ)R(IRARG1)=R(IRARG2) IF(II.EQ.JJ)DUM3(II)=R(IRARG2) IF(II.GT.JJ)R(IRARG3)=R(IRARG2) 40 CONTINUE DO50I=1,K IRARG=(I-1)*K+I R(IRARG)=DUM3(I) 50 CONTINUE RETURN END