Changeset 3856


Ignore:
Timestamp:
Mar 5, 2015, 10:58:06 AM (7 years ago)
Author:
guite
Message:

matrix: using permutation and LU decomposition for determinant and inverse computing

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/lol/math/matrix.h

    r3855 r3856  
    548548T determinant(mat_t<T, N, N> const &m)
    549549{
    550 #if 0 /* XXX: LU decomposition is currently broken */
    551550    mat_t<T, N, N> L, U;
    552     lu_decomposition(m, L, U);
     551    vec_t<int, N> P = p_vector(m);
     552    lu_decomposition(permute_rows(m, P), L, U);
    553553
    554554    T det = 1;
     
    556556        det *= U[i][i];
    557557
    558     return det;
    559 #else
    560     T ret = (T)0;
    561     for (int i = 0; i < N; ++i)
    562         ret += m[i][0] * cofactor(m, i, 0);
    563     return ret;
    564 #endif
     558    return permutation_det(P) * det;
    565559}
    566560
     
    593587        for (int j = 0 ; j < N ; ++j)
    594588        {
    595             if (!used[j] && (index_max == -1 || m[i][j] > m[i][index_max]))
     589            if (!used[j] && (index_max == -1 || abs(m[i][j]) > abs(m[i][index_max])))
    596590                index_max = j;
    597591        }
     
    627621}
    628622
     623template<typename T, int N>
     624mat_t<T, N, N> permute_cols(mat_t<T, N, N> const & m, vec_t<int, N> const & permutation)
     625{
     626    mat_t<T, N, N> result;
     627
     628    for (int i = 0 ; i < N ; ++i)
     629    {
     630        for (int j = 0 ; j < N ; ++j)
     631        {
     632            result[i][j] = m[permutation[i]][j];
     633        }
     634    }
     635
     636    return result;
     637}
     638
     639template<int N>
     640vec_t<int, N> p_transpose(vec_t<int, N> P)
     641{
     642    vec_t<int, N> result;
     643
     644    for (int i = 0 ; i < N ; ++i)
     645        result[P[i]] = i;
     646
     647    return result;
     648}
     649
    629650/*
    630651 * Compute the determinant of a permutation square matrix corresponding to the permutation vector
     
    707728mat_t<T, N, N> inverse(mat_t<T, N, N> const &m)
    708729{
    709 #if 0 /* XXX: LU decomposition is currently broken */
    710730    mat_t<T, N, N> L, U;
    711     lu_decomposition(m, L, U);
    712 
    713     mat_t<T, N, N> ret = u_inverse(U) * l_inverse(L);
     731    vec_t<int, N> P = p_vector(m);
     732
     733    lu_decomposition(permute_rows(m, P), L, U);
     734
     735    mat_t<T, N, N> ret = permute_cols(u_inverse(U) * l_inverse(L), p_transpose(P));
    714736
    715737    return ret;
    716 #else
    717     mat_t<T, N, N> ret;
    718 
    719     T d = determinant(m);
    720     if (d)
    721     {
    722         /* Transpose result matrix on the fly */
    723         for (int i = 0; i < N; ++i)
    724             for (int j = 0; j < N; ++j)
    725                 ret[i][j] = cofactor(m, j, i) / d;
    726     }
    727     return ret;
    728 #endif
    729738}
    730739
  • trunk/src/t/math/matrix.cpp

    r3855 r3856  
    137137    }
    138138
    139 #if 0 /* XXX: LU decomposition is currently broken */
    140139    lolunit_declare_test(lu_decomposition_3x3)
    141140    {
     
    204203               vec4(0,  0, -1,  1));
    205204        mat4 L, U;
    206         lu_decomposition(m, L, U);
    207         mat4 m2 = L * U;
     205        vec_t<int, 4> P = p_vector(m);
     206        mat4 permuted = permute_rows(m, P);
     207
     208        lu_decomposition(permute_rows(m, P), L, U);
     209        mat4 m2 = permute_rows(L * U, p_transpose(P));
    208210
    209211        for (int j = 0; j < 4; ++j)
     
    288290            lolunit_assert_doubles_equal(m2[i][j], mat4(1.f)[i][j], 1e-5);
    289291    }
    290 #endif
    291292
    292293    lolunit_declare_test(inverse_3x3)
Note: See TracChangeset for help on using the changeset viewer.