Я конвертирую программу на Фортране в C #. Это должно быть сделано по частям, с подтверждением концепций на этом пути.
Один из этих начальных шагов - репликация функций IMSL, которые он использует. К счастью, он использует только несколько избранных: некоторые тривиальные генерации случайных чисел, некоторые тривиальные инверсии нормального распределения, но также и один не совсем тривиальный: MCHOL.
Из документации :
Вычисляет верхнюю треугольную факторизацию вещественной симметричной матрицы A плюс диагональную матрицу D, где D определяется последовательно во время факторизации Холецкого, чтобы сделать A + D неотрицательным определенным.
Рутинный MCHOL вычисляет факторизацию Холецкого, RTR, для A + D, где A симметрична, а D является диагональной матрицей с достаточно большими диагональными элементами, так что A + D неотрицательно определена. Процедура аналогична описанной Джиллом, Мюрреем и Райтом (1981, стр. 108-111). Здесь, однако, мы допускаем, чтобы A + D было единичным.
(Более подробная информация и образец по ссылке).
Для подтверждения концепции мне нужно иметь возможность воспроизвести результаты, как указано в примере для документации MCHOL: передача в эту матрицу из образца:
DATA (A(1,J),J=1,N)/36.0, 12.0, 30.0, 6.0, 18.0/
DATA (A(2,J),J=1,N)/12.0, 20.0, 2.0, 10.0, 22.0/
DATA (A(3,J),J=1,N)/30.0, 2.0, 29.0, 1.0, 7.0/
DATA (A(4,J),J=1,N)/6.0, 10.0, 1.0, 14.0, 20.0/
DATA (A(5,J),J=1,N)/8.0, 22.0, 7.0, 20.0, 40.0/
И получить взамен следующее:
6.000 2.000 5.000 1.000 3.000
0.000 4.000 -2.000 2.000 4.000
0.000 0.000 0.000 0.000 0.000
0.000 0.000 0.000 3.000 3.000
0.000 0.000 0.000 0.000 2.449
До сих пор я пытался использовать Math.NET , но он не будет работать на этом примере матрицы, потому что он не определен положительно.
Я также пробовал части ALGLIB, в частности spdmatrixcholesky . Вроде работает, но только для части матрицы:
6 2 5 1 3
12 4 -2 2 4
30 2 0 1 7
6 10 1 14 20
8 22 7 20 40
Кто-нибудь знает, что я здесь не так делаю? Нужно ли здесь вызывать другую функцию?
Поскольку в карточках, кажется, нет быстрого ответа, было бы лучше, если бы я понял основную математику, чтобы хотя бы попытаться выяснить, что здесь происходит. Любое теоретическое обоснование или указатели также приветствуются.