詳細講解矩陣求逆的快速演算法(轉)

詳細講解矩陣求逆的快速演算法(轉)[@more@]

　　演算法介紹

　　矩陣求逆在3D程式中很常見，主要應用於求Billboard矩陣。按照定義的計算方法乘法運算，嚴重影響了效能。在需要大量Billboard矩陣運算時，矩陣求逆的最佳化能極大提高效能。這裡要介紹的矩陣求逆演算法稱為全選主元高斯-約旦法。

　　

　　高斯-約旦法（全選主元）求逆的步驟如下：

　　

　　首先，對於 k 從 0 到 n - 1 作如下幾步：

　　

　　從第 k 行、第 k 列開始的右下角子陣中選取絕對值最大的元素，並記住次元素所在的行號和列號，在透過行交換和列交換將它交換到主元素位置上。這一步稱為全選主元。

　　m(k, k) = 1 / m(k, k)

　　m(k, j) = m(k, j) * m(k, k)，j = 0, 1, ..., n-1；j != k

　　m(i, j) = m(i, j) - m(i, k) * m(k, j)，i, j = 0, 1, ..., n-1；i, j != k

　　m(i, k) = -m(i, k) * m(k, k)，i = 0, 1, ..., n-1；i != k

　　最後，根據在全選主元過程中所記錄的行、列交換的資訊進行恢復，恢復的原則如下：在全選主元過程中，先交換的行（列）後進行恢復；原來的行（列）交換用列（行）交換來恢復。

　　

　　實現(4階矩陣)

　　float Inverse(CLAYMATRIX& mOut, const CLAYMATRIX& rhs)

　　{

　　CLAYMATRIX m(rhs);

　　DWORD is[4];

　　DWORD js[4];

　　float fDet = 1.0f;

　　int f = 1;

　　

　　for (int k = 0; k < 4; k ++)

　　{

　　// 第一步，全選主元

　　float fMax = 0.0f;

　　for (DWORD i = k; i < 4; i ++)

　　{

　　for (DWORD j = k; j < 4; j ++)

　　{

　　const float f = Abs(m(i, j));

　　if (f > fMax)

　　{

　　fMax = f;

　　is[k] = i;

　　js[k] = j;

　　}

　　}

　　}

　　if (Abs(fMax) < 0.0001f)

　　return 0;

　　

　　if (is[k] != k)

　　{

　　f = -f;

　　swap(m(k, 0), m(is[k], 0));

　　swap(m(k, 1), m(is[k], 1));

　　swap(m(k, 2), m(is[k], 2));

　　swap(m(k, 3), m(is[k], 3));

　　}

　　if (js[k] != k)

　　{

　　f = -f;

　　swap(m(0, k), m(0, js[k]));

　　swap(m(1, k), m(1, js[k]));

　　swap(m(2, k), m(2, js[k]));

　　swap(m(3, k), m(3, js[k]));

　　}

　　

　　// 計算行列值

　　fDet *= m(k, k);

　　

　　// 計算逆矩陣

　　

　　// 第二步

　　m(k, k) = 1.0f / m(k, k);

　　// 第三步

　　for (DWORD j = 0; j < 4; j ++)

　　{

　　if (j != k)

　　m(k, j) *= m(k, k);

　　}

　　// 第四步

　　for (DWORD i = 0; i < 4; i ++)

　　{

　　if (i != k)

　　{

　　for (j = 0; j < 4; j ++)

　　{

　　if (j != k)

　　m(i, j) = m(i, j) - m(i, k) * m(k, j);

　　}

　　}

　　}

　　// 第五步

　　for (i = 0; i < 4; i ++)

　　{

　　if (i != k)

　　m(i, k) *= -m(k, k);

　　}

　　}

　　

　　for (k = 3; k >= 0; k --)

　　{

　　if (js[k] != k)

　　{

　　swap(m(k, 0), m(js[k], 0));

　　swap(m(k, 1), m(js[k], 1));

　　swap(m(k, 2), m(js[k], 2));

　　swap(m(k, 3), m(js[k], 3));

　　}

　　if (is[k] != k)

　　{

　　swap(m(0, k), m(0, is[k]));

　　swap(m(1, k), m(1, is[k]));

　　swap(m(2, k), m(2, is[k]));

　　swap(m(3, k), m(3, is[k]));

　　}

　　}

　　

　　mOut = m;

　　return fDet * f;

　　}

　　

　　比較

　　原演算法　 原演算法　　(經過高度最佳化)　　 新演算法

　　加法次數　 103　　　　 61　　　　　　　 39

　　乘法次數　 170　　　　 116　　　　　　　69

　　需要額外空間 16 * sizeof(float) 34 * sizeof(float) 25 * sizeof(float)

　　結果不言而喻