无法在三对角矩阵求解器中提取对角线作为向量答案

【问题标题】：Can't extract diagonals as vectors in a tridiagonal matrix solver无法在三对角矩阵求解器中提取对角线作为向量
【发布时间】：2019-07-20 10:48:01
【问题描述】：

我正在尝试编写 Thomas 算法来解决 AX=b 形式的三对角矩阵问题，其中 A 是三对角矩阵，X 是未知向量，b 是独立项向量。

当我尝试提取主对角线A(i,i)、左右对角线A(i,i+1) 和A(i,i-1) 时，程序无法这样做。真正尴尬的是，如果我将 print *, 放在 do 循环的中间，它会起作用。

我知道在单独的向量中获取对角线并不是绝对必要的，但为了解释的目的，我试图尽可能明确地做到这一点。

有人可以帮忙吗？

提前致谢

这是求解向量 X 的子程序（AM 是矩阵）

subroutine LU(N,AM,D,X)
implicit none
integer(kind=4) i,N
real(kind=4),dimension(N,N) :: AM
real(kind=4),dimension(N) :: G,X,D
real(kind=4),dimension(N) :: A,B,C,L,U

A(1)=AM(1,1)
B(1)=AM(1,2)
C(1)=0
A(N)=AM(N,N)
B(N)=0
C(N)=AM(N,N-1)

L(1)=A(1)
U(1)=B(1)/A(1)
G(1)=D(1)/L(1)

do i=2,N-1,+1
    C(i)=AM(i,i-1)
    A(i)=AM(i,i)
    print *,      !THIS IS THE "MAGICAL" print *,. REMOVE IT AND IT WON`T WORK
    B(i)=AM(i,i+1)

    L(i)=A(i)-C(i)*U(i-1)
    U(i)=B(i)/L(i)
    G(i)=(D(i)-C(i)*G(i-1))/L(i)
end do

i=N

L(i)=A(i)-C(i)*U(i-1)
U(i)=B(i)/L(i)
G(i)=(D(i)-C(i)*G(i-1))/L(i)

X(N)=G(N)

do i=N-1,1,-1
    X(i)=G(i)-U(i)*X(i+1)
end do

end subroutine LU

【问题讨论】：

你的意思是“程序没有这样做”？我将 -0.5 用于非对角线，将 +1 用于AM 的对角线，然后使用d=(/0.9, 1.1, 1.3, 1.4, 1.5/) 并收到5.7000003 9.6000004 11.300000 10.400000 6.6999998 作为带有和不带有print 语句的输出。我在 Ubuntu 12.04 上并使用 gfortran 4.6.3
如果插入一个不必要的print 语句“修复”了一个问题，那么代码就存在严重错误，可能访问数组边界之外的元素。我无法立即看到代码的数组索引错误，但我建议您重新编译并重新运行，并打开数组边界检查 - 请查看编译器的文档以了解如何执行此操作。
再次查看您的代码不应该是向量B 和C，我理解它们是非主对角线，大小为N-1？
@HighPerformanceMark 我认为他为此将C(1) 和B(N) 分配给0，因此实际上只使用了每个N-1 元素。

标签： fortran linear-algebra

【解决方案1】：

我无法使用 ifort 12.0 或 gfortran 4.8.1 重现该错误。您使用的是哪个编译器/版本？

我还对照 LAPACK 检查了您的例程，它按预期工作！

subroutine LU_lapack(A, b, x, stat)
  implicit none
  ! Arguments
  real(kind=4),intent(inout) :: A(:,:)
  real(kind=4),intent(in)    :: b(:,:)
  real(kind=4),intent(out)   :: x(:,:)
  integer,intent(out)             :: stat
  ! Local variables
  integer, allocatable            :: IPIV(:) ! Pivot indices for LAPACK
  integer                         :: N, NRHS

  x = b

  allocate( IPIV(size(x,1)),stat=stat )
  if (stat.ne.0) stop 'solve_LU_double: Cannot allocate memory!'

  N    = size(b,1)
  NRHS = size(b,2)

  ! Perform the LU-decomposition
  call SGETRF(N,N,A,N,IPIV,stat)
  if ( stat.ne.0 ) then
    deallocate(IPIV)
    return
  endif

  ! Solve system
  call SGETRS( 'N', N, NRHS, A, N, IPIV, x, N, stat )

  deallocate(IPIV)
end subroutine

根据你报告的奇怪行为，我猜在你调用你的例程之前有问题 - 可能是越界或错位的指针...... 你应该通过 valgrind 运行你的代码来发现这样的异常！

编辑：正如之前在其中一篇文章中提到的，使用编译器标志来检查数组边界等也可能有助于定位问题......首先尝试一下，它通常会快得多！对于 gfortran，我通常使用

gfortran -g -fimplicit-none -ffpe-trap=invalid,zero,overflow \
         -fbounds-check -Wall -Wextra

为了我使用的ifort

ifort -g -check noarg_temp_created -implicitnone -fpe0

【讨论】：