DGESV提供了错误的解决方案

Question

全部。我正在尝试通过调用拉帕克（Lapack）调用DGESV来解决Fortran中的AX = B系统。在哪里，

a = [

 1     0     1     0     1     0     1
 0     1     0     1     0     1     0
 0     0     4     0    32     0   108
 0     0     0    24     0   120     0
 0     0     0     0    48     0   192
 0     0     0     0     0    80     0
 0     0     0     0     0     0   120]

b = [ 0 0 2 6 0 0 0]

x可以简单地在matlab中计算为x = a b，它给出x =

-0.5000 -0.2500 0.5000 0.2500 0 0 0

在我在Fortran中做同样的事情时，它给了我完全不同的价值观。我的代码有什么问题，还是我在调用DGESV时犯了错误？这是我的fortran代码：

  program GG
  implicit none
  integer, parameter :: N=7
  integer :: i,j,ipiv(N),ok
  real(8), dimension(1:N,1:N) :: M
  real(8), dimension(1,1:N) :: b

  M(:,1)=(/1.,0.,1.,0.,1.,0.,1./)*1.d0
  M(:,2)=(/0.,1.,0.,1.,0.,1.,0./)*1.d0
  M(:,3)=(/0.,0.,2.**3/2,0.,4.**3/2,0.,6.**3/2/)*1.d0
  M(:,4)=(/0.,0.,0.,3.*(3.**2-1.),0.,5.*(5.**2-1.),0./)*1.d0
  M(:,5)=(/0.,0.,0.,0.,4.*(4.**2-2.**2),0.,6.*(6.**2-2.**2)/)*1.d0
  M(:,6)=(/0.,0.,0.,0.,0.,5.*(5.**2-3.**2),0./)*1.d0
  M(:,7)=(/0.,0.,0.,0.,0.,0.,6.*(6.**2-4.**2)/)*1.d0
  b=reshape((/0,0,2,6,0,0,0/)*1.d0,shape(b))

  call DGESV(N,1,M,N,ipiv,b,N,ok)
  write(*,*), b
  end program GG

此代码给出的结果是：

0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000来零0.2500000000000000000000000000000000000.33333333333333333337-0.37500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000来度

谢谢。

Answer 1

您在fortran中创建的矩阵似乎是您在MATLAB中定义的矩阵（行和列）的转置（行和列）。注意（在MATLAB中）：

>> A = [1     0     1     0     1     0     1
        0     1     0     1     0     1     0
        0     0     4     0    32     0   108
        0     0     0    24     0   120     0
        0     0     0     0    48     0   192
        0     0     0     0     0    80     0
        0     0     0     0     0     0   120];
>> B = [0 0 2 6 0 0 0];
>> x = AB.'
x =
   -0.5000
   -0.2500
    0.5000
    0.2500
         0
         0
         0
>> x = (A.')B.'  % A is transposed, and you get your Fortran result
x =
         0
         0
    0.5000
    0.2500
   -0.3333
   -0.3750
    0.0833

换句话说，在fortran中执行M(:,1)=...填充第一个列，而不是第一个 row 。如果将它们翻转到M(1,:)=...等等，我认为它应该与MATLAB的结果相匹配。