Как проверить диагональ матрицы? - PullRequest
Новая
       14

Как проверить диагональ матрицы?

7 голосов
/ 02 июня 2009

Мне нужно проверить, является ли одна матрица отклонений диагональной. Если нет, я сделаю разложение холесского ЛПНП. Но мне было интересно, какой самый надежный и быстрый способ проверить диагональ матрицы? Я использую Фортран.

Первое, что приходит мне в голову, это взять сумму всех элементов матрицы и вычесть диагональные элементы из этой суммы. Если ответ 0, матрица диагональна. Есть идеи получше?

В Фортране я напишу 

!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

Ответы [ 2 ]

10 голосов
/ 02 июня 2009

Было бы намного лучше просто обойти все недиагональные элементы и проверить, близки ли они к нулю (сравнение числа с плавающей точкой для неравенства склонно к ошибкам округления и может привести к ошибочным результатам).

Во-первых, как только вы найдете какой-либо нарушающий элемент, вы можете немедленно прекратить обход, и это может позволить значительно сократить время, если типичны нарушающие матрицы.

Во-вторых, это потенциально позволило бы улучшить развертывание цикла компилятором (компиляторы Fortran известны хорошими стратегиями оптимизации) и ускорить выполнение на кристалле благодаря меньшему количеству межинструкционных зависимостей.

Добавьте к этому тот факт, что предложенный вами алгоритм склонен к переполнению и накоплению ошибок, а алгоритм "traverse-and-test" - нет.

0 голосов
/ 16 июня 2009

Поиск в матрице ненулевых значений 

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

Чтобы использовать процедуры LAPACK двойной точности, замените первые «s» на «d». Так spotrf становится dpotrf

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

...