Jak sprawdzić, czy macierz jest przekątna?

Question

Muszę sprawdzić, czy jedna macierz wariancji jest przekątna. Jeśli nie, zrobię Cholesky LDL dekompozycji. Ale zastanawiałem się, który byłby najbardziej wiarygodny i najszybszy sposób na sprawdzenie przekątnej macierzy? Używam Fortran.

Pierwszą rzeczą, która przychodzi mi do głowy, to zsumowanie wszystkich elementów macierzy i odjęcie elementów diagonalnych od tej sumy. Jeśli odpowiedź wynosi 0, macierz jest przekątna. Jakieś lepsze pomysły?

W Fortran będę pisać

!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

wyszukiwania matrycę niezerowej wartości

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 

Aby użyć podwójnej precyzji procedury Lapack zastąpienia pierwszego „S” z „d”. Więc spotrf staje dpotrf

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

Answer 2

Byłoby o wiele lepiej po prostu przechodzić wszystkich elementów niediagonalnych i testu, jeśli są bliskie zeru (porównując liczbę zmiennoprzecinkową na nierówność jest podatny na błędy zaokrąglania i mogą prowadzić do błędnych wyników).

Po pierwsze, po znalezieniu dowolnego elementu naruszającego można natychmiast zatrzymać przechodzenie i może to pozwolić na znaczny spadek czasu, jeśli typowe są naruszenia macierzy.

Po drugie, potencjalnie pozwoliłoby to na lepsze rozwijanie pętli przez kompilator (kompilatory Fortran znane są z dobrych strategii optymalizacji) i szybsze wykonywanie mikroukładów ze względu na mniej zależności między instrukcjami.

Dodaj do tego fakt, że sugerowany algorytm jest podatny na przepełnienia i kumulację błędów, a algorytm "przechodź i testuj" nie jest.