dgtsv - LAPACK not returning answer

Question

I'm trying to solve a simple tridiagonal system of equations using LAPACK library. The code below explains it all.

I'm getting an array full of zeros (initialized ones), not the correct answer.

I checked the inputs, tried to compile with two compilers and everything seems fine. What is wrong?

The compilation line is:

ifort -L/usr/local/lib/ -llapack -lblas tLapack.f90 -o tlapack  
gfortran -L/usr/local/lib/ -llapack -lblas tLapack.f90 -o tlapack

the code is:

program lapackT

  implicit none

  ! dgtsv( integer(4) :: N,
  !        integer(4) :: NRHS,
  !        real(8)    :: DL[],
  !        real(8)    :: D [],
  !        real(8)    :: DU[],
  !        real(8)    :: B [],
  !        integer(4) :: LDB ,
  !        integer(4) :: info )

  ! [A][x] = [b]
  ! N    - The order of matrix [A] 
  ! NRHS - Number of coluns in [b]
  ! DL   - Array with the subdiag.
  ! D    - Main diagonal.
  ! DU   - Upper Diagonal.
  ! B    - Answer !!
  ! LDB  - length of array [B].
  ! INFO - If = 0 .. Uhul !!.

  real(8), dimension(3) :: mainDiag
  real(8), dimension(2) :: lowerDiag
  real(8), dimension(2) :: upperDiag
  real(8), dimension(3) :: unknow
  real(8), dimension(3) :: equalty

  integer(4) :: info = 0
  integer(4) :: i = 0

  integer(4) :: N    = 3
  integer(4) :: NRHS = 1
  integer(4) :: LDB  = 3

  mainDiag(1) = 2.0d0
  mainDiag(2) = 2.0d0
  mainDiag(3) = 2.0d0

  upperDiag(1) = 3.0d0
  upperDiag(2) = 3.0d0

  lowerDiag(1) = 1.0d0
  lowerDiag(2) = 1.0d0

  equalty(1) = 1.0d0
  equalty(2) = 1.0d0
  equalty(3) = 1.0d0

  unknow = 0.0d0 ! answer

  call dgtsv(N,NRHS,lowerDiag,mainDiag,upperDiag,equalty,LDB,info)



  write(*,*) info

  do i = 1,size(unknow)
    write(*,*) unknow(i)
  end do

  ! Correct answer: unknow = (/-1,1,0/)    ! real(8) values
  ! Answer Im getting: unknow = (/0,0,0/)  ! real(8) values


end program lapackT

score 3 · Answer 1 · edited May 23 '17 at 10:33

If you look at the documentation you'll see the answer is returned in your equalty argument (i.e. it overwrites the RHS) - as unknow is not passed how can it be affected by the call? And I agree, I'm bot convinced this is the greatest design ever ...

While I'm here please learn about kinds, some of what you are doing is so quarter of a century ago. Please look at Fortran 90 kind parameter . Anyway here's how I would write your program (which some would say it also slightly out of date nowadays), and the answer it gives:

ian-admin@agon ~/test/stack $ cat dt.f90
program lapackT

  implicit none

  Integer, Parameter :: wp = Selected_real_kind( 13, 70 )

  real(wp), dimension(3) :: mainDiag
  real(wp), dimension(2) :: lowerDiag
  real(wp), dimension(2) :: upperDiag
  real(wp), dimension(3) :: unknow
  real(wp), dimension(3) :: equalty

  integer :: info = 0
  integer :: i = 0

  integer :: N    = 3
  integer :: NRHS = 1
  integer :: LDB  = 3

  mainDiag(1) = 2.0_wp
  mainDiag(2) = 2.0_wp
  mainDiag(3) = 2.0_wp

  upperDiag(1) = 3.0_wp
  upperDiag(2) = 3.0_wp

  lowerDiag(1) = 1.0_wp
  lowerDiag(2) = 1.0_wp

  equalty(1) = 1.0_wp
  equalty(2) = 1.0_wp
  equalty(3) = 1.0_wp

  unknow = 0.0_wp ! answer

  call dgtsv(N,NRHS,lowerDiag,mainDiag,upperDiag,equalty,LDB,info)



  write(*,*) info

  do i = 1,size(unknow)
    write(*,*) equalty(i)
  end do

  ! Correct answer: unknow = (/-1,1,0/)    ! real(8) values
  ! Answer Im getting: unknow = (/0,0,0/)  ! real(8) values


end program lapackT
ian-admin@agon ~/test/stack $ gfortran -Wall -Wextra -fcheck=all -O -std=f95 dt.f90  -llapack
ian-admin@agon ~/test/stack $ ./a.out
           0
  -1.0000000000000000     
   1.0000000000000000     
   0.0000000000000000E+000
ian-admin@agon ~/test/stack $

score 2 · Answer 2 · answered Apr 23 '16 at 08:41

Unless dgtsv operates through side-effects, this sequence of statements (your code, without empty lines):

  unknow = 0.0d0 ! answer
  call dgtsv(N,NRHS,lowerDiag,mainDiag,upperDiag,equalty,LDB,info)
  write(*,*) info
  do i = 1,size(unknow)
    write(*,*) unknow(i)
  end do

does not update unknow. How could it not be all 0.0s ?

Isn't the result returned through equalty in your call to dgtsv ?

dgtsv - LAPACK not returning answer

2 Answers2