我正在寫一個庫,我想要有一些基本的NxN矩陣功能,沒有任何依賴關係,它是一個學習項目。我將我的表現與Eigen進行比較。我已經能夠相當平等,甚至可以在SSE2的前面擊敗它的表現,並且AVX2在很多方面都打敗了它(它只使用SSE2,所以不會超級意外)。計算行列式的最快方法?
我的問題是我使用高斯消去來創建一個上對角化矩陣,然後乘以對角線得到行列式。我擊敗Eigen爲N < 300,但之後,Eigen吹走了我,並作爲矩陣變得更糟變得更大。鑑於所有的內存是按順序訪問的,編譯器的反彙編看起來並不可怕,我不認爲這是一個優化問題。還有更多優化可以完成,但時序看起來更像是算法時序複雜性問題,或者我沒有看到主要的SSE優勢。簡單地展開這些循環對於我來說沒有太大的幫助。
是否有更好的算法來計算行列式?
標量代碼
/*
Warning: Creates Temporaries!
*/
template<typename T, int ROW, int COLUMN> MML_INLINE T matrix<T, ROW, COLUMN>::determinant(void) const
{
/*
This method assumes square matrix
*/
assert(row() == col());
/*
We need to create a temporary
*/
matrix<T, ROW, COLUMN> temp(*this);
/*We convert the temporary to upper triangular form*/
uint N = row();
T det = T(1);
for (uint c = 0; c < N; ++c)
{
det = det*temp(c,c);
for (uint r = c + 1; r < N; ++r)
{
T ratio = temp(r, c)/temp(c, c);
for (uint k = c; k < N; k++)
{
temp(r, k) = temp(r, k) - ratio * temp(c, k);
}
}
}
return det;
}
AVX2
template<> float matrix<float>::determinant(void) const
{
/*
This method assumes square matrix
*/
assert(row() == col());
/*
We need to create a temporary
*/
matrix<float> temp(*this);
/*We convert the temporary to upper triangular form*/
float det = 1.0f;
const uint N = row();
const uint Nm8 = N - 8;
const uint Nm4 = N - 4;
uint c = 0;
for (; c < Nm8; ++c)
{
det *= temp(c, c);
float8 Diagonal = _mm256_set1_ps(temp(c, c));
for (uint r = c + 1; r < N;++r)
{
float8 ratio1 = _mm256_div_ps(_mm256_set1_ps(temp(r,c)), Diagonal);
uint k = c + 1;
for (; k < Nm8; k += 8)
{
float8 ref = _mm256_loadu_ps(temp._v + c*N + k);
float8 r0 = _mm256_loadu_ps(temp._v + r*N + k);
_mm256_storeu_ps(temp._v + r*N + k, _mm256_fmsub_ps(ratio1, ref, r0));
}
/*We go Scalar for the last few elements to handle non-multiples of 8*/
for (; k < N; ++k)
{
_mm_store_ss(temp._v + index(r, k), _mm_sub_ss(_mm_set_ss(temp(r, k)), _mm_mul_ss(_mm256_castps256_ps128(ratio1),_mm_set_ss(temp(c, k)))));
}
}
}
for (; c < Nm4; ++c)
{
det *= temp(c, c);
float4 Diagonal = _mm_set1_ps(temp(c, c));
for (uint r = c + 1; r < N; ++r)
{
float4 ratio = _mm_div_ps(_mm_set1_ps(temp[r*N + c]), Diagonal);
uint k = c + 1;
for (; k < Nm4; k += 4)
{
float4 ref = _mm_loadu_ps(temp._v + c*N + k);
float4 r0 = _mm_loadu_ps(temp._v + r*N + k);
_mm_storeu_ps(temp._v + r*N + k, _mm_sub_ps(r0, _mm_mul_ps(ref, ratio)));
}
float fratio = _mm_cvtss_f32(ratio);
for (; k < N; ++k)
{
temp(r, k) = temp(r, k) - fratio*temp(c, k);
}
}
}
for (; c < N; ++c)
{
det *= temp(c, c);
float Diagonal = temp(c, c);
for (uint r = c + 1; r < N; ++r)
{
float ratio = temp[r*N + c]/Diagonal;
for (uint k = c+1; k < N;++k)
{
temp(r, k) = temp(r, k) - ratio*temp(c, k);
}
}
}
return det;
}
鑑於您已經準備好深入細節,並且Eigen是開源的... [爲什麼不看它的功能](https://github.com/RLovelett/eigen/search?utf8 =%E2%9C%93&q =行列式)或者通過它......? – HostileFork
他們的方法對我來說沒有多大意義,因此很難適應我圖書館的運作方式。如果我理解了它背後的數學推理,我就可以很容易地適應它。我認爲這與我正在開始研究的部分旋轉有關。其他方法對我來說是有意義的,但這是我第一次無法理解它背後的方法。一般來說,只是好奇另一個大腦是否有'最好的方式'的想法。即使在發佈這個問題之後,我仍然非常關注它,並會在找到更好的方法時發佈我的代碼。 – user2927848
[this](http://stackoverflow.com/questions/22479258/cholesky-decomposition-with-openmp)可能對你很有意思。 –