DGBSV: incorrect numerical results with -L/usr/lib64/atlas ?

DGBSV: incorrect numerical results with -L/usr/lib64/atlas ?

[Jorge D'Elia]

Hi all,

I do not know if this forum is the most appropriate, but as it based 
using gfortran I will try to start here. 

I usually build the ATLAS library in some beta version without any 
numerical problem that I remember. In a circumstantial way, I have to 
use the system ATLAS library and, today, I am getting some incorrect 
numerical results with the DGBSV lapack fortran subroutine, which is 
described below with a simple test. 

Perhaps, I am not using the correct sintaxis in the interface, maybe 
in the compiler flags or the link ones? Please, any ideas? 
Thanks in advance.

Regards.
Jorge.
--
A simple test: it is built for the General Band Matrix Factorization
and Multiple Right-Hand Side Solver DGBSV lapack fortran subroutine,
and it is based from:
  IBM Documentation
  SGBSV, DGBSV, CGBSV, and ZGBSV
  (General Band Matrix Factorization and Multiple Right-Hand Side Solve)
  https://www.ibm.com/docs/en/essl/6.2?topic=blaes-sgbsv-dgbsv-cgbsv-zgbsv-general-band-matrix-factorization-multiple-right-hand-side-solve

$ cat /proc/version
Linux version 5.13.19-200.fc34.x86_64 (mockbuild@bkernel01.iad2.fedoraproject.org) (gcc (GCC) 11.2.1 20210728 (Red Hat 11.2.1-1), GNU ld version 2.35.2-5.fc34) #1 SMP Sat Sep 18 16:32:24 UTC 2021

$ sudo dnf --best reinstall atlas*

$ which gfortran
/usr/beta/gcc-trunk/bin/gfortran

$ gfortran --version
GNU Fortran (GCC) 12.0.0 20210914 (experimental)

$ gfortran -march=native -mtune=native -fall-intrinsics -fcheck=all -fimplicit-none -fmax-errors=4 -fPIC -std=f2018 -Wall -Waliasing -Warray-temporaries -Wcharacter-truncation -Werror -Wextra -Wimplicit-interface -Wimplicit-procedure -Wintrinsic-shadow -Wline-truncation -Wuse-without-only -Wrealloc-lhs-all -Wsurprising -Wtabs -Wuninitialized -Wunused-parameter -Og -o trouble.exe trouble.f90  -L/usr/lib64/atlas -llapack -lf77blas -lcblas -ltatlas

(Note: the same results either -ltatlas, -lsatlas, -flto=auto flags).

$ ./trouble.exe

(snip)

info =  1   # <-- i.e. the solution has not been computed.

$ cat trouble.f90

module m_interfase
  implicit none
  private
  public :: idp, iin, dgbsv
  integer, parameter :: idp = kind (1.0d0)
  integer, parameter :: iin = kind (1)    
  interface
    subroutine dgbsv (n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info)
      integer           info, kl, ku, ldab, ldb, n, nrhs
      integer           ipiv (*)
      double precision  ab (ldab,*), b (ldb,*)
      external          dgbtrf, dgbtrs, xerbla
      intrinsic         max
    end subroutine dgbsv
  end interface
end module m_interfase
!
program trouble
  use  m_interfase, only: iin, idp, dgbsv
  implicit none
  real (idp), parameter :: cero = 0.0_idp
  real (idp), parameter :: uno  = 1.0_idp
  real (idp), parameter :: dos  = 2.0_idp
  real (idp), parameter :: tres = 3.0_idp
  real (idp), parameter :: cuat = 4.0_idp
  real (idp), parameter :: cinc = 5.0_idp
  real (idp), parameter :: seis = 6.0_idp
  real (idp), parameter :: siet = 7.0_idp
  real (idp), parameter :: ocho = 8.0_idp
  real (idp), parameter :: nuev = 9.0_idp
  real (idp), dimension (9,9), parameter :: &
    ccc = reshape ( [ &
      uno,  dos, tres, cero, cero, cero, cero, cero, cero, & ! 1
      dos,  uno,  dos, tres, cero, cero, cero, cero, cero, & ! 2
     tres,  dos,  uno,  dos, tres, cero, cero, cero, cero, & ! 3
     cuat, tres,  dos,  uno,  dos, tres, cero, cero, cero, & ! 4
     cero, cuat, tres,  dos,  uno,  dos, tres, cero, cero, & ! 5
     cero, cero, cuat, tres,  dos,  uno,  dos, tres, cero, & ! 6
     cero, cero, cero, cuat, tres,  dos,  uno,  dos, tres, & ! 7
     cero, cero, cero, cero, cuat, tres,  dos,  uno,  dos, & ! 8
     cero, cero, cero, cero, cero, cuat, tres,  dos,  uno  & ! 9
     ], [9,9])
  real (idp), dimension (9,3), parameter :: &
    vvv = reshape ( [ &
      uno,  uno,  uno,  uno,  uno,  uno,  uno,  uno,  uno, & ! 1
      uno, -uno,  uno, -uno,  uno, -uno,  uno, -uno,  uno, & ! 2
      uno,  dos, tres, cuat, cinc, seis, siet, ocho, nuev  & ! 3
      ], [9,3])
  real    (idp), dimension (9,9) :: aaa
  real    (idp), dimension (9,3) :: bbb
  integer (iin), dimension (9)   :: ipiv
  integer (iin) :: nn, klo, kup, nrhs, lda, ldb, info
  integer (iin) :: i
  !
  nn   = size (aaa,1)
  klo  = 2
  kup  = 3
  nrhs = size (bbb,2)
  lda  = 2 * klo + kup + 1
  ldb  = nn
  aaa  = ccc
  bbb  = vvv
  ipiv = 0
  !
  write (*,*)
  write (*,*)" dgbsv (General Band Side) test "
  write (*,*)
  write (*,*)" system matrix: aaa  "
  do  i = 1, nn
  write (*,*) real (aaa (i,:))
  end do
  write (*,*)
  write (*,*)" source: bbb  "
  do  i = 1, nn
  write (*,*) real (bbb (i,:))
  end do
  write (*,*)" nn   = ", nn
  write (*,*)" klo  = ", klo
  write (*,*)" kup  = ", kup
  write (*,*)" nrhs = ", nrhs
  write (*,*)" lda  = ", lda
  write (*,*)" ldb  = ", ldb
  !           N  KL   KU   NRHS   A   LDA   IPIV  B    LDB  INFO)
  !           |   |    |    |     |    |    |     |      |     |
 !call dgbsv (9 , 2  , 3  , 3   , aaa, 8  , ipiv, bbb,   9, info)
  call dgbsv (nn, klo, kup, nrhs, aaa, lda, ipiv, bbb, ldb, info)
  !
  write (*,*)
  write (*,*)" output aaa: "
  do  i = 1, nn
  write (*,*) real (aaa (i,:))
  end do
  write (*,*)
  write (*,*)" output bbb: "
  do  i = 1, nn
  write (*,*) real (bbb (i,:))
  end do 
  write (*,*)
  write (*,*)" output ipiv: ", ipiv
  write (*,*)
  write (*,*)" info = ", info
