2015-08-14 49 views
0

我試圖使用FORTRAN 95.我求解線性方程系統的運行與LAPACK庫一個簡單的程序:Ax=B分割錯誤例程

A = [4 -2 3] 
    [1 3 -4] 
    [3 1 2] 

B=[ 1 
    -7 
    5] 

x被解向量

解是

x = [-1 
    2 
    3] 

這是代碼。我使用兩個子例程:SGETRFSGETRS。第一個函數SGETRF計算矩陣的LU分解,第二個子程序求解方程組。

program main 
implicit none 

integer :: i,j,info1,info2 
integer :: neqn ! number of equations 
real,dimension(3,3) :: coeff 
real,dimension (3,1) :: lhs 
real,dimension (3,1) :: ipiv 

neqn=3 

coeff = reshape((/4,1,3,-2,3,1,3,-4,2/),(/3,3/)) 
lhs = reshape ((/1,-7,5/),(/3,1/)) 

call SGETRF (neqn,1,coeff,neqn,ipiv,infO1) 
     if (info1==0) then 
      call SGETRS ('N',neqn,1,coeff,neqn,ipiv,lhs,neqn,info2) !Error 
     else 
     end if 

write (*,*) 'Answer: ' 
     do j=1,neqn,1 
      write (*,100) lhs(j,1) 
      100 format (F12.5,' ,') 
     end do 

     write (*,100) (lhs) 

end program 

按LAPACK文檔SGETRF,於我而言,M=neqn=3, N=1, A=coeff, LDA=3 我編譯的程序作爲gfortran main.f95 -o main -fcheck=all -llapack 而我得到的錯誤:

Program received signal SIGSEGV: Segmentation fault - invalid memory reference. 

Backtrace for this error: 
#0 0x7F758C3B3777 
#1 0x7F758C3B3D7E 
#2 0x7F758C00BD3F 
#3 0x7F758CA2F3EF 
#4 0x7F758C9BE8ED 
#5 0x400AE0 in MAIN__ at main.f95:19 
Segmentation fault (core dumped) 

線19 call SGETRS ('N',neqn,1,coeff,neqn,ipiv,lhs,neqn,info2) 我不明白爲什麼會出現是錯誤的。任何想法或意見?

回答

1

您的錯誤是由第二個參數SGETRF引起的。該參數是coeff的第二維,因此應爲3neqn

0

爲了詳細說明Stefan的正確答案,這裏是對代碼的修改。我相信,錯誤一些可能是由對LAPACK規格小心刪除編程(例如,ipiv陣列應該是等級1 Integer),並通過避免這麼多字面常量:

Program main 
    Implicit None 
    Integer, Parameter :: neqn = 3, nrhs = 1 
    Integer :: info 
    Real :: coeff(neqn, neqn) 
    Real :: lhs(neqn, nrhs) 
    Integer :: ipiv(neqn) 
    coeff = reshape([4,1,3,-2,3,1,3,-4,2], shape(coeff)) 
    lhs = reshape([1,-7,5], shape(lhs)) 
    Call sgetrf(neqn, size(coeff,2), coeff, size(coeff,1), ipiv, info) 
    If (info==0) Then 
    Call sgetrs('N', neqn, size(lhs,2), coeff, size(coeff,1), ipiv, lhs, & 
     size(lhs,1), info) 
    If (info==0) Then 
     Write (*, *) 'Answer: ' 
     Write (*, *)(lhs) 
    End If 
    End If 
End Program