对称上三角稀疏矩阵非零元素的一维存储及LU分解法求解的C#代码

以下是对称上三角稀疏矩阵非零元素的一维存储及LU分解法求解的C#代码示例：

using System;

namespace SymmetricUpperTriangularSparseMatrix
{
    public class Matrix
    {
        int n;
        int[] irow, jcol;
        double[] val;

        public Matrix(int n, int[] irow, int[] jcol, double[] val)
        {
            this.n = n;
            this.irow = irow;
            this.jcol = jcol;
            this.val = val;
        }

        public void LUDecompose(out Matrix L, out Matrix U)
        {
            int[] il = new int[n + 1];
            int[] iu = new int[n + 1];
            int[] jl = new int[val.Length];
            int[] ju = new int[val.Length];
            double[] al = new double[val.Length];
            double[] au = new double[val.Length];

            for (int i = 0; i < n + 1; i++)
            {
                il[i] = -1;
                iu[i] = -1;
            }

            for (int k = 0; k < n; k++)
            {
                il[k] = k;
                iu[k] = k;
                int ll = -1, lu = -1;
                for (int p = irow[k]; p < irow[k + 1]; p++)
                {
                    int j = jcol[p];
                    if (j < k)
                    {
                        if (il[k] == k) il[k] = p;
                        jl[++ll] = j;
                        al[ll] = val[p];
                    }
                    else if (j == k)
                    {
                        if (il[k] == k) il[k] = p;
                        if (iu[k] == k) iu[k] = p;
                        al[k] = val[p];
                        au[k] = 1.0;
                    }
                    else
                    {
                        ju[++lu] = j;
                        au[lu] = val[p];
                    }
                }
                for (int p = 0; p <= ll; p++)
                {
                    int i = jl[p];
                    double lki = al[p] / al[k];
                    for (int q = iu[i]; q <= irow[i + 1] - 1; q++)
                    {
                        int j = jcol[q];
                        if (j > k) break;
                        int r = j == k ? k : il[k];
                        if (il[j] < 0)
                        {
                            jl[++ll] = j;
                            al[ll] = -lki * val[q];
                            il[j] = ll;
                        }
                        else
                        {
                            al[il[j]] -= lki * val[q];
                        }
                        if (il[k] < 0)
                        {
                            il[k] = lu + 1;
                            iu[k + 1] = lu + 1;
                            ju[++lu] = k;
                            au[lu] = al[k];
                        }
                        au[r] -= lki * au[q];
                    }
                }
            }

            int[] ilu = new int[n + 1];
            int[] jlu = new int[val.Length];
            double[] valL = new double[val.Length];
            double[] valU = new double[val.Length];

            for (int i = 0; i < n + 1; i++)
            {
                ilu[i] = -1;
            }

            int cntL = 0, cntU = 0;
            for (int i = 0; i < n; i++)
            {
                for (int p = il[i]; p <= il[i + 1] - 1; p++)
                {
                    int j = jl[p];
                    if (j < i)
                    {
                        jlu[cntL] = j;
                        valL[cntL] = al[p];
                        if (ilu[i] < 0) ilu[i] = cntL;
                        cntL++;
                    }
                    else if (j == i)
                    {
                        valL[cntL] = 1.0;
                        if (ilu[i] < 0) ilu[i] = cntL;
                        cntL++;
                    }
                    else
                    {
                        jlu[cntU] = j;
                        valU[cntU] = au[p];
                        if (ilu[i] < 0) ilu[i] = cntU;
                        cntU++;
                    }
                }
                for (int p = iu[i]; p <= iu[i + 1] - 1; p++)
                {
                    int j = ju[p];
                    if (j < i)
                    {
                        jlu[cntL] = j;
                        valL[cntL] = al[p];
                        if (ilu[i] < 0) ilu[i] = cntL;
                        cntL++;
                    }
                    else if (j == i)
                    {
                        valL[cntL] = 1.0;
                        if (ilu[i] < 0) ilu[i] = cntL;
                        cntL++;
                    }
                    else
                    {
                        jlu[cntU] = j;
                        valU[cntU] = au[p];
                        if (ilu[i] < 0) ilu[i] = cntU;
                        cntU++;
                    }
                }
            }

            L = new Matrix(n, ilu, jlu, valL);
            U = new Matrix(n, ilu, jlu, valU);
        }

        public double[] LUSolve(double[] b)
        {
            Matrix L, U;
            LUDecompose(out L, out U);

            double[] y = new double[n];
            for (int i = 0; i < n; i++)
            {
                y[i] = b[i];
                for (int p = L.irow[i]; p < L.irow[i + 1]; p++)
                {
                    int j = L.jcol[p];
                    if (j < i) y[i] -= L.val[p] * y[j];
                }
            }

            double[] x = new double[n];
            for (int i = n - 1; i >= 0; i--)
            {
                x[i] = y[i];
                for (int p = U.irow[i]; p < U.irow[i + 1]; p++)
                {
                    int j = U.jcol[p];
                    if (j > i) x[i] -= U.val[p] * x[j];
                }
                x[i] /= U.val[U.irow[i]];
            }

            return x;
        }
    }

    class Program
    {
        static void Main(string[] args)
        {
            int n = 3;
            int[] irow = { 0, 1, 2, 3 };
            int[] jcol = { 0, 1, 2, 1, 2, 2 };
            double[] val = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
            Matrix A = new Matrix(n, irow, jcol, val);

            double[] b = { 1.0, 2.0, 3.0 };
            double[] x = A.LUSolve(b);

            Console.WriteLine("Solution:");
            for (int i = 0; i < n; i++)
            {
                Console.WriteLine(x[i]);
            }
        }
    }
}

该代码中定义了一个 Matrix 类，其构造函数输入稀疏矩阵的维数、非零元素的行列下标及值。LUDecompose 方法实现了对称上三角稀疏矩阵的 LU 分解，并输出分解后的 L 和 U 矩阵。LUSolve 方法实现了使用 LU 分解求解线性方程组。最后的 Main 方法给出了一个样例，其输入的稀疏矩阵是一个 $3 \times 3$ 的矩阵，用 LU 分解求解线性方程组。