end program trouble
--

Re: DGBSV: incorrect numerical results with -L/usr/lib64/atlas ?

[Thomas Koenig]

Hi Jorge,

> I do not know if this forum is the most appropriate, but as it based
> using gfortran I will try to start here.
> 
> I usually build the ATLAS library in some beta version without any
> numerical problem that I remember. In a circumstantial way, I have to
> use the system ATLAS library and, today, I am getting some incorrect
> numerical results with the DGBSV lapack fortran subroutine, which is
> described below with a simple test.
> 
> Perhaps, I am not using the correct sintaxis in the interface, maybe
> in the compiler flags or the link ones? Please, any ideas?
> Thanks in advance.

I tried your code against the reference LAPACK and got the result below,
which also has info=1 at the end.  There is no argument mismatch
(at least) in the called routine, I checked that by including the
refblas routine into the source.

So, I don't know what's wrong.  Possibly something with the numerics,
or a wrong size of an argument?  It would be possible to compile the
reference LAPACK and BLAS with debugging and check with a debugger
where the info flag is set.

If I may guess, it's probably an error in one of the array sizes -
LAPACK is notorious that way.  I wish they would use modern Fortran's
array passing mechanism.

Best regards

	Thomas

Re: DGBSV: incorrect numerical results with -L/usr/lib64/atlas ?

