21#include <ql/math/matrixutilities/choleskydecomposition.hpp>
22#include <ql/math/comparison.hpp>
27 Size i, j, size =
S.rows();
29 QL_REQUIRE(size ==
S.columns(),
30 "input matrix is not a square matrix");
31 #if defined(QL_EXTRA_SAFETY_CHECKS)
32 for (i=0; i<
S.rows(); i++)
34 QL_REQUIRE(
S[i][j] ==
S[j][i],
35 "input matrix is not symmetric");
38 Matrix result(size, size, 0.0);
40 for (i=0; i<size; i++) {
41 for (j=i; j<size; j++) {
44 sum -= result[i][k]*result[j][k];
47 QL_REQUIRE(flexible || sum > 0.0,
48 "input matrix is not positive definite");
51 result[i][i] = std::sqrt(std::max<Real>(sum, 0.0));
58 :
Real(sum / result[i][i]);
Matrix used in linear algebra.
QL_INTEGER Integer
integer number
std::size_t Size
size of a container
Matrix CholeskyDecomposition(const Matrix &S, bool flexible)
bool close_enough(const Quantity &m1, const Quantity &m2, Size n)