0

I am trying to run a simple program with LAPACK library using fortran 95. I am solving linear system of equations as: Ax=B

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

B=[ 1
   -7
   5]

x is solution vector

Solution is

x = [-1
     2
     3]

Here is the code. I am using two subroutines: SGETRF and SGETRS. First function SGETRF computes LU decomposition of the matrix and second subroutine solves the system of equations.

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

As per LAPACK documentation SGETRF, in my case, M=neqn=3, N=1, A=coeff, LDA=3 I compiled the program as gfortran main.f95 -o main -fcheck=all -llapack And I get the error:

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)

Line 19 is call SGETRS ('N',neqn,1,coeff,neqn,ipiv,lhs,neqn,info2) I do not understand why there is error. Any idea or comments?

2 Answers 2

1

Your error is caused by the second parameter to SGETRF. This parameter is the second dimension of coeff and should thus be 3 or neqn.

Sign up to request clarification or add additional context in comments.

Comments

0

To elaborate on Stefan's correct answer, here is a reworking of your code. I believe some potential for error is removed by programming carefully against the LAPACK spec (for example, the ipiv array should be rank-1 Integer) and by avoiding so many literal constants:

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

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.