[Jorge D'Elia]

Hi Thomas,

----- Mensaje original -----
> De: "Thomas Koenig" 
> Para: "Jorge D'Elia" , "Gfortran List" <fortran@gcc.gnu.org>
> CC: "Jorge D'Elia" 
> Enviado: Miércoles, 29 de Septiembre 2021 3:02:46
> Asunto: Re: DGBSV: incorrect numerical results with -L/usr/lib64/atlas ?
>
> Hi Jorge,
> 
>> I do not know if this forum is the most appropriate, but as it based
>> using gfortran I will try to start here.
>> 
>> I usually build the ATLAS library in some beta version without any
>> numerical problem that I remember. In a circumstantial way, I have to
>> use the system ATLAS library and, today, I am getting some incorrect
>> numerical results with the DGBSV lapack fortran subroutine, which is
>> described below with a simple test.
>> 
>> Perhaps, I am not using the correct sintaxis in the interface, maybe
>> in the compiler flags or the link ones? Please, any ideas?
>> Thanks in advance.
> 
> I tried your code against the reference LAPACK and got the result below,
> which also has info=1 at the end.  There is no argument mismatch
> (at least) in the called routine, I checked that by including the
> refblas routine into the source.

Ok. Thanks for those checks.
 
> So, I don't know what's wrong.  Possibly something with the numerics,
> or a wrong size of an argument?  It would be possible to compile the
> reference LAPACK and BLAS with debugging and check with a debugger
> where the info flag is set.
> 
> If I may guess, it's probably an error in one of the array sizes -

Good catch! (see below).

> LAPACK is notorious that way.  I wish they would use modern 
> Fortran's array passing mechanism.

I agree. That update would be very welcome looking to the future.

> Best regards
> 
>         Thomas

Well. After a couple of coffees, it turns out that it 
was a false alarm because the end user wrongly used 
the lapack' subroutines. 

As often happens it was an oversight: in the test I forgot
to insert a call to some routine that converts the usual
representation of an array to the BLAS-General-Band Storage
format. Below I copy the corrected test and everything 
returned to normal. 

Sorry for the noise. Thus, atlas and gfortran work 
very well, as long as they are used correctly ...

Many thanks for your time in addressing this issue.


Best regards. 
Jorge.

!begin (fixed) test
! BLAS-general-band storage mode:
!   the storage mode packs the matrix elements by columns into
!   a two-dimensional array, such that each diagonal of the
!   matrix appears as a row in the packed array.
module m_interfase
  implicit none
  private
  public :: idp, iin, dgbsv, conversor_agb
  integer, parameter :: idp = kind (1.0d0)
  integer, parameter :: iin = kind (1)    
  !
  interface
    subroutine dgbsv (n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info)
      integer           info, kl, ku, ldab, ldb, n, nrhs
      integer           ipiv (*)
      double precision  ab (ldab,*), b (ldb,*)
      external          dgbtrf, dgbtrs, xerbla
      intrinsic         max
    end subroutine dgbsv
  end interface
  !
contains
  !
  ! ................................................................
  ! A conversor from the general band matrix AAA (mm,nn) to
  ! the BLAS-general-band storage mode AGB (lda,nn) for dgbsv,
  ! where AGB has klo extra rows (i.e. lda = (ku + 1 + klo) + klo)
  ! needed for the dgbsv solution. Ref.
  !   https://www.ibm.com/docs/en/essl/
  !   /6.2?topic=representation-blas-general-band-storage-mode
  ! ................................................................
  subroutine conversor_agb (mm,nn,lda,kup,klo,aaa,agb)
    integer (iin), intent (in)  :: mm, nn, lda, kup, klo
    real    (idp), intent (in)  :: aaa (mm,nn)
    real    (idp), intent (out) :: agb (lda,nn)
    integer (iin) :: i, j, k, p
    agb = 0.0_idp
    do  j = 1, nn
      k = kup + 1 - j
      do  i = max (1,j-kup), min (mm,j+klo)
        p = klo + (k + i)
        agb (p,j) = aaa (i,j)
      end do
    end do
  end subroutine conversor_agb
  !
end module m_interfase
!
program extrouble
  use  m_interfase, only: iin, idp, dgbsv, conversor_agb
  implicit none
  !
  real (idp), parameter :: cero = 0.0_idp
  real (idp), parameter :: uno  = 1.0_idp
  real (idp), parameter :: dos  = 2.0_idp
  real (idp), parameter :: tres = 3.0_idp
  real (idp), parameter :: cuat = 4.0_idp
  real (idp), parameter :: cinc = 5.0_idp
  real (idp), parameter :: seis = 6.0_idp
  real (idp), parameter :: siet = 7.0_idp
  real (idp), parameter :: ocho = 8.0_idp
  real (idp), parameter :: nuev = 9.0_idp
  real (idp), dimension (9,9), parameter :: &
    ccc = reshape ( [ &
      uno,  dos, tres, cero, cero, cero, cero, cero, cero, & ! 1
      dos,  uno,  dos, tres, cero, cero, cero, cero, cero, & ! 2
     tres,  dos,  uno,  dos, tres, cero, cero, cero, cero, & ! 3
     cuat, tres,  dos,  uno,  dos, tres, cero, cero, cero, & ! 4
     cero, cuat, tres,  dos,  uno,  dos, tres, cero, cero, & ! 5
     cero, cero, cuat, tres,  dos,  uno,  dos, tres, cero, & ! 6
     cero, cero, cero, cuat, tres,  dos,  uno,  dos, tres, & ! 7
     cero, cero, cero, cero, cuat, tres,  dos,  uno,  dos, & ! 8
     cero, cero, cero, cero, cero, cuat, tres,  dos,  uno  & ! 9
     ], [9,9])
  !
  real (idp), dimension (9,3), parameter :: &
    vvv = reshape ( [ &
      uno,  uno,  uno,  uno,  uno,  uno,  uno,  uno,  uno, & ! 1
      uno, -uno,  uno, -uno,  uno, -uno,  uno, -uno,  uno, & ! 2
      uno,  dos, tres, cuat, cinc, seis, siet, ocho, nuev  & ! 3
      ], [9,3])
  !
  real    (idp), dimension (8,9) :: aaa
  real    (idp), dimension (9,3) :: bbb
  integer (iin), dimension (9)   :: ipiv
  integer (iin) :: mm, nn, klo, kup, nrhs, lda, ldb, info
  integer (iin) :: i
  !
  nn   = size (aaa,2)
  mm   = nn
  klo  = 2
  kup  = 3
  nrhs = size (bbb,2)
  lda  = 2 * klo + kup + 1
  ldb  = nn
  !
  call conversor_agb (mm,nn,lda,kup,klo,ccc,aaa)
  bbb  = vvv
  ipiv = 0
  !
  write (*,*)
  write (*,*)" dgbsv (General Band Side) test "
  write (*,*)
  write (*,*)" system matrix: aaa (lda,n) "
  do  i = 1, lda
  write (*,*) real (aaa (i,:))
  end do
  write (*,*)
  write (*,*)" source: bbb (ldb,nrhs) "
  do  i = 1, ldb
  write (*,*) real (bbb (i,:))
  end do
  write (*,*)" nn   = ", nn
  write (*,*)" klo  = ", klo
  write (*,*)" kup  = ", kup
  write (*,*)" nrhs = ", nrhs
  write (*,*)" lda  = ", lda
  write (*,*)" ldb  = ", ldb
  !           N  KL   KU   NRHS   A   LDA   IPIV  B    LDB  INFO)
  !           |   |    |    |     |    |    |     |      |     |
 !call dgbsv (9 , 2  , 3  , 3   , aaa, 8  , ipiv, bbb,   9, info)
  call dgbsv (nn, klo, kup, nrhs, aaa, lda, ipiv, bbb, ldb, info)
  !
  write (*,*)
  write (*,*)" output aaa: "
  do  i = 1, lda
  write (*,*) real (aaa (i,:))
  end do
  write (*,*)
  write (*,*)" output bbb: "
  do  i = 1, ldb
  write (*,*) real (bbb (i,:))
  end do 
  write (*,*)
  write (*,*)" output ipiv: ", ipiv
  write (*,*)
  write (*,*)" info = ", info
end program extrouble
!end fixed test