C.hr DSWEEP C@ C....*...1.........2.........3.........4.........5.........6.........7.* C DSWEEP 10/4/72 C C PURPOSE C INVERT A POSITIVE DEFINITE MATRIX IN PLACE BY A DIAGONAL SWEEP. C C USAGE C CALL DSWEEP(A,N,EPS,IER) C C ARGUMENTS C A - SYMMETRIC POSITIVE DEFINITE N BY N MATRIX STORED COLUMNWISE C (STORAGE MODE OF 0). ON RETURN CONTAINS THE INVERSE OF A C STORED COLUMNWISE. C REAL*8 C N - NUMBER OF ROWS AND COLUMNS OF A. C INTEGER C EPS - INPUT CONSTANT USED AS A RELATIVE TOLERANCE IN TESTING FOR C DEGENERATE RANK. A REASONABLE VALUE FOR EPS IS 1.D-13. C REAL*8 C IER - ERROR PARAMETER CODED AS FOLLOWS C IER=0 NO ERROR, RANK OF A IS N. C IER.GT.0 A IS SINGULAR, RANK OF A IS N-IER. C INTEGER C C REMARK C IF IER.GT.0 THEN DSWEEP RETURNS A G-INVERSE. C C REFERENCE C SCHATZOFF, M. ET AL. EFFICIENT CALCULATION OF ALL POSSIBLE REG- C RESSIONS. TECHNOMETRICS, 10. 769-79 (NOVEMBER 1968) C C PROGRAMMER C DR. A. RONALD GALLANT C DEPARTMANT OF STATISTICS C NORTH CAROLINA STATE UNIVERSITY C RALEIGH, NORTH CAROLINA 27695-8203 C C C SUBROUTINE DSWEEP(A,N,EPS,IER) IMPLICIT REAL*8 (A-H,O-Z) save REAL*8 A(N*N) TOL=0.D0 DO 5 I=1,N II=-N+N*I+I TEST=A(II) 5 IF(TEST.GT.TOL) TOL=TEST TOL=TOL*EPS IER=0 DO 50 K=1,N KK=-N+N*K+K AKK=A(KK) IF(AKK.GT.TOL) GO TO 20 DO 10 J=K,N KJ=-N+N*J+K 10 A(KJ)=0.D0 IF(K.EQ.1) GO TO 16 KLESS1=K-1 DO 15 I=1,KLESS1 IK=KK-I 15 A(IK)=0.D0 16 IER=IER+1 GO TO 50 20 D=1.D0/AKK DO 25 I=1,N DO 25 J=I,N IF((I.EQ.K).OR.(J.EQ.K)) GO TO 25 IJ=N*(J-1)+I IF(I.LT.K) AIK=A(-N+N*K+I) IF(I.GT.K) AIK=A(-N+N*I+K) IF(K.LT.J) AKJ=A(-N+N*J+K) IF(K.GT.J) AKJ=-A(-N+N*K+J) A(IJ)=A(IJ)-AIK*AKJ*D 25 CONTINUE DO 30 J=K,N KJ=-N+N*J+K 30 A(KJ)=A(KJ)*D IF(K.EQ.1) GO TO 36 KLESS1=K-1 DO 35 I=1,KLESS1 IK=KK-I 35 A(IK)=-A(IK)*D 36 A(KK)=D 50 CONTINUE DO 55 I=1,N DO 55 J=I,N IF(I.EQ.J) GO TO 55 IJ=-N+N*J+I JI=-N+N*I+J A(JI)=A(IJ) 55 CONTINUE RETURN END