TNT::Array2D<double> A(r,r);
    TNT::Array2D<double> L;
    TNT::Array2D<double> U;
    TNT::Array1D<double> X;
    TNT::Array1D<double> B1(r);
    TNT::Array1D<double> B2(r);
    TNT::Array1D<double> B3(r);
    TNT::Array1D<double> B4(r);
    TNT::Array1D<double> Z1;
    TNT::Array1D<double> Z2;
    TNT::Array1D<double> Z3;
    TNT::Array1D<double> Z4;

    A[0][0] = -116.66654;
    A[1][1] = -116.66654;
    A[2][2] = 133.33383;
    A[3][3] = 50.000125;
    A[4][4] = 50.000125;

    A[0][1] = 583.33346;
    A[0][2] = -333.33308;
    A[0][3] = 100.00012;
    A[0][4] = 100.00012;

    A[1][0] = 583.33346;
    A[1][2] = -333.33308;
    A[1][3] = 100.00012;
    A[1][4] = 100.00012;

    A[2][0] = -333.33308;
    A[2][1] = -333.33308;
    A[2][3] = 200.00025;
    A[2][4] = 200.00025;

    A[3][0] = 100.00012;
    A[3][1] = 100.00012;
    A[3][2] = 200.00025;
    A[3][4] = -649.99988;

    A[4][0] = 100.00012;
    A[4][1] = 100.00012;
    A[4][2] = 200.00025;
    A[4][3] = -649.99988;

    B1[0] = -0.33388066;
    B1[1] = 1.08033290;
    B1[2] = -0.98559856;
    B1[3] = 1.31947922;
    B1[4] = -0.09473435;

    B2[0] = -0.33388066;
    B2[1] = 1.08033290;
    B2[2] = -0.98559855;
    B2[3] = 1.32655028;
    B2[4] = -0.10180541;

    B3[0] = 0.72677951;
    B3[1] = 0.72677951;
    B3[2] = -0.27849178;
    B3[3] = 0.96592583;
    B3[4] = 0.96592583;

    B4[0] = 0.73031505;
    B4[1] = 0.73031505;
    B4[2] = -0.27142071;
    B4[3] = 0.96946136;
    B4[4] = 0.96946136;

    JAMA::LU<double> lu(A);

    L = lu.getL();
    U = lu.getU();

    Z1 = lu.solve(B1);
    Z2 = lu.solve(B2);
    Z3 = lu.solve(B3);
    Z4 = lu.solve(B4);