The dgetri subroutine from Lapack return segment failure

!My test program is listed as follow: !----------------------------------------------------------------- program main implicit none integer, parameter :: n =3D 3 integer, parameter :: lda =3D n real, dimension( 1 : n, 1 : n ) :: a real, dimension( 1 : n, 1 : n ) :: b integer, dimension( 1 : n ) :: ipiv integer, parameter :: lwork =3D 1000 real(kind=3D8), dimension( lwork ) :: work integer :: info integer :: i, j a( 1, 1 ) =3D 5 a( 1, 2 ) =3D 0 a( 1, 3 ) =3D 0 a( 2, 1 ) =3D 0 a( 2, 2 ) =3D 3 a( 2, 3 ) =3D 1 a( 3, 1 ) =3D 0 a( 3, 2 ) =3D 2 a( 3, 3 ) =3D 1 b =3D a print "('old',/,3(f6.2,1x))", ((a(i, j), j =3D 1, n), i =3D 1, n) call dgetri(n, a, lda, ipiv, work, lwork, info) print "('inv',/,3(f6.2,1x))", ((b(i, j), j =3D 1, n), i =3D 1, n) b =3D matmul(a, b) print "('mul',/,3(f6.2,1x))", ((b(i, j), j =3D 1, n), i =3D 1, n) end program !The result is: !old ! 5.00 0.00 0.00 ! 0.00 3.00 1.00 ! 0.00 2.00 1.00 !inv ! 5.00 0.00 NaN ! 0.00 NaN -0.00 ! 0.00 0.00 2.14 !mul ! NaN NaN NaN ! NaN NaN NaN ! 0.00 NaN NaN !=E6=AE=B5=E9=94=99=E8=AF=AF I don't know how to find out the reason caused this kind error. I apprecaite any suggestion.

On 2009-09-22 09:31:48 -0300, leed <leed2005@gmail.com> said: > !My test program is listed as follow: > !----------------------------------------------------------------- > program main > implicit none > integer, parameter :: n = 3 > integer, parameter :: lda = n > real, dimension( 1 : n, 1 : n ) :: a > real, dimension( 1 : n, 1 : n ) :: b a and b are single > integer, dimension( 1 : n ) :: ipiv > integer, parameter :: lwork = 1000 > real(kind=8), dimension( lwork ) :: work work is double assuming kind=8 is double > integer :: info > integer :: i, j > > a( 1, 1 ) = 5 > a( 1, 2 ) = 0 > a( 1, 3 ) = 0 > a( 2, 1 ) = 0 > a( 2, 2 ) = 3 > a( 2, 3 ) = 1 > a( 3, 1 ) = 0 > a( 3, 2 ) = 2 > a( 3, 3 ) = 1 > > b = a > print "('old',/,3(f6.2,1x))", ((a(i, j), j = 1, n), i = 1, n) > call dgetri(n, a, lda, ipiv, work, lwork, info) double version with single argument. type mismatch! > print "('inv',/,3(f6.2,1x))", ((b(i, j), j = 1, n), i = 1, n) > b = matmul(a, b) > print "('mul',/,3(f6.2,1x))", ((b(i, j), j = 1, n), i = 1, n) > end program > > !The result is: > !old > ! 5.00 0.00 0.00 > ! 0.00 3.00 1.00 > ! 0.00 2.00 1.00 > !inv > ! 5.00 0.00 NaN > ! 0.00 NaN -0.00 > ! 0.00 0.00 2.14 > !mul > ! NaN NaN NaN > ! NaN NaN NaN > ! 0.00 NaN NaN > !段错误 > > I don't know how to find out the reason caused this kind error. > I apprecaite any suggestion.

On Sep 22, 8:31=A0am, leed <leed2...@gmail.com> wrote: > !My test program is listed as follow: > !----------------------------------------------------------------- > program main > implicit none > integer, parameter :: n =3D 3 > integer, parameter :: lda =3D n > real, dimension( 1 : n, 1 : n ) :: a > real, dimension( 1 : n, 1 : n ) :: b > integer, dimension( 1 : n ) :: ipiv > integer, parameter :: lwork =3D 1000 > real(kind=3D8), dimension( lwork ) :: work > integer :: info > integer :: i, j > > a( 1, 1 ) =3D 5 > a( 1, 2 ) =3D 0 > a( 1, 3 ) =3D 0 > a( 2, 1 ) =3D 0 > a( 2, 2 ) =3D 3 > a( 2, 3 ) =3D 1 > a( 3, 1 ) =3D 0 > a( 3, 2 ) =3D 2 > a( 3, 3 ) =3D 1 > > b =3D a > print "('old',/,3(f6.2,1x))", ((a(i, j), j =3D 1, n), i =3D 1, n) > call dgetri(n, a, lda, ipiv, work, lwork, info) at the top of the source code from NETLIB, this is listed: SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) [snip] * .. Scalar Arguments .. INTEGER INFO, LDA, LWORK, N * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ), WORK( * ) * .. The routine you called expects A to be of type DOUBLE PRECISION but your array is SINGLE PRECISION (REAL). An argument mismtach is a frequent cause of troubles in programs posted here! Also note that real(kind=3D8) is not portable. It MAY be DOUBLE PRECISION on your compiler, on mine the kind number might be 2 or something else. The kind number is not necessarily the number of bytes occupied by the data item. IMO it's better to use DOUBLE PRECISION or one of the more modern methods of specifying kind. > print "('inv',/,3(f6.2,1x))", ((b(i, j), j =3D 1, n), i =3D 1, n) > b =3D matmul(a, b) > print "('mul',/,3(f6.2,1x))", ((b(i, j), j =3D 1, n), i =3D 1, n) > end program > > !The result is: > !old > ! =A05.00 =A0 0.00 =A0 0.00 > ! =A00.00 =A0 3.00 =A0 1.00 > ! =A00.00 =A0 2.00 =A0 1.00 > !inv > ! =A05.00 =A0 0.00 =A0 =A0NaN > ! =A00.00 =A0 =A0NaN =A0-0.00 > ! =A00.00 =A0 0.00 =A0 2.14 The garbage returned here is a good clue that you have written data somewhere it was not meant to go. So your called routine will read garbage because A occupies only half as much space as it should AND the called routine will write beyond A because it is only half as large as it should be. In a larger program you might see various other types of unexpected behavior. > !mul > ! =A0 NaN =A0 =A0NaN =A0 =A0NaN > ! =A0 NaN =A0 =A0NaN =A0 =A0NaN > ! =A00.00 =A0 =A0NaN =A0 =A0NaN GIGO. Garbage in - garbage out. -- e

