如何测试矩阵是否对角线？

Question

我需要测试一个方差矩阵是否是对角线.如果没有,我会做Cholesky LDL分解.但我想知道哪种最可靠,最快的测试方法是矩阵对角线？我正在使用Fortran.

我想到的第一件事是获取矩阵的所有元素的总和,并从该总和中减去对角线元素.如果答案是0,则矩阵是对角线的.有更好的想法吗？

在Fortran我会写

!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

只是遍历所有非对角线元素并测试它们是否接近于零(比较不等式的浮点数易于舍入误差并且可能导致错误结果)会好得多.

首先,一旦找到任何违规元素,您可以立即停止遍历,如果违反矩阵是典型的,这可能会显着缩短时间.

其次,它可能允许编译器更好地循环展开(Fortran编译器以良好的优化策略而闻名),并且由于较少的指令间依赖性,可以实现更快的片上执行.

除此之外,您建议的算法容易出现溢出和错误累积,而"遍历和测试"算法则不然.

Answer 2

在矩阵中搜索非零值

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/