克洛脱(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());
}
}
}
运行测试