Segmentation Fault in Gauss-Legendre Subroutine
Fatemeh
I'm trying to use gauleg (Gauss-Legendre) subroutine of numerical
recipes(Fortran90).
this subroutine makes a bug in the following lines:
x(1:m)=xm-xl*z !Scale the root to the
desired interval,
x(n:n-m+1:-1)=xm+xl*z !and put in its
symmetric counterpart.
w(1:m)=2.0_dp*xl/((1.0_dp-z**2)*pp**2) !Compute the weight
w(n:n-m+1:-1)=w(1:m) !and its symmetric counterpart.
I call this subroutine such as
lower= -PI ; upper=DATAN(2.D0*PI*KB*Temprature*fermipoles/radius)
call gauleg(lower,upper,XX,WW)
and
real(kind=8),allocatable::XX(:),WW(:)
integer,parameter:: Nparameter=20
allocate(XX(Nparameter),WW(Nparameter))
so the dimension of XX , WW is 20.
but in subroutine the program prints for me "n"=1 where n is the size
(x) , size(w)
finally I have segmentation fault at the mentioned lines.
Is there anyone can guide me what is the element which I didn't notice
correctly?
====================
SUBROUTINE gauleg(x1,x2,x,w)
USE nrtype; USE nrutil, ONLY :arth,assert_eq,nrerror
IMPLICIT NONE
REAL(kind=8), INTENT(IN) :: x1,x2
REAL(kind=8), DIMENSION(:), INTENT(OUT) :: x,w
REAL(DP), PARAMETER :: EPS=3.0e-14_dp
!Given the lower and upper limits of integration x1 and x2, this
routine returns arrays x and w
!of length N containing the abscissas and weights of the Gauss-
Legendre N-point quadrature
!formula. The parameter EPS is the relative precision. Note that
internal computations are
!done in double precision.
INTEGER(I4B) :: its,j,m,n
INTEGER(I4B), PARAMETER :: MAXIT=10
REAL(DP) :: xl,xm
REAL(DP), DIMENSION((size(x)+1)/2) :: p1,p2,p3,pp,z,z1
LOGICAL(LGT), DIMENSION((size(x)+1)/2) :: unfinished
n=assert_eq(size(x),size(w),'gauleg')
m=(n+1)/2 !The roots are symmetric in the interval,
!so we only
xm=0.5_dp*(x2+x1) !have to find half of them.
xl=0.5_dp*(x2-x1)
z=cos(PI_D*(arth(1,1,m)-0.25_dp)/(n+0.5_dp)) !Initial
approximations to the roots.
unfinished=.true.
do its=1,MAXIT !Newtonâs method carried out simultane
where (unfinished) !ously on the roots.
p1=1.0
p2=0.0
end where
do j=1,n !Loop up the recurrence relation to get
!the Legendre polynomials evaluated
!at z.
where (unfinished)
p3=p2
p2=p1
p1=((2.0_dp*j-1.0_dp)*z*p2-(j-1.0_dp)*p3)/j
end where
end do
!p1 now contains the desired Legendre polynomials. We next
compute pp, the derivatives,
!by a standard relation involving also p2, the polynomials of
one lower order.
where (unfinished)
pp=n*(z*p1-p2)/(z*z-1.0_dp)
z1=z
z=z1-p1/pp !Newtonâs method.
unfinished=(abs(z-z1) > EPS)
end where
if (.not. any(unfinished)) exit
end do
print *, n, m, its, maxit
if (its == MAXIT+1) call nrerror('too many iterations in
gauleg')
x(1:m)=xm-xl*z !Scale the root to the desired interval,
x(n:n-m+1:-1)=xm+xl*z !and put in its symmetric counterpart.
w(1:m)=2.0_dp*xl/((1.0_dp-z**2)*pp**2) !Compute the weight
w(n:n-m+1:-1)=w(1:m) !and its symmetric counterpart.
END SUBROUTINE gauleg
=====================
Thanks in advance,
Fatemeh
Dick Hendrickson
In order to use a ":" as a dimension, GAULEG needs to
be in a module. Is it? Anytime you want a subroutine
argument to have some hidden property (like it's dimensions)
passed on the call, you'll need to put the subroutine in a
module and USE the module in the main program. Also,
you'll need to compile the module before you compile any
routines that use the module.
On problems where something bad happens with a dummy argument
in a subroutine you should show all of the declarations in
the program that calls the subroutine, otherwise it's hard to
guess what went wrong.
You should get in the habit of putting ALL of your subroutines
and functions in one or more modules. That will prevent most
problems involving dummy arguments.
Dick Hendrickson
PS. Technically, you don't need to put the routine in a
module, you can write something called an interface block,
but that's a really really bad way to do it in this
simple case.
> REAL(DP), PARAMETER :: EPS=3.0e-14_dp
> !Given the lower and upper limits of integration x1 and x2, this
> routine returns arrays x and w
> !of length N containing the abscissas and weights of the Gauss-
> Legendre N-point quadrature
> !formula. The parameter EPS is the relative precision. Note that
> internal computations are
> !done in double precision.
> INTEGER(I4B) :: its,j,m,n
> INTEGER(I4B), PARAMETER :: MAXIT=10
> REAL(DP) :: xl,xm
> REAL(DP), DIMENSION((size(x)+1)/2) :: p1,p2,p3,pp,z,z1
> LOGICAL(LGT), DIMENSION((size(x)+1)/2) :: unfinished
> n=assert_eq(size(x),size(w),'gauleg')
> m=(n+1)/2 !The roots are symmetric in the interval,
> !so we only
> xm=0.5_dp*(x2+x1) !have to find half of them.
> xl=0.5_dp*(x2-x1)
> z=cos(PI_D*(arth(1,1,m)-0.25_dp)/(n+0.5_dp)) !Initial
> approximations to the roots.
> unfinished=.true.
> do its=1,MAXIT !Newtonās method carried out simultane
> where (unfinished) !ously on the roots.
> p1=1.0
> p2=0.0
> end where
> do j=1,n !Loop up the recurrence relation to get
> !the Legendre polynomials evaluated
> !at z.
> where (unfinished)
> p3=p2
> p2=p1
> p1=((2.0_dp*j-1.0_dp)*z*p2-(j-1.0_dp)*p3)/j
> end where
> end do
> !p1 now contains the desired Legendre polynomials. We next
> compute pp, the derivatives,
> !by a standard relation involving also p2, the polynomials of
> one lower order.
> where (unfinished)
> pp=n*(z*p1-p2)/(z*z-1.0_dp)
> z1=z
> z=z1-p1/pp !Newtonās method.
Fatemeh
> > where (unfinished) !ously on the roots.
> > p1=1.0
> > p2=0.0
> > end where
> > do j=1,n !Loop up the recurrence relation to get
> > !the Legendre polynomials evaluated
> > !at z.
> > where (unfinished)
> > p3=p2
> > p2=p1
> > p1=((2.0_dp*j-1.0_dp)*z*p2-(j-1.0_dp)*p3)/j
> > end where
> > end do
> > !p1 now contains the desired Legendre polynomials. We next
> > compute pp, the derivatives,
> > !by a standard relation involving also p2, the polynomials of
> > one lower order.
> > where (unfinished)
> > pp=n*(z*p1-p2)/(z*z-1.0_dp)
> > z1=z
> > unfinished=(abs(z-z1) > EPS)
> > end where
> > if (.not. any(unfinished)) exit
> > end do
> > print *, n, m, its, maxit
> > if (its == MAXIT+1) call nrerror('too many iterations in
> > gauleg')
>
> > x(1:m)=xm-xl*z !Scale the root to the desired interval,
> > x(n:n-m+1:-1)=xm+xl*z !and put in its symmetric counterpart.
> > w(1:m)=2.0_dp*xl/((1.0_dp-z**2)*pp**2) !Compute the weight
> > w(n:n-m+1:-1)=w(1:m) !and its symmetric counterpart.
> > END SUBROUTINE gauleg
> > =====================
>
> > Thanks in advance,
> > Fatemeh
>
>
I have put the subroutine after my program:
program
.
.
.
end program
subroutine gauleg(......)
.
.
end subroutine gauleg
so I don't have any module.
the surprising point for me is that this subroutine doesn't have any
input as a dimension of X , W.
Also according to n=assert_eq(size(x),size(w),'gauleg'), the
subroutine takes n equal to 1.
So I don't know really does this subroutine works correctly or no?
I tried to add an input to subroutine as N
and wrote :
integer,intent(in)::n
REAL(kind=8), dimension(n), INTENT(OUT) :: x(:),w(:)
and also in nr.f90 I made some corrections:
INTERFACE
SUBROUTINE gauleg(x1,x2,x,w,n)
USE nrtype
REAL(kind=8), INTENT(IN) :: x1,x2
integer,intent(in)::n
REAL(kind=8), DIMENSION(n), INTENT(OUT) :: x,w
END SUBROUTINE gauleg
END INTERFACE
and deleted the line n=assert_eq(size(x),size(w),'gauleg')
but it doesn't work and it gives me :
*** glibc detected *** ./a.out: munmap_chunk(): invalid pointer:
0x0830c400 ***
*** glibc detected *** ./a.out: malloc(): memory corruption:
0x0830c430 ***
======= Backtrace: =========
/lib/libc.so.6[0xb7d6d4b6]
/lib/libc.so.6[0xb7d6f701]
/lib/libc.so.6(__libc_malloc+0x97)[0xb7d70d07]
/lib/ld-linux.so.2[0xb7f2e2a7]
/lib/ld-linux.so.2[0xb7f337a4]
/lib/ld-linux.so.2[0xb7f2f786]
/lib/ld-linux.so.2[0xb7f3314c]
/lib/libc.so.6[0xb7dfd4b2]
/lib/ld-linux.so.2[0xb7f2f786]
/lib/libc.so.6[0xb7dfd5b1]
/lib/libc.so.6(__libc_dlopen_mode+0x3b)[0xb7dfd6db]
/lib/libc.so.6[0xb7ddb069]
/lib/libc.so.6(backtrace+0xda)[0xb7ddb1ea]
/lib/libc.so.6[0xb7d67822]
/lib/libc.so.6[0xb7d6d4b6]
/usr/lib/libgfortran.so.2(_gfortran_internal_free+0x21)[0xb7e7a2a1]
./a.out[0x804e29b]
./a.out[0x804ad0b]
./a.out[0x8178a37]
/lib/libc.so.6(__libc_start_main+0xe0)[0xb7d1cfe0]
./a.out[0x8048fb1]
======= Memory map: ========
08048000-0817a000 r-xp 00000000 08:03 3211376 /data/fmirjani/
fortran/LDA/exact-HM/l2/100/a.out
0817a000-0817b000 r-xp 00132000 08:03 3211376 /data/fmirjani/
fortran/LDA/exact-HM/l2/100/a.out
0817b000-0817c000 rwxp 00133000 08:03 3211376 /data/fmirjani/
fortran/LDA/exact-HM/l2/100/a.out
0817c000-0832b000 rwxp 0817c000 00:00 0 [heap]
b6200000-b6221000 rwxp b6200000 00:00 0
b6221000-b6300000 ---p b6221000 00:00 0
b63d8000-b7d07000 rwxp b63d8000 00:00 0
b7d07000-b7e34000 r-xp 00000000 08:02 533 /lib/libc-2.6.1.so
b7e34000-b7e35000 r-xp 0012c000 08:02 533 /lib/libc-2.6.1.so
b7e35000-b7e37000 rwxp 0012d000 08:02 533 /lib/libc-2.6.1.so
b7e37000-b7e3a000 rwxp b7e37000 00:00 0
b7e3a000-b7e44000 r-xp 00000000 08:02 609 /lib/libgcc_s.so.1
b7e44000-b7e46000 rwxp 00009000 08:02 609 /lib/libgcc_s.so.1
b7e46000-b7e69000 r-xp 00000000 08:02 541 /lib/libm-2.6.1.so
b7e69000-b7e6b000 rwxp 00022000 08:02 541 /lib/libm-2.6.1.so
b7e6b000-b7f00000 r-xp 00000000 08:02 114899 /usr/lib/
libgfortran.so.2.0.0
b7f00000-b7f02000 rwxp 00095000 08:02 114899 /usr/lib/
libgfortran.so.2.0.0
b7f02000-b7f03000 rwxp b7f02000 00:00 0
b7f22000-b7f3c000 r-xp 00000000 08:02 524 /lib/ld-2.6.1.so
b7f3c000-b7f3e000 rwxp 0001a000 08:02 524 /lib/ld-2.6.1.so
bfac7000-bfadc000 rwxp bfac7000 00:00 0 [stack]
ffffe000-fffff000 r-xp 00000000 00:00 0 [vdso]
Aborted
Is there anyone can help me? I've got really confused with this
subroutine!
Best,
Fatemeh
Richard Maine
> On Apr 11, 5:26 pm, Dick Hendrickson <dick.hendrick...@att.net> wrote:
> > In order to use a ":" as a dimension, GAULEG needs to
> > be in a module.
..
> > PS. Technically, you don't need to put the routine in a
> > module, you can write something called an interface block,
> > but that's a really really bad way to do it in this
> > simple case.
> so I don't have any module.
Well, then that's your problem. As Dick said, you need one.
> the surprising point for me is that this subroutine doesn't have any
> input as a dimension of X , W.
> Also according to n=assert_eq(size(x),size(w),'gauleg'), the
> subroutine takes n equal to 1.
> So I don't know really does this subroutine works correctly or no?
Dick is a pretty smart guy. Take his advice instead of wandering off and
trying other more complicated things. These questions are all related.
The (:) used for the dimension in the dummy arguments means assumed
shape. With assumed shape, the compiletr will take care of passing the
dimensions for you. But that requires an explicit interface (which is
most easily done by putting the routine an a module).
> I tried to add an input to subroutine as N
That's what I mean by wandering off and trying other more complicated
things.
> and also in nr.f90 I made some corrections:
> INTERFACE
> SUBROUTINE gauleg(x1,x2,x,w,n)
> USE nrtype
> REAL(kind=8), INTENT(IN) :: x1,x2
> integer,intent(in)::n
> REAL(kind=8), DIMENSION(n), INTENT(OUT) :: x,w
> END SUBROUTINE gauleg
> END INTERFACE
You didn't mention this before. It is *extremely* relevant. As Dick also
mentioned, having an interface body (such as this) is an alternative to
putting the subroutine in a module. But that only works if the interface
body is in the calling routine or accessible there, possibly by being in
a module that is USEd there.
You are omitting all of the critical context. What is this interface
body in in nr.f90? It can't just stand alone; that won't compile (unless
perhaps it is in a file meant to be INCLUDEd in your source instead of
compiled separately).
I'd suggest reading whatever documentation comes with the nr routines.
Presumably you also have the book (as I thought I recalled it was
illegal to distribut ethe code separately from the book). I didn't drag
it out to check, but it ought to have directions on how to use the
routines. Those directions, in turn ought to say something about USEing
the appropriate module, or perhaps INCLUDEing a file. I would be really
quite shocked if such directions weren't there.
I suggest following those directions, which probably match Dick's
comment.
> Is there anyone can help me? I've got really confused with this
> subroutine!
Stop trying to modify the subroutine. Instead concentrate on the
directions about how to correctly use it, likley by USEing an
appropriate module.
--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
robin
news:a4babd34-c21a-457d...@n1g2000vba.googlegroups.com...
>I have put the subroutine after my program:
>program
>.
>.
>end program
>subroutine gauleg(......)
>.
>.
>end subroutine gauleg
>so I don't have any module.
>the surprising point for me is that this subroutine doesn't have any
>input as a dimension of X , W.
>Also according to n=assert_eq(size(x),size(w),'gauleg'), the
>subroutine takes n equal to 1.
>So I don't know really does this subroutine works correctly or no?
You won't know until you put the subroutine in a module.
>I tried to add an input to subroutine as N
>and wrote :
You need a module.
Please don't post again until you have done that.
>Is there anyone can help me? I've got really confused with this
>subroutine!
Dick Hendrickson told you what to do. Now, DO it.
Or, as they say in the movies, "Do it NOW". :)