C#，码海拾贝（31）——约化“对称矩阵“为“对称三对角阵“的“豪斯荷尔德Householder变换法“之C#源代码，《C#数值计算算法编程》源代码升级改进版

csharp
using System;

namespace HouseholderTransformation
{
    class Program
    {
        static void Main(string[] args)
        {
            // 定义对称矩阵
            double[] symmetricMatrix = {
                { 4 1 -2 }
                { 1 2 3 }
                { -2 3 5 }
            };

            // 获取矩阵的阶数
            int n = symmetricMatrix.GetLength(0);

            // 初始化单位矩阵
            double[] identityMatrix = new double[n n];
            for (int i = 0; i < n; i++)
            {
                identityMatrix[i i] = 1;
            }

            // 对称矩阵约化为对称三对角阵
            for (int k = 0; k < n - 2; k++)
            {
                double[] x = new double[n - k];
                double[] v = new double[n - k];
                double[] w = new double[n - k];

                for (int i = k; i < n; i++)
                {
                    x[i - k] = symmetricMatrix[i k];
                }

                double normX = 0;
                for (int i = 0; i < n - k; i++)
                {
                    normX += x[i] * x[i];
                }
                normX = Math.Sqrt(normX);

                v[0] = x[0] + Math.Sign(x[0]) * normX;
                for (int i = 1; i < n - k; i++)
                {
                    v[i] = x[i];
                }

                double normV = 0;
                for (int i = 0; i < n - k; i++)
                {
                    normV += v[i] * v[i];
                }
                normV = Math.Sqrt(normV);

                for (int i = 0; i < n - k; i++)
                {
                    v[i] /= normV;
                }

                for (int i = 0; i < n - k; i++)
                {
                    for (int j = 0; j < n - k; j++)
                    {
                        w[i] += v[j] * symmetricMatrix[k + j k + i];
                    }
                }

                for (int i = 0; i < n - k; i++)
                {
                    for (int j = 0; j < n - k; j++)
                    {
                        symmetricMatrix[k + i k + j] -= 2 * v[i] * w[j];
                    }
                }

                for (int i = 0; i < n - k; i++)
                {
                    for (int j = 0; j < n; j++)
                    {
                        w[i] = 0;
                        for (int l = 0; l < n - k; l++)
                        {
                            w[i] += v[l] * identityMatrix[k + l j];
                        }
                    }
                }

                for (int i = 0; i < n - k; i++)
                {
                    for (int j = 0; j < n; j++)
                    {
                        identityMatrix[k + i j] -= 2 * v[i] * w[j];
                    }
                }
            }

            // 输出对称三对角阵
            Console.WriteLine(对称三对角阵：);
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    Console.Write(symmetricMatrix[i j] +  );
                }
                Console.WriteLine();
            }
        }
    }
}