CODEWITHSUNDEEP
fortranlapack
Alex
Alex

Reputation: 63

Question regarding the syntax of SGESV in Fortran

I am a bit confused on this subroutine. I have read the documentation but I am a bit confused what exactly the IPIV vector does and how exactly I set my leading dimension. I read that the leading dimension helps to find the starting point for the matrix elements in each successive column of the array. For example lets say we want to solve

Ax = B

where

integer, parameter :: sp = selected_real_kind(6,37)
real(kind=sp),dimension(:,:),intent(inout) :: A
real(kind=sp),dimension(:),intent(inout) :: B
integer, dimension(10) :: IPIV

where sp is for single precision which I have set in my main program and the dimensions are

A(10,10)
B(10)

which are set in my main program and passed to this subroutine Should I set my subroutine as

integer :: n,INFO
n = size(A,1)


IPIV = 0
call SGESV(n,n,A,2*n,IPIV,B,2*n,INFO)

or

call SGESV(n,n,A,n,IPIV,B,n,INFO)

and for IPIV I should just create a vector of size 10 and initialize it with zeros?

edit : I have used

call sgesv(n, n, A, n, ipiv, B, n, INFO)

as proposed as well but I get a segmentation error Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

I have printed the matrix sizes and they are correct which are the size of the matrix A is 100 and the size of the vector is 10

Edit2 : So in my main I have a loop which inside my loop it calculates a matrix of A (10,10) and a vector B(10) at each iteration. Then I call my subroutine to solve the system

call solver(A,B)

However I get the segmentation error which I do not understand since the dimensions are correct. (To check it I printed the size of the matrix and the vector and commented out the call to my subroutine and they are 100 and 10)

Perhaps I should make my matrices allocatable? But I do not see a problem with that since at each iteration I calculate the matrix and the vector though a series of calculations and overwrite them.

Basically I declare the matrix and the vector as follows

   real(sp) , dimension (10) :: B
   real(sp) , dimension (10,10) :: A

then inside my loop a series of calculations are performed to fill them with values and then I call my subroutine and then repeat with new values

Upvotes: 1

Views: 417

Answers (2)

jalex
jalex

Reputation: 1524

Below is an example of calling gesv the generic Lapack95 equivalent (and much simpler) of sgesv and dgesv.

subroutine test_lapack95(n)
use BLAS95
use LAPACK95
use f95_precision
implicit none

integer, intent(in) :: n
real(float), allocatable :: A(:,:), LU(:,:)
real(float), allocatable :: b(:), x(:)
integer, allocatable :: ipiv(:)

allocate(A(n,n))
allocate(b(n))
allocate(ipiv(n))

! Fill values in A and b
call prepare_values(n, A, b)
        
LU = A
x = b
call gesv(LU,x,ipiv)
! Solve to A*x=b, for x

end subroutine

don't worry about the helper function prepare_values, it just fill in A and b.

Upvotes: 0

jack
jack

Reputation: 1790

You are using an old interface to lapack. Note my lower answer for the modern/generic routine.

Old interface You would call it like

call sgesv(n, n, A, n, ipiv, B, n, info)

Reasoning:

  • leading dimensions are n and not 2n
  • ipiv is an output variable s.t. you dont need to initialize it with 0

Modern interface: LAPACK95 It is alot easier to just use the modern interfaces which provide generic calls as such

call gesv(A, B, ipiv=ipiv, info=info)

You dont need to specify the data types (e.g. no more sgesv) nor matrix dimensions.

Make sure that you need to use the appropriate module

use lapack95

Upvotes: 2

Related Questions