7

我需要测试一个方差矩阵是否是对角线。如果没有，我会做 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
4

2 回答 2

11

最好只遍历所有非对角线元素并测试它们是否接近零（比较浮点数的不等式容易出现舍入错误并可能导致错误结果）。

首先，一旦发现任何违规元素，您可以立即停止遍历，如果违规矩阵是典型的，这可能会显着减少时间。

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

除此之外，您建议的算法容易发生溢出和错误累积，而“遍历和测试”算法则不会。

于 2009-06-02T05:37:28.340 回答
0

在矩阵中搜索非零值

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/

于 2009-06-16T09:49:30.733 回答