{ Copyright (C) 1986 Adam Fritz, 133 Main St., Afton, N.Y. 13730 }

procedure PanReif ( var a,b,r,x : MATRIX ;
                        lda,n : integer ;
                    var iterx : integer ;
                        residx : real ) ;
     {                                                     }
     { PROGRAM:                                            }
     {                                                     }
     {    Single precision, real matrix inversion.         }
     {    Pan-Reif method, Reference 1.                    }
     {                                                     }
     { VERSION:                  DATE:                     }
     {                                                     }
     {    TURBO Pascal 1.1          April 25, 1986         }
     {                                                     }
     { DESCRIPTION:                                        }
     {                                                     }
     {    Initialize inverse                               }
     {    Iteration                                        }
     {       Compute residual matrix                       }
     {       Update inverse                                }
     {                                                     }
     { REFERENCE:                                          }
     {                                                     }
     {    The Inversion of Large Matrices                  }
     {    Thomas E. Phipps, Jr.                            }
     {    BYTE Magazine                                    }
     {    April, 1986                                      }
     {                                                     }
     { AUTHORS:                                            }
     {                                                     }
     {    Adam Fritz (TURBO Pascal)                        }
     {                                                     }

var
     i,j            : integer ;

     rs,cs          : real ;    { row and column sums }
     rx,cx          : real ;    { maximum row and column sums }
     t              : real ;    { temporary variable }
     iter           : integer ; { iteration count }
     resid          : real ;    { residual }

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

procedure Initialize ;

begin
                                { initialize inverse }
   rx := 0.0 ;
   cx := 0.0 ;
   for i := 1 to n do begin
      rs := SASUM(n,a[i,1],1) ;
      if rs > rx then
         rx := rs ;
      cs := SASUM(n,a[1,i],lda) ;
      if cs > cx then
         cx := cs
   end ;
   t := 1.0/(rx*cx) ;

   for i := 1 to n do begin
      SCOPY(n,a[i,1],1,b[1,i],lda) ;
      SSCAL(n,t,b[1,i],lda) ;
   end

end ;

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

function MatrixResidual : real ;

begin
   rx := 0.0 ;
   for i := 1 to n do begin
      for j := 1 to n do begin
         t := SDOT(n,b[i,1],1,a[1,j],lda) ;
         r[i,j] := -t ;
         t := abs(t) ;
         if i <> j then
            if t > rx then
               rx := t
      end ;
      t := 1.0 + r[i,i] ;
      r[i,i] := t ;
      t := abs(t) ;
      if t > rx then
         rx := t
   end ;
   MatrixResidual := rx
end ;

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

procedure UpdateInverse ;

begin
                                { update resid matrix }
   for i := 1 to n do begin
      t := r[i,i] ;
{      t := 1.0+(1.0+(1.0+(1.0+t)*t)*t)*t ; }
      t := 1.0 + t ;
      SSCAL(n,t,r[i,1],1) ;
      r[i,i] := 1.0 + r[i,i]
   end ;
                                { compute new inverse }
   for i := 1 to n do
      for j := 1 to n do
         x[i,j] := SDOT(n,r[i,1],1,b[1,j],lda) ;
                                { update inverse }
   for i := 1 to n do
      SCOPY(n,x[i,1],1,b[i,1],1)

end ;

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

begin
                                { initialize }
   Initialize ;
                                { iteration with computation of }
                                { residual matrix and inverse   }
                                { update                        }
   iter := 0 ;
   while (MatrixResidual > residx) and (iter < iterx) do begin
      iter := iter + 1 ;
      UpdateInverse
   end ;
                                { return iteration count }
   iterx := iter

end ;

{ Copyright (C) 1986 Adam Fritz, 133 Main St., Afton, N.Y. 13730 }