Hinzufügen von Pivoting-Funktionalität zu einem Doolittle-Algorithmus

Question

Bisher habe ich diesen Code für eine LU-Dekomposition. Es nimmt ein Eingabearray auf und gibt die untere und obere Dreiecksmatrix zurück.

void LUFactorization (int d, const double*S, double*L, double*U) 
{ 
    for(int k = 0; k < d; ++k){ 
     if (
     for(int j = k; j < d; ++j){ 
     double sum = 0.; 
     for(int p = 0; p < k; ++p) { 
      sum+=L[k*d+p]*L[p*d+j]; 
      cout << L[k*d+p] << endl; 
     } 
     sum = S[k*d+j] - sum; 
     L[k*d+j]=sum;   
     U[k*d+j]=sum; 
     } 
     for(int i = k + 1; i < d; ++i){ 
     double sum=0.; 
     for(int p = 0; p < k; ++p) sum+=L[i*d+p]*L[p*d+k]; 
     L[i*d+k]=(S[i*d+k]-sum)/L[k*d+k]; 
     } 
    } 

    for(int k = 0; k < d; ++k){ 
     for(int j = k; j < d; ++j){ 
     if (k < j) L[k*d+j]=0; 
     else if (k == j) L[k*d+j]=1; 
     } 
    } 
} 

Gibt es eine Möglichkeit, dies anzupassen, um Reihenaustausch durchzuführen? Wenn nicht, gibt es einen anderen Algorithmus, auf den ich ausgerichtet sein könnte?

Dank

Answer 1

Der üblicher Ansatz für LU mit schwenkbarem Zersetzungen eine Permutation Array P zu speichern, die als die Identität Permutation initialisiert (P={0,1,2,...,d - 1}) und dann Swapping Einträge in P anstelle von Swapping Reihen in S

Wenn Sie dieses Permutationsarray haben, muss jeder Zugriff auf SP[i] anstelle von i für die Zeilennummer verwenden.

Beachten Sie, dass P selbst ist Teil der Ausgabe, da sie die Matrix Permutation darstellt, so dass
P*A = L*U, also wenn Sie es verwenden möchten, lineare Gleichungssysteme zu lösen, werden Sie P auf der rechten Seite beantragen Vor der Anwendung der Vorwärts- und Rückwärtssubstitution