This repository has been archived by the owner on May 11, 2023. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 7
/
LUD.cs
116 lines (109 loc) · 3.43 KB
/
LUD.cs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
using RCNet.Extensions;
using RCNet.MathTools.VectorMath;
using System;
namespace RCNet.MathTools.MatrixMath
{
/// <summary>
/// Implements the LU (Lowed-Upper) decomposition of a square matrix.
/// </summary>
/// <remarks>
/// <para>
/// This class is based on a class from the public domain JAMA package.
/// </para>
/// <para>
/// http://math.nist.gov/javanumerics/jama/
/// </para>
/// </remarks>
public class LUD
{
//Attributes
private readonly int _size;
private readonly double[][] _lu;
//Constructor
/// <summary>
/// Creates an initialized instance.
/// </summary>
/// <param name="matrix">The square matrix.</param>
public LUD(Matrix matrix)
{
//Check the matrix is the square matrix
if (!matrix.IsSquareMatrix)
{
throw new ArgumentException("Matrix is not square.", "matrix");
}
//Initialization
_size = matrix.NumOfRows;
//LU data
_lu = new double[_size][];
for (int i = 0; i < _size; i++)
{
_lu[i] = new double[_size];
_lu[i].Populate(0);
}
//Decomposition
for (int i = 0; i < _size; i++)
{
for (int j = i; j < _size; j++)
{
double sum = 0;
for (int k = 0; k < i; k++)
{
sum += _lu[i][k] * _lu[k][j];
}
_lu[i][j] = matrix.Data[i][j] - sum;
}
for (int j = i + 1; j < _size; j++)
{
double sum = 0;
for (int k = 0; k < i; k++)
{
sum += _lu[j][k] * _lu[k][i];
}
_lu[j][i] = (1d / _lu[i][i]) * (matrix.Data[j][i] - sum);
}
}
return;
}
//Properties
/// <summary>
/// LU data
/// </summary>
public double[][] LU { get { return _lu; } }
//Methods
/// <summary>
/// Solves system of linear equations A*x = b.
/// </summary>
/// <param name="b">Vector of desired results (b).</param>
/// <returns>Vector of computed linear coefficients (x).</returns>
public Vector Solve(Vector b)
{
if (_size != b.Length)
{
throw new ArgumentException($"Invalid length ({b.Length}) of specified vector. Expected length {_size}.", "b");
}
double[] bData = b.Data;
//Find solution
double[] y = new double[_size];
for (int i = 0; i < _size; i++)
{
double sum = 0;
for (int k = 0; k < i; k++)
{
sum += _lu[i][k] * y[k];
}
y[i] = bData[i] - sum;
}
double[] x = new double[_size];
for (int i = _size - 1; i >= 0; i--)
{
double sum = 0;
for (int k = i + 1; k < _size; k++)
{
sum += _lu[i][k] * x[k];
}
x[i] = (1d / _lu[i][i]) * (y[i] - sum);
}
return new Vector(x, false);
}
}//LUD
}//Namespace