2

我正在尝试编写 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
4

1 回答 1

0

我无法使用 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 -g -check noarg_temp_created -implicitnone -fpe0
于 2013-09-03T18:35:00.847 回答