鸟语天空

克洛脱(Crout)LU分解——C#实现

作者：追风剑情 发布于：2018-10-14 15:10 分类：Algorithms

克洛脱(Crout)LU分解原理参见 克洛脱(Crout)矩阵分解——LU分解

using System;
using System.Text;

namespace ConsoleApp1
{
    class Program
    {
        static void Main(string[] args)
        {
            Console.WriteLine("克洛脱分解测试");
            //测试数据
            const int N = 4;
            float[] A = new float[N*N]
            {
                2, 9, 0, 0,
                6, 2, 2, 0,
                0, 8, 2, 3,
                0, 0, 10, 2
            };
            float[] L;
            float[] U;

            CroutLU(N, A, out L, out U);

            Console.WriteLine("A:");
            PrintM(N, A);
            Console.WriteLine("L:");
            PrintM(N, L);
            Console.WriteLine("U:");
            PrintM(N, U);

            Console.WriteLine("注: 因为这里用的是float型数据,计算时会有一定的精度丢失.");

            Console.Read();
        }

        /**
         * 克洛脱(Crout)LU分解
         * @param N NxN矩阵
         * @param A 原矩阵
         * @param L L矩阵
         * @param U U矩阵
         */
        static void CroutLU(int N, float[] A, out float[] L, out float[] U)
        {
            L = new float[N*N];
            U = new float[N*N];
            float[] a_arr = new float[N];
            float[] b_arr = new float[N];
            float[] c_arr = new float[N];
            float[] l_arr = new float[N];
            float[] m_arr = new float[N];
            float[] u_arr = new float[N];

            //a_arr[]存放原矩阵主对角线数据
            for (int i=0; i<N; i++)
            {
                a_arr[i] = A[i * N + i];
            }

            //b_arr[]存放原矩阵下对角线数据
            for(int i=1; i<N; i++)
            {
                b_arr[i] = A[i * N + i - 1];
            }

            //c_arr[]存放原矩阵上对角线数据
            for (int i=0; i<N-1; i++)
            {
                c_arr[i] = A[i * N + i + 1];
            }

            l_arr[0] = a_arr[0];
            u_arr[0] = c_arr[0] / l_arr[0];

            //根据公式求出其余l_arr[]、m_arr[]、u_arr[]的值
            for(int i=1; i<N; i++)
            {
                m_arr[i] = b_arr[i];
                l_arr[i] = a_arr[i] - m_arr[i] * u_arr[i - 1];
                if (i < N - 1)
                    u_arr[i] = c_arr[i] / l_arr[i];
            }

            //给L矩阵赋值
            L[0] = l_arr[0];
            for(int i=1; i<N; i++)
            {
                //给主对角线赋值
                L[i * N + i] = l_arr[i];
                //给下对角线赋值
                L[i * N + i - 1] = m_arr[i];
            }
            //给U矩阵赋值
            for (int i = 0; i < N-1; i++)
            {
                //给主对角线赋值
                U[i * N + i] = 1;
                //给上对角线赋值
                U[i * N + i + 1] = u_arr[i];
            }
            U[N*N - 1] = 1;
        }

        static void PrintM(int N, float[] M)
        {
            StringBuilder sb = new StringBuilder();
            for(int i=0; i<N; i++)
            {
                for (int j = 0; j < N; j++)
                    sb.Append(M[i * N + j] +"   ");
                sb.AppendLine();
            }
            Console.WriteLine(sb.ToString());
        }
    }
}

运行测试

标签: Algorithms