Пример реализации LDLT-разложения Холецкого

C#

///
///<summary>
///Преобразование матрицы А в LDLT вид
///Автор: Михайловский Е.А. 07.04.2014
///</summary>
///<param name="A">Исходная матрица</param>
///<param name="L">Нижняя треугольная матрица с 1 на диагонали</param>
///<param name="D">Массив диагональных элементов</param>

void LDLT(double[][] A, out double[,] L, out double[] D)
{
 int n = A.GetLength(0);
 D = new double[n];
 L = new double[n, n];
 int i, j, k;
 double sum = 0;
 for (i = 0; i < n; i++)//столбцы
  for (j = i; j < n; j++)//строки
  {
    sum = A[j][i];//значение вычисляемого элемента
    for (k = 0; k < i; k++)//вычитание элементов строки из вычисляемого элемента
     sum = sum - L[i, k] * D[k] * L[j, k];
    if (i == j)//диагональный элемент
    {
     if (sum <= 0) new Exception("A is not positive deﬁnite");
     D[i] = sum;//диагональный элемент
     L[i, i] = 1;//диагональ
    }
    else L[j, i] = sum / D[i];//внедиагональный элемент
  }
}