この関数は，最初は FORTRAN で書いたものを C に移して，さらに AWK に移植したのですが，AWK の変な制限のためにひねくれたコードになってしまった箇所があります。

　この実対称行列の固有値・固有ベクトルを求めるアルゴリズムはかなり高速なので，AWK を使っても見劣りしないかもしれません。10 変数 47 ケースのデータを使って，EPSON PC-486 NOTE AS で，１０秒で結果が出ました。

　固有値・固有ベクトルを求める関数があれば，数量化II類や正準判別分析も書けますね。あと，逆行列を求める関数があれば，鬼に金棒か(^_^;)

# 実対称行列の，固有値(e)，固有ベクトル(v)を求める関数 eigen.awk

function eigen2(a, n, ne, eps, e, v,     sw, i, j, k, k1, m, nea, nm1,
nm2, nv, nva, a1, ajk, d, ei, eps1, eps2, f, q, r, s, sr, t, wj1, w, lw)
{
   nv = nva = nea = abs(ne)
   if (nea == 0 || nea > n) return 1
   nm1 = n-1; nm2 = n-2
   if (eps <= 0) {
       eps = 1e-16
   }
   if (n == 1) {
       e[0] = a[0,0]
       if (nv != 0) {
           v[0,0] = 1
       }
       return 0
   }
   if (n == 2) {
       w[0,0] = a[0,1]; t = 0.5*(a[0,0]+a[1,1])
       r = a[0,0]*a[1,1]-a[0,1]^2
       d = t*t-r
       q = abs(t)+sqrt(d)
       if (t < 0) {
           q = -q
       }
       t *= ne
       if (t >= 0) {
           e[0] = q
           if (nea == 2) {
               e[1] = r/q
           }
       }
       else {
           e[0] = r/q
           if (nea == 2) {
               e[1] = q
           }
       }
   }
   else {
       for (k = 1; k <= nm2; k++) {
           k1 = k-1
           s = 0
           for (i = k; i < n; i++) {
               s += a[k1,i]^2
           }
           w[0,k1] = 0
           if (s != 0) {
               sr = sqrt(s)
               a1 = a[k1,k]
               if (a1 < 0) {
                   sr = -sr
               }
               w[0,k1] = -sr
               r = 1/(s+a1*sr)
               a[k1,k] = a1+sr
               for (i = k; i < n; i++) {
                   s = 0
                   for (j = k; j <= i; j++) {
                       s += a[j,i]*a[k1,j]
                   }
                   if (i != nm1) {
                       for (j = i+1; j < n; j++) {
                           s += a[i,j]*a[k1,j]
                       }
                   }
                   w[0,i] = s*r
               }
               s = 0
               for (i = k; i < n; i++) {
                   s += a[k1,i]*w[0,i]
               }
               t = 0.5*s*r
               for (i = k; i < n; i++) {
                   w[0,i] -= t*a[k1,i]
               }
               for (j = k; j < n; j++) {
                   wj1 = w[0,j]
                   ajk = a[k1,j]
                   for (i = j; i < n; i++) {
                       a[j,i] -= a[k1,i]*wj1+w[0,i]*ajk
                   }
               }
           }
       }
       w[0,nm2] = a[nm2,nm1]
       for (i = 0; i < n; i++) {
           w[5,i] = a[i,i]
       }
       r = abs(w[5,0])+abs(w[0,0])
       f = abs(w[0,nm2])+abs(w[5,nm1])
       r = max(r, f)
       for (i = 1; i < nm1; i++) {
           f = abs(w[0,i-1])+abs(w[5,i])+abs(w[0,i])
           r = max(r, f)
       }
       eps1 = r*1e-16
       eps2 = r*eps
       for (i = 0; i < nm1; i++) {
           w[1,i] = w[0,i]^2
       }
       if (ne < 0) {
           r = -r
       }
       f = r
       for (i = 0; i < nea; i++) {
           e[i] = -r
       }
       pass2 = 0
       for (k = 1; k <= nea; k++) {
           if (!pass2) {
               k1 = k-1
               d = e[k1]
           }
           pass2 = 0
           t = 0.5*(d+f)
           if (abs(d-f) > eps2 && t != d && t != f) {
               j = 0
               i = 1
               for (;;) {
                   if (!pass1) {
                       q = w[5,i-1]-t
                   }
                   pass1 = 0
                   if (q >= 0) {
                       j++
                   }
                   if (q != 0) {
                       if (++i > n) {
                           break
                       }
                       q = w[5,i-1]-t-w[1,i-2]/q
                       pass1 = 1
                       continue
                   }
                   else {
                       i += 2
                   }
                   if (i > n) {
                       break
                   }
               }
               if (ne < 0) {
                   j = n-j
               }
               if (j < k) {
                   f = t
                   pass2 = 1
                   k--
                   continue
               }
               d = t
               m = min(j, nea)
               for (i = k-1; i < m; i++) {
                   e[i] = t
               }
               pass2 = 1
               k--
               continue
           }
           e[k1] = t
       }
   }
   if (nv != 0) {
       if (n == 2) {
           w[5,0] = a[0,0]
           w[5,1] = a[1,1]
       }
       w[0,nm1] = 0
       for (i = 0; i < nva; i++) {
           ei = e[i]
           for (j = 0; j < n; j++) {
               w[1,j] = w[5,j]-ei
               w[2,j] = w[0,j]
               v[i,j] = 1
           }
           sw = 0
           for (j = 0; j < nm1; j++) {
               if (abs(w[1,j]) >= abs(w[0,j])) {
                   if (w[1,j] == 0) {
                       w[1,j] = 1e-30
                   }
                   w[4,j] = w[0,j]/w[1,j]
                   lw[j] = 0
                   w[1,j+1] -= w[4,j]*w[2,j]
                   w[3,j] = 0
               }
               else {
                   w[4,j] = w[1,j]/w[0,j]
                   lw[j] = 1
                   w[1,j] = w[0,j]
                   t = w[2,j]
                   w[2,j] = w[1,j+1]
                   w[3,j] = w[2,j+1]
                   w[1,j+1] = t-w[4,j]*w[2,j]
                   w[2,j+1] = -w[4,j]*w[3,j]
               }
           }
           if (w[1,nm1] == 0) {
               w[1,nm1] = 1e-30
           }
           if (i > 0 && abs(e[i]-e[i-1]) < eps1) {
               for (j = 0; j < n; j++) {
                   v[i,j] = rand()
               }
           }
           for (;;) {
               t = v[i,nm1]
               r = v[i,nm2]
               v[i,nm1] = t/w[1,nm1]
               v[i,nm2] = (r-w[2,nm2]*v[i,nm1])/w[1,nm2]
               for (k = nm2-1; k >= 0; k--) {
                   v[i,k] = (v[i,k]-w[2,k]*v[i,k+1]-w[3,k]*v[i,k+2])/w[1,k]
               }
               if (sw == 0) {
                   sw = 1
                   for (j = 0; j < nm1; j++) {
                       if (lw[j] == 1) {
                           t = v[i,j]
                           v[i,j] = v[i,j+1]
                           v[i,j+1] = t-w[4,j]*v[i,j+1]
                       }
                       else {
                           v[i,j+1] -= w[4,j]*v[i,j]
                       }
                   }
                   continue
               }
               else
               break
           }
       }
       if (n > 2) {
           for (i = 0; i < nm2; i++) {
               w[0,i] *= -a[i,i+1]
           }
           for (i = 0; i < nva; i++) {
               for (k = nm2-1; k >= 0; k--) {
                   r = w[0,k]
                   if (r != 0) {
                       r = 1/r
                       s = 0
                       for (j = k+1; j < n; j++) {
                           s += a[k,j]*v[i,j]
                       }
                       r *= s
                       for (j = k+1; j < n; j++) {
                           v[i,j] -= r*a[k,j]
                       }
                   }
               }
           }
       }
       for (i = 0; i < nva; i++) {
           if (i == 0) {
               m = 1
           }
           else {
               if (abs(e[i]-e[i-1]) >= eps1) {
                   m = i+1
               }
               else {
                   for (j = m-1; j < i; j++) {
                       s = 0
                       for (k = 0; k < n; k++) {
                           s += v[j,k]*v[i,k]
                       }
                       for (k = 0; k < n; k++) {
                           v[i,k] -= s*v[j,k]
                       }
                   }
               }
           }
           t = s = 0
           for (j = 0; j < n; j++) {
               s += v[i,j]^2
           }
           if (s != 0) {
               t = sqrt(1/s)
           }
           for (j = 0; j < n; j++) {
               v[i,j] *= t
           }
       }
   }
   return 0
}