home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / eigenlib / elmhes.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-11-18  |  2.5 KB  |  99 lines
  1. elmhes(nm,n,low,igh,a,intch)
  2.  
  3. int n,nm,*igh,*low;
  4. double **a;
  5. int *intch;
  6. {
  7. int i,j,m,la,kp1,mm1,mp1;
  8. double x,y,fabs();
  9.  
  10. /*    this subroutine is a translation of the algol procedure elmhes,
  11.       num. math. 12, 349-368(1968) by martin and wilkinson.
  12.       handbook for auto. comp., vol.ii-linear algebra, 339-358(1971).
  13.  
  14.       given a real general matrix, this subroutine
  15.       reduces a submatrix situated in rows and columns
  16.       low through igh to upper hessenberg form by
  17.       stabilized elementary similarity transformations.
  18.  
  19.       on input
  20.  
  21.          nm must be set to the row dimension of two-dimensional
  22.            array parameters as declared in the calling program
  23.            dimension statement.
  24.  
  25.          n is the order of the matrix.
  26.  
  27.          low and igh are integers determined by the balancing
  28.            subroutine  balanc.  if  balanc  has not been used,
  29.            set low=1, igh=n.
  30.  
  31.          a contains the input matrix.
  32.  
  33.       on output
  34.  
  35.          a contains the hessenberg matrix.  the multipliers
  36.            which were used in the reduction are stored in the
  37.            remaining triangle under the hessenberg matrix.
  38.  
  39.          int contains information on the rows and columns
  40.            interchanged in the reduction.
  41.            only elements low through igh are used.
  42.  
  43.       this routine is a C-translation of the FORTRAN 77 source code
  44.       written by the mathematics and computer science division,
  45.       argonne national laboratory
  46.       last change :   september 1989.
  47.  
  48.       mark myers
  49.       Center for Applied Mathematics 
  50.       Cornell University    (607) 255-4195
  51.  
  52.       ----------------------------------------------------------    */
  53.  
  54.  la = *igh - 1;
  55.  kp1 = *low + 1;
  56.  if (la < kp1)
  57.    return;
  58.  for (m=kp1;m<=la;m++) 
  59.     {mm1 = m - 1;
  60.      x = 0.e0;
  61.      i = m;
  62.  
  63.      for (j=m;j<=*igh;j++) 
  64.        if (fabs(a[j][mm1]) > fabs(x)) 
  65.           {x = a[j][mm1];
  66.            i = j;}
  67.  
  68.      intch[m] = i;
  69.      if (i !=  m) 
  70.        {for (j=mm1;j<=n;j++)    /* interchange rows & col of a */
  71.           {y = a[i][j];
  72.            a[i][j] = a[m][j];
  73.            a[m][j] = y; }
  74.  
  75.         for (j=1;j<=*igh;j++) 
  76.           {y = a[j][i];
  77.            a[j][i] = a[j][m];
  78.            a[j][m] = y;}
  79.        }
  80.  
  81.      if (x == 0.e0) 
  82.        continue;
  83.      mp1 = m + 1;
  84.  
  85.      for (i=mp1;i<=*igh;i++) 
  86.        {y = a[i][mm1];
  87.         if (y == 0.e0)
  88.       continue;
  89.         y = y / x;
  90.         a[i][mm1] = y;
  91.  
  92.         for (j=m;j<=n;j++)
  93.           a[i][j] = a[i][j] - y * a[m][j]; 
  94.  
  95.         for (j=1;j<=*igh;j++)
  96.           a[j][m] = a[j][m] + y * a[j][i];
  97.      } }     
  98. }
  99.