LU.cs

usingSystem;

namespaceAlgorithms.Numeric.Decomposition;

/// <summary>
/// LU-decomposition factors the "source" matrix as the product of lower triangular matrix
/// and upper triangular matrix.
/// </summary>
publicstaticclassLu
{
/// <summary>
/// Performs LU-decomposition on "source" matrix.
/// Lower and upper matrices have same shapes as source matrix.
/// Note: Decomposition can be applied only to square matrices.
/// </summary>
/// <param name="source">Square matrix to decompose.</param>
/// <returns>Tuple of lower and upper matrix.</returns>
/// <exception cref="ArgumentException">Source matrix is not square shaped.</exception>
publicstatic(double[,]L,double[,]U)Decompose(double[,]source)
{
if(source.GetLength(0)!=source.GetLength(1))
{
thrownewArgumentException("Source matrix is not square shaped.");
}

varpivot=source.GetLength(0);
varlower=newdouble[pivot,pivot];
varupper=newdouble[pivot,pivot];

for(vari=0;i<pivot;i++)
{
for(vark=i;k<pivot;k++)
{
doublesum=0;

for(varj=0;j<i;j++)
{
sum+=lower[i,j]*upper[j,k];
}

upper[i,k]=source[i,k]-sum;
}

for(vark=i;k<pivot;k++)
{
if(i==k)
{
lower[i,i]=1;
}
else
{
doublesum=0;

for(varj=0;j<i;j++)
{
sum+=lower[k,j]*upper[j,i];
}

lower[k,i]=(source[k,i]-sum)/upper[i,i];
}
}
}

return(L:lower,U:upper);
}

/// <summary>
/// Eliminates linear equations system represented as A*x=b, using LU-decomposition,
/// where A - matrix of equation coefficients, b - vector of absolute terms of equations.
/// </summary>
/// <param name="matrix">Matrix of equation coefficients.</param>
/// <param name="coefficients">Vector of absolute terms of equations.</param>
/// <returns>Vector-solution for linear equations system.</returns>
/// <exception cref="ArgumentException">Matrix of equation coefficients is not square shaped.</exception>
publicstaticdouble[]Eliminate(double[,]matrix,double[]coefficients)
{
if(matrix.GetLength(0)!=matrix.GetLength(1))
{
thrownewArgumentException("Matrix of equation coefficients is not square shaped.");
}

varpivot=matrix.GetLength(0);
varupperTransform=newdouble[pivot,1];// U * upperTransform = coefficients
varsolution=newdouble[pivot];// L * solution = upperTransform
(double[,]l,double[,]u)=Decompose(matrix);

for(vari=0;i<pivot;i++)
{
doublepivotPointSum=0;

for(varj=0;j<i;j++)
{
pivotPointSum+=upperTransform[j,0]*l[i,j];
}

upperTransform[i,0]=(coefficients[i]-pivotPointSum)/l[i,i];
}

for(vari=pivot-1;i>=0;i--)
{
doublepivotPointSum=0;

for(varj=i;j<pivot;j++)
{
pivotPointSum+=solution[j]*u[i,j];
}

solution[i]=(upperTransform[i,0]-pivotPointSum)/u[i,i];
}

returnsolution;
}
}