Wie testen, ob die Matrix diagonal ist?

Question

Ich muss testen, ob eine Varianzmatrix diagonal ist. Wenn nicht, mache ich Cholesky LDL Zerlegung. Aber ich fragte mich, welcher der zuverlässigste und schnellste Weg zum Testen die Matrixdiagonale wäre? Ich benutze Fortran.

Das erste, was mir in den Sinn kommt, ist die Summe aller Elemente der Matrix und subtrahieren diagonale Elemente von dieser Summe. Wenn die Antwort 0 ist, ist die Matrix diagonal. Irgendwelche besseren Ideen?

In Fortran werde ich

!A is my matrix 
k=0.0d0 
do i in 1:n #n is the number of rows/colums 
k = k + A(i,i) 
end do 

if(abs(sum(A)-k) < epsilon(k)*sum(A)) then 
#do cholesky LDL, which I have to write myself, haven't found any subroutines for that in Lapack or anywhere else 
end if 

Answer 1

Suche die Matrix für die Nichtnullwerte

logical :: not_diag 
integer :: i, j 

not_diag = .false. 

outer: do i = 2, size(A,1) 
    do j = i, size(A, 2) 
    if (A(i,j) > PRECISION) then 
     not_diag = .true. 
     exit outer 
    end if 
    end 
end outer 

if (not_diag) then 
    ! DO LDL' decomposition 
end if 

Um die doppelte Genauigkeit LAPACK Routinen ersetzen die ersten ‚s‘ mit ‚d‘ zu verwenden. So spotrf wird dpotrf

http://www.netlib.org/lapack/double/

Answer 2

schreiben Es wäre viel besser, nur die außerhalb der Diagonalen alle Elemente und Tests zu durchlaufen, wenn sie in der Nähe von Null sind (eine Gleitkommazahl für Ungleichheit zu vergleichen ist anfällig für Rundungsfehler und kann zu fehlerhaften Ergebnissen führen).

Erstens, sobald Sie ein verletzendes Element finden, können Sie sofort mit dem Traversieren aufhören und dies kann eine signifikante Zeitabnahme ermöglichen, wenn die Verletzung von Matrizen typisch ist.

Zweitens, würde es möglicherweise für eine bessere Schleife entrolling durch den Compiler (Fortran Compiler sind bekannt für gute Optimierungsstrategien) und für eine schnellere On-Chip-Ausführung aufgrund weniger Inter-Befehl Abhängigkeiten.

Hinzu kommt die Tatsache, dass Ihr vorgeschlagener Algorithmus anfällig für Überläufe und Fehleranhäufung ist und der "traverse-and-test" -Algorithmus nicht ist.