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
