using System;

namespace Zhou.CSharp.Algorithm
{
    /// <summary>
    /// 矩阵类
    /// 作者：周长发
    /// 改进：深度混淆
    /// https://blog.csdn.net/beijinghorn
    /// </summary>
    public partial class Matrix
    {

        /// <summary>
        /// 求赫申伯格矩阵全部特征值的QR方法
        /// </summary>
        /// <param name="src">源矩阵</param>
        /// <param name="dblU">一维数组，长度为矩阵的阶数，返回时存放特征值的实部</param>
        /// <param name="dblV">一维数组，长度为矩阵的阶数，返回时存放特征值的虚部</param>
        /// <param name="nMaxIt">迭代次数</param>
        /// <param name="eps">计算精度</param>
        /// <returns>求解是否成功</returns>
        public static bool ComputeEvHBerg(Matrix src, out double[] dblU, out double[] dblV, int nMaxIt = 100, double eps = 1.0E-7)
        {
            int m, it, i, j, k, ii, jj, kk;
            double b, c, w, g, xy, p, q, r, x, s, e, f, z, y;

            int n = src.Columns;
            dblU = new double[n];
            dblV = new double[n];

            it = 0;
            m = n;
            while (m != 0)
            {
                int u = m - 1;
                while ((u > 0) && (Math.Abs(src[u * n + u - 1]) >
                    eps * (Math.Abs(src[(u - 1) * n + u - 1]) + Math.Abs(src[u * n + u]))))
                {
                    u = u - 1;
                }
                ii = (m - 1) * n + m - 1;
                jj = (m - 1) * n + m - 2;
                kk = (m - 2) * n + m - 1;
                int uu = (m - 2) * n + m - 2;
                if (u == m - 1)
                {
                    dblU[m - 1] = src[(m - 1) * n + m - 1];
                    dblV[m - 1] = 0.0;
                    m = m - 1;
                    it = 0;
                }
                else if (u == m - 2)
                {
                    b = -(src[ii] + src[uu]);
                    c = src[ii] * src[uu] - src[jj] * src[kk];
                    w = b * b - 4.0 * c;
                    if (Math.Abs(w) < float.Epsilon)
                    {
                        return false;
                    }
                    y = Math.Sqrt(Math.Abs(w));
                    if (w > 0.0)
                    {
                        xy = 1.0;
                        if (b < 0.0)
                        {
                            xy = -1.0;
                        }
                        dblU[m - 1] = (-b - xy * y) / 2.0;
                        dblU[m - 2] = c / dblU[m - 1];
                        dblV[m - 1] = 0.0; dblV[m - 2] = 0.0;
                    }
                    else
                    {
                        dblU[m - 1] = -b / 2.0;
                        dblU[m - 2] = dblU[m - 1];
                        dblV[m - 1] = y / 2.0;
                        dblV[m - 2] = -dblV[m - 1];
                    }

                    m = m - 2;
                    it = 0;
                }
                else
                {
                    if (it >= nMaxIt)
                    {
                        return false;
                    }
                    it = it + 1;
                    for (j = u + 2; j <= m - 1; j++)
                    {
                        src[j * n + j - 2] = 0.0;
                    }
                    for (j = u + 3; j <= m - 1; j++)
                    {
                        src[j * n + j - 3] = 0.0;
                    }
                    for (k = u; k <= m - 2; k++)
                    {
                        if (k != u)
                        {
                            p = src[k * n + k - 1];
                            q = src[(k + 1) * n + k - 1];
                            r = 0.0;
                            if (k != m - 2)
                            {
                                r = src[(k + 2) * n + k - 1];
                            }
                        }
                        else
                        {
                            x = src[ii] + src[uu];
                            y = src[uu] * src[ii] - src[kk] * src[jj];
                            ii = u * n + u;
                            jj = u * n + u + 1;
                            kk = (u + 1) * n + u;
                            uu = (u + 1) * n + u + 1;
                            p = src[ii] * (src[ii] - x) + src[jj] * src[kk] + y;
                            q = src[kk] * (src[ii] + src[uu] - x);
                            r = src[kk] * src[(u + 2) * n + u + 1];
                        }

                        if ((Math.Abs(p) + Math.Abs(q) + Math.Abs(r)) > float.Epsilon)//  != 0.0)
                        {
                            xy = 1.0;
                            if (p < 0.0)
                            {
                                xy = -1.0;
                            }
                            s = xy * Math.Sqrt(p * p + q * q + r * r);
                            if (k != u)
                            {
                                src[k * n + k - 1] = -s;
                            }
                            e = -q / s;
                            f = -r / s;
                            x = -p / s;
                            y = -x - f * r / (p + s);
                            g = e * r / (p + s);
                            z = -x - e * q / (p + s);
                            for (j = k; j <= m - 1; j++)
                            {
                                ii = k * n + j;
                                jj = (k + 1) * n + j;
                                p = x * src[ii] + e * src[jj];
                                q = e * src[ii] + y * src[jj];
                                r = f * src[ii] + g * src[jj];
                                if (k != m - 2)
                                {
                                    kk = (k + 2) * n + j;
                                    p = p + f * src[kk];
                                    q = q + g * src[kk];
                                    r = r + z * src[kk];
                                    src[kk] = r;
                                }

                                src[jj] = q; src[ii] = p;
                            }

                            j = k + 3;
                            if (j >= m - 1)
                            {
                                j = m - 1;
                            }
                            for (i = u; i <= j; i++)
                            {
                                ii = i * n + k;
                                jj = i * n + k + 1;
                                p = x * src[ii] + e * src[jj];
                                q = e * src[ii] + y * src[jj];
                                r = f * src[ii] + g * src[jj];
                                if (k != m - 2)
                                {
                                    kk = i * n + k + 2;
                                    p = p + f * src[kk];
                                    q = q + g * src[kk];
                                    r = r + z * src[kk];
                                    src[kk] = r;
                                }

                                src[jj] = q;
                                src[ii] = p;
                            }
                        }
                    }
                }
            }

            return true;
        }
 

    }
}