COPYRIGHT (C) 1986 ADAM FRITZ, 133 MAIN ST., AFTON, N.Y. 13730
C
      SUBROUTINE PR ( A, B, R, X, LDA, N, ITERX, RESIDX )
C
      REAL A(LDA,1), B(LDA,1), R(LDA,1), X(LDA,1)
      INTEGER I, J, LDA, N, ITER, ITERX
      REAL RESID, RESIDX, RS, CS, RX, CX, T
      LOGICAL DONE
C                                                      
C PROGRAM:                                             
C                                                      
C    SINGLE PRECISION, REAL MATRIX INVERSION.          
C    PAN-REIF METHOD, REFERENCE 1.                     
C                                                      
C VERSION:                  DATE:                      
C                                                      
C    MICROSOFT FORTRAN         APRIL 25, 1986          
C                                                      
C DESCRIPTION:                                         
C                                                      
C    INITIALIZE INVERSE                                
C    ITERATION                                         
C       COMPUTE RESIDUAL MATRIX                        
C       UPDATE INVERSE                                 
C                                                      
C REFERENCE:                                           
C                                                      
C    THE INVERSION OF LARGE MATRICES                   
C    THOMAS E. PHIPPS, JR.                             
C    BYTE MAGAZINE                                     
C    APRIL, 1986                                       
C                                                      
C AUTHORS:                                             
C                                                      
C    ADAM FRITZ 
C                                                      
C     INITIALIZE
C
      ITER = 0
C
C     INITIALIZE INVERSE  
C
      RX = 0.0  
      CX = 0.0  
      DO 111 I = 1, N  
         RS = SASUM(N,A(I,1),LDA)  
         IF (RS .GT. RX)
     *      RX = RS  
         CS = SASUM(N,A(1,I),1)  
         IF (CS .GT. CX)
     *      CX = CS
  111 CONTINUE
      T = 1.0/(RX*CX)  

      DO 113 I = 1, N  
         CALL SCOPY(N,A(I,1),LDA,B(1,I),1)  
         CALL SSCAL(N,T,B(1,I),1)  
  113 CONTINUE
C
C     ITERATION  
C
  115 CONTINUE
C
C        COMPUTE RESIDUAL MATRIX  
C
         RESID = 0.0
         DO 119 I = 1, N  
            DO 117 J = 1, N  
               T = SDOT(N,B(I,1),LDA,A(1,J),1)  
               R(I,J) = -T  
               T = ABS(T)
               IF (I .NE. J .AND. T .GT. RESID) RESID = T
  117       CONTINUE  
            T = 1.0 + R(I,I)  
            R(I,I) = T
            T = ABS(T)  
            IF (T .GT. RESID) RESID = T
  119    CONTINUE  
C
C        DONE?
C
         IF (RESID .LE. RESIDX .OR. ITER .GE. ITERX) GO TO 127
C
            ITER = ITER + 1  
C
C           UPDATE RESIDUAL MATRIX  
C  
            DO 121 I = 1, N  
               T = R(I,I)  
               T = 1.0 + T  
C              T = 1.0+(1.0+(1.0+(1.0+T)*T)*T)*T    
               CALL SSCAL(N,T,R(I,1),LDA)  
               R(I,I) = 1.0 + R(I,I)
  121       CONTINUE  
C
C           COMPUTE NEW INVERSE  
C
            DO 123 I = 1, N
               DO 123 J = 1, N
                  X(I,J) = SDOT(N,R(I,1),LDA,B(1,J),1)  
  123       CONTINUE
C
C           UPDATE INVERSE  
C
            DO 125 I = 1, N
               CALL SCOPY(N,X(I,1),LDA,B(I,1),LDA)
  125       CONTINUE
C
C        LOOP
C
         GO TO 115
C
C        EXIT
C
  127 CONTINUE
C
      ITERX = ITER
      RESIDX = RESID
C
      RETURN
      END
C
COPYRIGHT (C) 1986 ADAM FRITZ, 133 MAIN ST., AFTON, N.Y. 13730