e p chandler wrote: > On Sep 22, 8:31 am, leed <leed2...@gmail.com> wrote: >> !My test program is listed as follow: .... >> real, dimension( 1 : n, 1 : n ) :: a .... >> call dgetri(n, a, lda, ipiv, work, lwork, info) .... > The routine you called expects A to be of type DOUBLE PRECISION but > your array is SINGLE PRECISION (REAL). > An argument mismtach is a frequent cause of troubles in programs > posted here! Maybe one could add that lapack implements quite a regular naming convention for its subprograms: initial letter --- subprogram works on s real d double precision c complex z double precision complex In addition, lapack95 implements generic subprograms. Giorgio

Thank you all for your nice reply. I test sgetri and dgetri subroutines respectively according above suggestion. First test sgetri ---------------------------------------------------------------------------= -------------- program main implicit none integer, parameter :: n =3D 3 integer, parameter :: lda =3D n real, dimension( 1 : n, 1 : n ) :: a real, dimension( 1 : n, 1 : n ) :: b integer, dimension( 1 : n ) :: ipiv integer, parameter :: lwork =3D 1000 real, dimension( lwork ) :: work integer :: info integer :: i, j a( 1, 1 ) =3D 5 a( 1, 2 ) =3D 0 a( 1, 3 ) =3D 0 a( 2, 1 ) =3D 0 a( 2, 2 ) =3D 3 a( 2, 3 ) =3D 1 a( 3, 1 ) =3D 0 a( 3, 2 ) =3D 2 a( 3, 3 ) =3D 1 b =3D a print "('old',/,3(f6.2,1x))", ((a(i, j), j =3D 1, n), i =3D 1, n) call sgetri(n, a, lda, ipiv, work, lwork, info) print "('info=3D',I3)", info print "('inv',/,3(f6.2,1x))", ((a(i, j), j =3D 1, n), i =3D 1, n) b =3D matmul(a, b) print "('mul',/,3(f6.2,1x))", ((b(i, j), j =3D 1, n), i =3D 1, n) end program Result is incrrect, and is listed as follow: old 5.00 0.00 0.00 0.00 3.00 1.00 0.00 2.00 1.00 info=3D 0 inv 0.00 0.00 -0.00 1.00 1.00 -0.33 -2.00 1.00 1.00 mul 0.00 0.00 0.00 5.00 2.33 0.20 -10.00 5.00 -0.40 ---------------------------------------------------------------------------= -------------- Second, the dgetri ---------------------------------------------------------------------------= -------------- program main implicit none integer, parameter :: n =3D 3 integer, parameter :: lda =3D n double precision, dimension( 1 : n, 1 : n ) :: a double precision, dimension( 1 : n, 1 : n ) :: b integer, dimension( 1 : n ) :: ipiv integer, parameter :: lwork =3D 1000 double precision, dimension( lwork ) :: work integer :: info integer :: i, j a( 1, 1 ) =3D 5 a( 1, 2 ) =3D 0 a( 1, 3 ) =3D 0 a( 2, 1 ) =3D 0 a( 2, 2 ) =3D 3 a( 2, 3 ) =3D 1 a( 3, 1 ) =3D 0 a( 3, 2 ) =3D 2 a( 3, 3 ) =3D 1 b =3D a print "('old',/,3(f6.2,1x))", ((a(i, j), j =3D 1, n), i =3D 1, n) call dgetri(n, a, lda, ipiv, work, lwork, info) print "('info=3D',I3)", info print "('inv',/,3(f6.2,1x))", ((a(i, j), j =3D 1, n), i =3D 1, n) b =3D matmul(a, b) print "('mul',/,3(f6.2,1x))", ((b(i, j), j =3D 1, n), i =3D 1, n) end program The result is listed as follow: old 5.00 0.00 0.00 0.00 3.00 1.00 0.00 2.00 1.00 =E6=AE=B5=E9=94=99=E8=AF=AF(segment failure) ---------------------------------------------------------------------------= --------------