#/***********************************************************
#    eigen.c -- QR法
#***********************************************************/
require "matutil.rb"      # 行列用小道具集 

EPS = 1E-6               # 非対角成分の許容誤差
MAX_ITER = 100           # 最大の繰返し数

def house( n, offset, x )  # Householder変換 
    y = x[offset,n]
    s = Math::sqrt(innerproduct(n, y, y));  # 内積の平方根, すなわち大きさ 
    if (s != 0) 
        if (x[offset] < 0); s = -s; end
        x[offset] += s;  t = 1 / Math::sqrt(x[offset] * s)
        for i in 0...n; x[i + offset] *= t; end
    end
    return -s
end

def tridiagonalize( n, a, d, ee) # 3重対角化 
    e = ee.dup; e.shift
    for k in 0...(n - 2)
        v = a[k];  d[k] = v[k]
        e[k] = house(n - k - 1, k + 1, v)
        if (e[k] == 0); next; end
        for i in (k + 1)...n
            s = 0
            for j in (k + 1)...i; s += a[j][i] * v[j]; end
            for j in i...n; s += a[i][j] * v[j]; end
            d[i] = s
        end
        t = innerproduct(n-k-1, v[k+1,n-k-1], d[k+1,n-k-1]) / 2
        (n - 1).downto(k + 1) do |i|
            p = v[i];  q = d[i] - t * p;  d[i] = q
            for j in i...n
                a[i][j] -= p * d[j] + q * v[j]
            end
        end
    end
    if (n >= 2);   d[n - 2] = a[n - 2][n - 2]
                   e[n - 2] = a[n - 2][n - 1];  end
    if (n >= 1);    d[n - 1] = a[n - 1][n - 1]; end
    (n - 1).downto(0) do |i| 
        v = a[k]
        if (k < n - 2) 
            for i in (k + 1)...n
                w = a[i] 
                t = innerproduct(n-k-1, v[k+1,n-k-1], w[k+1,n-k-1])
                for j in (k + 1)...n
                    w[j] -= t * v[j]
                end
            end
        end
        ee = [ee[0]] + e
        for i in 0...n; v[i] = 0; end
        v[k] = 1
    end
    return ee
end

def eigen( n, a, d, e )
    e = tridiagonalize(n, a, d, e)  # 3重対角化 
    e[0] = 0  # 番人 
    (n - 1).downto(1) do |h|
        j = h
        while (e[j].abs > EPS * (d[j - 1].abs + d[j].abs))
            j -= 1;  # $\mbox\tt e[$j$]end \ne 0$ のブロックの始点を見つける 
        end
        if (j == h); next; end
        iter = 0
        begin 
            if ((iter += 1) > MAX_ITER); return 1; end
            w = (d[h - 1] - d[h]) / 2
            t = e[h] * e[h]
            s = Math::sqrt(w * w + t);  if (w < 0); s = -s; end
            x = d[j] - d[h] + t / (w + s);  y = e[j + 1]
            for k in j...h
                if (x.abs >= y.abs) 
                    t = -y / x;  c = 1 / Math::sqrt(t * t + 1)
                    s = t * c
                else 
                    t = -x / y;  s = 1 / Math::sqrt(t * t + 1)
                    c = t * s
                end
                w = d[k] - d[k + 1]
                t = (w * s + 2 * c * e[k + 1]) * s
                d[k] -= t;  d[k + 1] += t
                if (k > j); e[k] = c * e[k] - s * y; end
                e[k + 1] += s * (c * w - 2 * s * e[k + 1])
                # 次の5行は固有ベクトルを求めないなら不要 
                for i in 0...n
                    x = a[k][i];  y = a[k + 1][i]
                    a[k    ][i] = c * x - s * y
                    a[k + 1][i] = s * x + c * y
                end
                if (k < h - 1) 
                    x = e[k + 1];  y = -s * e[k + 2]
                    e[k + 2] *= c
                end
            end
        end while (e[h].abs > \
                EPS * (d[h - 1].abs + d[h].abs))
    end
#/*
#  以下は固有値の降順に整列しているだけ. 必要なければ省く.
#  固有ベクトルを求めないなら固有ベクトルの整列はもちろん不要.
#  なお, \tt matutil.cend 中で行列の各行をベクトルへのポインタ
#  として定義しているので, 行交換はポインタのすげ替えだけで済む.
#*/
    for k in 0...n
        h = k;  t = d[h]
        for i in k...n
            if (d[i] > t);  h = i;  t = d[h];  end
        end
        d[h] = d[k];  d[k] = t
        v = a[h];  a[h] = a[k];  a[k] = v
    end
    return 0
end

#************ 以下はテスト用 ****************

ULONG_MAX = 4294967295
ULONG_MAXplus = 4294967296

def rnd  # 乱数  0 < rnd() < 1 
    $seed = $seed * 69069 % ULONG_MAXplus
    return $seed / (ULONG_MAX + 1.0)
end

printf("n = ");  n = gets.to_i
printf("乱数の種 (正の整数) = ")
$seed = gets.to_i;  seed |= 1
a = new_matrix(n, n);  u = new_matrix(n, n)
d = new_vector(n);  work = new_vector(n)
for i in 0...n
    for j in i...n
        u[i][j] = a[i][j] = a[j][i] = rnd() - rnd()
    end
end
printf("A:\n")
matprint(a, n, 7, "%10.6f")
if (eigen(n, u, d, work) == 1)
    printf("収束しません.\n")
end
printf("固有値:\n")
vecprint(d, n, 5, "% -14g")
s = 0
for i in 0...n
    for j in 0...n
        t = 0
        for k in 0...n
            t += u[k][i] * d[k] * u[k][j]
        end
        s += (a[i][j] - t) * (a[i][j] - t)
    end
end
printf("二乗平均誤差: %g\n", Math::sqrt(s) / n)

#exit 0