Comment tester si la matrice est diagonale?

Question

J'ai besoin de tester si une matrice de variance est diagonale. Sinon, je vais faire une décomposition de Cholesky LDL. Mais je me demandais quel serait le moyen le plus fiable et le plus rapide de tester est la diagonale de la matrice? J'utilise Fortran.

La première chose qui me vient à l'esprit est de prendre la somme de tous les éléments de la matrice, et de soustraire les éléments diagonaux de cette somme. Si la réponse est 0, la matrice est diagonale. De meilleures idées?

Dans Fortran Je vais écrire

!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

Rechercher la matrice pour les valeurs non nul

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 

Pour utiliser la double précision routines LAPACK remplacer le premier 'l' avec 'd'. Alors spotrf devient dpotrf

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

Answer 2

Il serait préférable de traverser juste tous les éléments hors diagonale et test si elles sont proches de zéro (la comparaison d'un nombre à virgule flottante pour l'inégalité est sujette à erreurs d'arrondi et peut conduire à des résultats erronés). Tout d'abord, une fois que vous avez trouvé un élément de violation, vous pouvez arrêter immédiatement la traversée, ce qui peut permettre une diminution significative du temps si les matrices de violation sont typiques. Deuxièmement, cela permettrait potentiellement un meilleur déroulement de la boucle par le compilateur (les compilateurs Fortran sont connus pour de bonnes stratégies d'optimisation) et pour une exécution sur puce plus rapide en raison de moins de dépendances inter-instructions. Ajoutez à cela le fait que votre algorithme suggéré est sujet aux débordements et à l'accumulation d'erreurs et que l'algorithme "traverse-et-test" ne l'est pas.