COMPGROUPS.NET | Search | Post Question | Groups | Stream | About | Register

### The dgetri subroutine from Lapack return segment failure

• Follow

```!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.
```
 0
Reply leed2005 (11) 9/22/2009 12:31:48 PM

```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.

```
 0
Reply g.sande (1183) 9/22/2009 1:31:34 PM

```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

```
 0
Reply epc8 (1259) 9/22/2009 2:19:49 PM

```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
```
 0
Reply pastgio1 (24) 9/22/2009 7:42:51 PM

```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)
---------------------------------------------------------------------------=
--------------

```
 0
Reply leed2005 (11) 9/23/2009 12:40:45 AM

```Finally, I found that why dgetri subroutine can't work.

The array a refer to dgetri should after LU decomposition by dgetrf.

```
 0
Reply leed2005 (11) 9/24/2009 12:31:35 AM

5 Replies
45 Views

Similiar Articles:

7/25/2012 2:42:50 AM