Ich schreibe eine Bibliothek, wo ich einige grundlegende NxN-Matrix-Funktionalität haben möchte, die keine Abhängigkeiten hat und es ist ein bisschen ein Lernprojekt. Ich vergleiche meine Leistung mit Eigen. Ich konnte in der Lage sein, ziemlich gleich zu sein und sogar seine Leistung auf einer Paarfront mit SSE2 zu übertreffen und mit AVX2 an einigen Fronten zu schlagen (es benutzt nur SSE2, also nicht super überraschend).Schnellster Weg Determinante zu berechnen?
Mein Problem ist, ich Gaußsche Eliminations bin mit einem Ober diagonalisierte Matrix dann die Diagonale Multiplikation erstellen die determinant.I schlagen Eigen für N < 300 zu bekommen, aber danach bläst Eigen mich weg und es wird nur noch schlimmer als die Matrizen größer werden. Angesichts der Tatsache, dass auf den gesamten Speicher sequenziell zugegriffen wird und die Zerlegung des Compilers nicht schrecklich aussieht, glaube ich nicht, dass es sich um ein Optimierungsproblem handelt. Es gibt mehr Optimierung, die durchgeführt werden kann, aber die Timings sehen viel mehr wie ein algorithmisches Timing-Komplexitäts-Problem aus oder es gibt einen großen SSE-Vorteil, den ich nicht sehe. Einfach die Loops ein wenig abzurollen hat mir nicht viel gebracht.
Gibt es einen besseren Algorithmus zur Berechnung von Determinanten?
Scalar Code
/*
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;
}
Vorausgesetzt, Sie sind klar bereit, in Details zu graben, und Eigen ist Open Source ... [Warum nicht schauen, was es tut] (https://github.com/RLovelett/eigen/search?utf8 =% E2% 9C% 93 & q = Determinante) oder durchschreiten ...? – HostileFork
Ihre Methode macht für mich keinen Sinn und macht es schwierig, mich an die Arbeitsweise meiner Bibliothek anzupassen. Wenn ich die mathematische Argumentation dahinter verstehen würde, wäre ich in der Lage, sie einfach genug anzupassen. Ich denke, es hat etwas mit teilweisem Schwenken zu tun, in das ich zu schauen beginne. Die anderen Methoden haben für mich Sinn ergeben, aber dies ist das erste Mal, dass ich die Methodik dahinter nicht verstehen konnte. Im Allgemeinen, nur neugierig, ob ein anderes Gehirn eine Vorstellung von "dem besten Weg" hat. Auch nachdem ich diese Frage gepostet habe, schaue ich immer noch sehr genau hin und poste meinen Code, wenn ich einen besseren Weg finde. – user2927848
[dies] (http://stackoverflow.com/questions/22479258/cholesky-decomposition-with-openmp) könnte für Sie interessant sein. –