using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;

namespace vcsJACOBI
{
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
        }
        //--------------------------------------------------------------------
        private void toolStripButton1_Click(object sender, EventArgs e)
        {
            //n：固有方程式の次数(配列宣言はn-1)
            //ev[k]：固有値(大きい順に出力)
            //vec[i=0〜n-1][k]：固有ベクトル
            string dat;
            int i, j, k;
            int n = 3;
            double[,] a = new double[n, n];
            double[] ev = new double[n];
            double[,] vec = new double[n, n];
            double[,] ar = new double[n, n];
            double[,] c = new double[n, n];

            //Case-1
            a[0, 0] = 4; a[0, 1] = 2; a[0, 2] = 3;
            a[1, 0] = 2; a[1, 1] = 6; a[1, 2] = 1;
            a[2, 0] = 3; a[2, 1] = 1; a[2, 2] = 5;

            for (i = 0; i <= n - 1; i++)
            {
                for (j = 0; j <= n - 1; j++)
                {
                    ar[i, j] = a[i, j];
                }
            }

            JACOBI(n, a, ev, vec);

            Console.WriteLine("*** Eigen value ***");
            dat = String.Format("{0,16:0.000000e+000}", ev[0]);
            for (j = 1; j <= n - 1; j++)
            {
                dat = dat + String.Format("{0,16:0.000000e+000}", ev[j]);
            }
            Console.WriteLine(dat);
            Console.WriteLine("*** Eigen vector ***");
            for (i = 0; i <= n - 1; i++)
            {
                dat = String.Format("{0,16:0.000000e+000}", vec[i, 0]);
                for (j = 1; j <= n - 1; j++)
                {
                    dat = dat + String.Format("{0,16:0.000000e+000}", vec[i, j]);
                }
                Console.WriteLine(dat);
            }

            //Case-1 出力結果
            //*** Eigen value ***
            //9.000000e+000   4.732051e+000   1.267949e+000
            //*** Eigen vector ***
            //5.773503e-001   2.113249e-001   7.886751e-001
            //5.773503e-001  -7.886751e-001  -2.113249e-001
            //5.773503e-001   5.773503e-001  -5.773503e-001

            //検算
            for (i = 0; i <= n - 1; i++)
            {
                for (j = 0; j <= n - 1; j++)
                {
                    c[i, j] = 0.0;
                    for (k = 0; k <= n - 1; k++)
                    {
                        c[i, j] = c[i, j] + ar[i, k] * vec[k, j];
                    }
                }
            }
            Console.WriteLine("*** 検算 ***");
            for (i = 0; i <= n - 1; i++)
            {
                dat = String.Format("{0,16:0.000000e+000}", c[i, 0]);
                for (j = 1; j <= n - 1; j++)
                {
                    dat = dat + String.Format("{0,16:0.000000e+000}", c[i, j]);
                }
                dat = dat + "  ";
                dat = dat + String.Format("{0,16:0.000000e+000}", vec[i, 0] * ev[0]);
                for (j = 1; j <= n - 1; j++)
                {
                    dat = dat + String.Format("{0,16:0.000000e+000}", vec[i, j] * ev[j]);
                }
                Console.WriteLine(dat);
            }

        }
        //--------------------------------------------------------------------
        private void JACOBI(int n, double[,] a, double[] ev, double[,] vec)
        {
            //Jacobi法による固有値・固有ベクトル計算（奥村先生）
            //固有値ev[k],固有ベクトルvec[i,k]
            int i, j, k, im = 0;
            double s, c, aoff, eps, t, aa = 0.0, aj, ak, vj, vk, vmax, bb;

            for (i = 0; i <= n - 1; i++)
            {
                for (j = 0; j <= n - 1; j++)
                {
                    vec[i, j] = 0.0;
                }
                vec[i, i] = 1.0;
            }
            aoff = 0.0;
            for (j = 0; j <= n - 2; j++)
            {
                for (k = j + 1; k <= n - 1; k++)
                {
                    aoff = aoff + a[j, k] * a[j, k];
                }
            }
            eps = aoff * 2.0;
            for (j = 0; j <= n - 1; j++)
            {
                eps = eps + a[j, j] * a[j, j];
            }
            eps = eps * 0.00000000005;

            while (aoff > eps)
            {
                for (j = 0; j <= n - 2; j++)
                {
                    for (k = j + 1; k <= n - 1; k++)
                    {
                        if (a[j, k] == 0.0)
                        {
                            t = 0.0;
                        }
                        else
                        {
                            aa = (a[k, k] - a[j, j]) / (2.0 * a[j, k]);
                        }
                        if (aa >= 0)
                        {
                            t = 1.0 / (aa + Math.Sqrt(1.0 + aa * aa));
                        }
                        else
                        {
                            t = 1.0 / (aa - Math.Sqrt(1.0 + aa * aa));
                        }
                        c = 1.0 / Math.Sqrt(1.0 + t * t);
                        s = t * c;
                        a[j, j] = a[j, j] - t * a[j, k];
                        a[k, k] = a[k, k] + t * a[j, k];
                        a[j, k] = 0.0;
                        for (i = 0; i <= j - 1; i++)
                        {
                            aj = a[i, j];
                            ak = a[i, k];
                            a[i, j] = aj * c - ak * s;
                            a[i, k] = aj * s + ak * c;
                        }
                        for (i = j + 1; i <= k - 1; i++)
                        {
                            aj = a[j, i];
                            ak = a[i, k];
                            a[j, i] = aj * c - ak * s;
                            a[i, k] = aj * s + ak * c;
                        }
                        for (i = k + 1; i <= n - 1; i++)
                        {
                            aj = a[j, i];
                            ak = a[k, i];
                            a[j, i] = aj * c - ak * s;
                            a[k, i] = aj * s + ak * c;
                        }
                        for (i = 0; i <= n - 1; i++)
                        {
                            vj = vec[i, j];
                            vk = vec[i, k];
                            vec[i, j] = vj * c - vk * s;
                            vec[i, k] = vj * s + vk * c;
                        }
                    }
                }
                aoff = 0.0;
                for (j = 0; j <= n - 2; j++)
                {
                    for (k = j + 1; k <= n - 1; k++)
                    {
                        aoff = aoff + a[j, k] * a[j, k];
                    }
                }
            }

            for (k = 0; k <= n - 1; k++)
            {
                ev[k] = a[k, k];
            }
            for (k = 0; k <= n - 2; k++)
            {
                vmax = -1.0E+30;
                for (i = k; i <= n - 1; i++)
                {
                    if (ev[i] >= vmax)
                    {
                        im = i;
                        vmax = ev[i];
                    }
                }
                aa = ev[k]; ev[k] = ev[im]; ev[im] = aa;
                for (i = 0; i <= n - 1; i++)
                {
                    bb = vec[i, k];
                    vec[i, k] = vec[i, im];
                    vec[i, im] = bb;
                }
            }
        }
        //--------------------------------------------------------------------
    }
}