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

     /* Program:                                            */
     /*                                                     */
     /*    PR - Computes inverse B of matrix A using        */
     /*         Pan-Reif algorithm.                         */
     /*                                                     */
     /* Version:                  Date:                     */
     /*                                                     */
     /*    Eco-C 1.1                 04/25/86               */
     /*                                                     */
     /* Description:                                        */
     /*                                                     */
     /*       lda   - Leading dimension of matrix  A.       */
     /*       n     - Order of matrix  A.                   */
     /*       Iterx - Maximum allowed iterations.           */
     /*       Residx - Convergence criterion.               */         
     /*                                                     */
     /* Subprograms:                                        */
     /*                                                     */
     /*    SDOT       from BLAS                             */
     /*    SSCAL      from BLAS                             */
     /*    SASUM      from BLAS                             */
     /*    SCOPY      from BLAS                             */
     /*                                                     */
     /* Reference:                                          */
     /*                                                     */
     /*    The Inversion of Large Matrices                  */
     /*    Thomas E. Phipps, Jr.                            */
     /*    BYTE Magazine                                    */
     /*    April, 1986                                      */
     /*                                                     */
     /* Authors:                                            */
     /*                                                     */
     /*    Adam Fritz                                       */
     /*    133 Main Street                                  */
     /*    Afton, New York 13730                            */

                                /* initialize inverse */
Initialize()
{
   int i ;
   float rs, cs ;               /* row and column sums */
   float rx, cx ;               /* maximum row and column sums */
   float SASUM() ;
   float t ;                    /* temporary variable */

                                /* compute scaling factor */
   rx = cx = 0.0 ;
   for ( i=0; i<n; i++ )
   {
      if ((rs = SASUM(n, &A[i][0], 1)) > rx) 
         rx = rs ;
      if ((cs = SASUM(n, &A[0][i], lda)) > cx)
         cx = cs ;
   }
   t = 1.0/(rx*cx) ;
Š                                /ª initializå inverså bù trans- */
                                /ª posinç anä scalinç originaì  */
                                /ª matriø                       */
   for ( i=0; i<n; i++ )
   {
      SCOPY(n, &A[i][0], 1, &B[0][i], lda) ;
      SSCAL(n, t, &B[0][i], lda) ;
   }
}

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

float MatrixResidual()
{
   int i, j ;
   float resid, t, SDOT(), Abs() ;

   resid = 0.0 ;
   for ( i=0; i<n; i++ )
   {
      for ( j=0; j<n; j++ )
      {
         t = SDOT(n, &B[i][0], 1, &A[0][j], lda) ;
         R[i][j] = -t ;
         t = Abs(t) ;
         if (i != j && t > resid)
            resid = t ;
      }
      t = 1.0 + R[i][i] ;
      R[i][i] = t ;
      t = Abs(t) ;
      if (t > resid)
         resid = t ;
   }
   return ( resid ) ;
}

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

UpdateInverse()
{
   int i, j ;
   float t, SDOT() ;
                                /* update residual matrix */
   for ( i=0; i<n; i++ )
   {
      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][0], 1) ;
      R[i][i] = 1.0 + R[i][i] ;
   }
                                /* compute new inverse */
   for ( i=0; i<n; i++ )
Š      for ( j=0; j<n; j++ )
         X[i][j] = SDOT(n, &R[i][0], 1, &B[0][j], lda) ;
                                /* update inverse */
   for ( i=0; i<n; i++ )
      SCOPY(n, &X[i][0], 1, &B[i][0], 1) ;
}

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

PR()
{
   int iter ;                   /* iteration counter */
   float MatrixResidual() ;       /* convergence flag */

                                /* initialize */
   Initialize() ;
                                /* compute residual matrix */
                                /* and update inverse      */
   iter = 0 ;
   while ((MatrixResidual() > Residx) && (iter++ < Iterx))
      UpdateInverse() ;
                                /* return iteration count */
   Iterx = iter ;
}

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