Dear all ;

After running my program, I see at the end of the execution some lines
such as the following :

*** glibc detected *** ./a.out: free(): invalid next size (fast):
0x08199df0 ***
======= Backtrace: =========
/lib/libc.so.6[0xb7d854b6]
/lib/libc.so.6(cfree+0x89)[0xb7d87179]
/usr/lib/libgfortran.so.2(_gfortran_deallocate+0x29)[0xb7e92239]
./a.out[0x805c034]
./a.out[0x818e407]
/lib/libc.so.6(__libc_start_main+0xe0)[0xb7d34fe0]
./a.out[0x8048e31]
======= Memory map: ========
08048000-08190000 r-xp 00000000 08:03 3196154    /data/fmirjani/
fortran/LDA/test/a.out
08190000-08191000 r-xp 00147000 08:03 3196154    /data/fmirjani/
fortran/LDA/test/a.out
08191000-08192000 rwxp 00148000 08:03 3196154    /data/fmirjani/
fortran/LDA/test/a.out
08192000-081b3000 rwxp 08192000 00:00 0          [heap]
b7c00000-b7c21000 rwxp b7c00000 00:00 0
b7c21000-b7d00000 ---p b7c21000 00:00 0
b7d1e000-b7d1f000 rwxp b7d1e000 00:00 0
b7d1f000-b7e4c000 r-xp 00000000 08:02 533        /lib/libc-2.6.1.so
b7e4c000-b7e4d000 r-xp 0012c000 08:02 533        /lib/libc-2.6.1.so
b7e4d000-b7e4f000 rwxp 0012d000 08:02 533        /lib/libc-2.6.1.so
b7e4f000-b7e52000 rwxp b7e4f000 00:00 0
b7e52000-b7e5c000 r-xp 00000000 08:02 609        /lib/libgcc_s.so.1
b7e5c000-b7e5e000 rwxp 00009000 08:02 609        /lib/libgcc_s.so.1
b7e5e000-b7e81000 r-xp 00000000 08:02 541        /lib/libm-2.6.1.so
b7e81000-b7e83000 rwxp 00022000 08:02 541        /lib/libm-2.6.1.so
b7e83000-b7f18000 r-xp 00000000 08:02 114899     /usr/lib/
libgfortran.so.2.0.0
b7f18000-b7f1a000 rwxp 00095000 08:02 114899     /usr/lib/
libgfortran.so.2.0.0
b7f1a000-b7f1b000 rwxp b7f1a000 00:00 0
b7f3a000-b7f54000 r-xp 00000000 08:02 524        /lib/ld-2.6.1.so
b7f54000-b7f56000 rwxp 0001a000 08:02 524        /lib/ld-2.6.1.so
bffc3000-bffd8000 rwxp bffc3000 00:00 0          [stack]
ffffe000-fffff000 r-xp 00000000 00:00 0          [vdso]
Aborted

Would you please explain me what are these lines telling ? (I use
gfortran for compiling.)

SWAG: you're trying to deallocate an already unallocated array or
you're trying to deallocate an array that was never allocated.

      program start
      !use nrtype
      !use nrutil
      implicit none

      real(kind=8)::betta,func1,func
      REAL(kind=8)::bessj0_s,bessj1_s,x,x1,x2,y,xacc
      REAL(kind=8)::sumn,frac
      REAL(kind=8)::rtsec
      REAL(kind=8)::z,answer
      real(kind=8)::u,energy
      real(kind=8)::integralarguman,t,hbar,a,mass
      real
(kind=8),parameter::PI=3.141592653589793238462643383279502884197
      real(kind=8),allocatable::XX(:),WW(:)

      double precision,allocatable::hamiltonian(:,:),amatrix(:,:),W
(:),WORK(:),nan(:,:),matr(:),Vxc(:)

      integer::lda,lwork,info
      integer::nelec,i,j,iteration,h,nn,NMAX,mcounter
      integer,parameter:: Nparameter=20
      integer:: iter,icnt
!=======================================reading the input from
terminal

      u=0.D0
      allocate(XX(Nparameter),WW(Nparameter))
      print*,"enter the dimension of matrix,nelec:"
      read*,nn,nelec
     !number of iteration 500
      allocate(nan(1:nn,1:500),Vxc(nn))   !500?
      LWORK=max(1,3*nn-1) ; LDA=max(1,nn)
      NMAX=nn  !check it later  !ASK IT,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
      allocate(W(NMAX),WORK(LWORK))

      allocate (hamiltonian(nn,nn),amatrix(LDA,NMAX))

!=========================================
      !first block
      do while(u.lt.20.D0)
         energy=0
        !hbar=6.58211899*E-16 !eVs
        ! hbar=1.054571628E-34 ! Js
        !mass=0.25*9.10938188E-31 !Kg
        !a=3.0E-10 !0.3nm
        !t=DBLE(hbar**2)/DBLE(2*mass*(a**2))
        t=1.D0
        iteration=2
        do i=1,nn  !initial value for n(i)
          nan(i,iteration)=DBLE(nelec)/DBLE(nn)
          nan(i,1)=0
        end do
!calculate the integral===================================
      frac=u/(2*t)
      if (u==0.D0) then
        betta=2.D0
      else
      !iter=30
      ! do icnt=1,iter
       call LAGZO(Nparameter,XX,WW)
       answer=0.D0
       do i=1,Nparameter!icnt
         answer=answer+Func(XX(i),frac)*WW(i)
       end do
      !end do
! finding betta(u/t)====================================
          y=1.9D0
          if ((func1(y,answer)).gt.0.D0)then
          x2=y
           do while((func1(y,answer)).gt.0.D0)
           y=y-0.01
           end do
          x1=y
          end if
          if((u==0.1D0).or.(u==0.2D0))print*,func1(1.9D0,answer),func1
(1.D0,answer)
         !print*,x1,x2
          xacc=3  !????????????????????????????????????????
          !x1=1.1D0;x2=1.999d0
         betta=rtsec(x1,x2,xacc,answer)
       end if! if(u==0.D0)
       !print*,u,betta
!===================================== start the loop , since finding
beta and integral answer is independent of the loop (only dependent to
u)
!this part calculates the Vxc(i)
!check the convergency criterion=======
       sumn=0
       do i=1,nn
        sumn=(nan(i,iteration)-nan(i,(iteration-1)))**2+sumn
       end do
       !if(iteration==2) print*,sumn,"ll"
       do while (sumn.gt.10E-05) !convergency criterion
       Vxc(i)=-2*t*cos(PI*nan(i,iteration)/betta)+(2*t*cos(PI*nan
(i,iteration)/2))-(u*nan(i,iteration)/2)
!=====================================

!constructing the hamiltonian and Diagonalization part

       do i=1,nn
        do j=1,nn
         if (i==j)   then
          hamiltonian(i,j)=2*t+(u/2)*nan(i,iteration)+Vxc(i)
          else if ((i==j-1).or.(i==j+1)) then
           hamiltonian(i,j)=-t
          else
           hamiltonian(i,j)=0
          end if
          !amatrix(i,j)=hamiltonian(i,j)
        end do !j
       end do!i
       do i=1,nn
        do j=i,nn
          amatrix(i,j)=hamiltonian(i,j)
        end do
       end do
      !print*,((amatrix(I,J),J=I,nn),I=1,nn)
      !diagonalization
      CALL DSYEV('V','U',nn,amatrix,LDA,W,WORK,LWORK,INFO)
      !  print*,"test"
      write(12,*)(W(J),J=1,nn) !eigenvalue
      write(13,*)((amatrix(I,J),J=1,nn),I=1,nn) !eigenvectors
!calculate ni from psi (eigenvectors)====== !nan(i)=?  !I suppose that
the nr. of eigenvectors is m
      iteration=iteration+1
      do i=1,nn ! for all of the sites: i ha ,we calculate nan(i)
       nan(i,iteration)=0
       do mcounter=1,nn  !nn= m  !is the number of m (eigen vectors)
equal to nn ? I think yes
         nan(i,iteration)=nan(i,iteration)+amatrix(i,mcounter)**2
       end do !m
      end do   ! i

      sumn=0
       do i=1,nn
        sumn=(nan(i,iteration)-nan(i,(iteration-1)))**2+sumn
       end do

      end do !related to convergency criterion
!===========================
!calculate energy for each u
      do i=1,nn
       energy=W(i)+energy
        ! print*,i,W(i)

      end do
      write(77,*)u,energy
!---------------------------------------------------
       u=u+0.1D0
      end do !u

      end program
      !don't forget to define the above arguments in main & define
xacc
!-------------------------------------------------------=======

       !  CONTAINS
     !SUBROUTINE NewTrapzd(a, b, Result, N, yy)
     !REAL(Kind=8) :: A, B, X, YY, Result, Del
     !INTEGER :: I, N
     !Result = 0.5D0*(func(A,YY) + func(B,YY))
     !Del = (B-A)/DBLE(N-1)
     !DO I=1, N-2
      ! X = A + I*Del
      ! Result = Result + Del*func(X,YY)
      !END DO
      !Result = Result * Del
      !END SUBROUTINE NewTrapzd
     !++++++++++++++++++++++++++++++++++++++++++++
     !       SUBROUTINE trapzd(a,b,s,n,yy)
      !      use nrtype
        !This routine computes the nth stage of refinement of an
extended trapezoidal rule. func is
        !input as the name of the function to be integrated between
limits a and b, also input. When
        !called with n=1, the routine returns as s the crudest
estimate of b
        !a f(x)dx. Subsequent
        !calls with n=2,3,... (in that sequential order) will improve
the accuracy of s by adding 2n-2
        !additional interior points. s should not be modified between
sequential calls
      ! IMPLICIT NONE
     !  REAL(SP), INTENT(IN) :: a,b,yy
     !  REAL(SP), INTENT(INOUT) :: s
     !  INTEGER(I4B), INTENT(IN) :: n
     !  real::func
       !INTERFACE
       !FUNCTION func(x)
       !USE nrtype
       !REAL(SP), DIMENSION(:), INTENT(IN) :: x
       !REAL(SP), DIMENSION(size(x)) :: func
       !END FUNCTION func
       !END INTERFACE
     !  REAL(SP) :: del,fsum
     !  INTEGER(I4B) :: it
     !  if (n == 1) then
     !  s=0.5_sp*(b-a)*(func(a,yy)+func(b,yy))
     !  else
     !  it=2**(n-2)
     !  del=(b-a)/it !This is the spacing of the points to be added.
     !  fsum=sum(func(arth(a+0.5_sp*del,del,it),yy))
     !  s=0.5_sp*(s+del*fsum) !This replaces s by its refined value.
     !  end if
     !  END SUBROUTINE trapzd
!-------------------------------------------------------
      SUBROUTINE LAGZO(N,X,W)

        do i=1,nn
         sumn=(nan(i,iteration)-nan(i,(iteration-1)))**2+sumn
        end do
        !if(iteration==2) print*,sumn,"ll"
        do while (sumn.gt.10E-05) !convergency criterion
          Vxc(i)=...nan(i,iteration)...

which is undertandable, because at the moment of evaluation, I has
value of NN+1 (leftover from the previous loop), but nan is allocated
to (NN,500).

When I move the problematic Vxc(i) line into the next loop (was this
the intent?), the program gets stuck in an endless loop. I don't have
a slightest idea what it is supposed to do.

Please turn on bounds checking, whichever the syntax is.

real (kind=8),parameter::PI=3.1415926

To have it be double precision you need to append "_8" to the
end of the constant.  As I think others have said, the use of
8 for a kind type isn't portable.  You'd be better off to do
something like
     integer, parameter ::  dp = selected_real_kind(what you want)
Or at least
     integer, parameter ::  dp = 8
(not portable, but you'll only have to change one thing when you
have a portability problem)
and then
real (kind=dp),parameter::PI=3.1415926.....197_dp

If you compile your program with "-g" (debugging symbols) and possibly
with "-O0" (no optimization), you will get a better output of the crash
location than:

Running your program through "valgrind" might give a better error
location as it might detect the error itself and not the later consequences.

The error looks as you somehow managed to corrupt a pointer, e.g. by
assigning values to an array outside its bounds or by assigning to an
unallocated/unassociated ALLOCATABLE/POINTER variable.

(Even if you don't have an explicit allocation, variables such as
automatic arrays need to get allocated at some point; the compiler
handles this automatically/internally by calling the relevant GLIBC
routines.)

As other have written, it makes sense to turn on all debugging options,
especially -fbounds-check.

Good luck in debugging the program.

Tobias

You seem to be confused about WHAT is being calculated in WHICH loop.
For example, the coefficients for Gauss-Laguerre integration are
calculated for EACH value of U, but they depend only on N (actually
Nparameter). So they are a "loop invariant".

Also, I repeat my question about your integrand used with Gauss
Laguerre integration. You have changed it from that in a different
thread, but are you really certain that Gauss-Laguerre is really the
proper technique for an equation of the form

  func=(besj0(x/yy)*besj1(x/yy))/(x*(1+exp(-x))) ?

I am ASSUMING that you have factored (exp(-x)) out of func BEFORE
subjecting it to G-L integration.

Next, you have a 4 way branch deal with singularities in your
integrand. I added PRINT *... and STOP statements to the three
singularity cases. None of these were invoked while running your
program.

 2. I am suspicious of the method used with func1.

         y=1.9D0
          if ((func1(y,answer)).gt.0.D0)then
          x2=y
           do while((func1(y,answer)).gt.0.D0)
           y=y-0.01
           end do
          x1=y
          end if
          if((u==0.1D0).or.(u==0.2D0))print*,func1
(1.9D0,answer),func1
(1.D0,answer)
         !print*,x1,x2
          xacc=3  !????????????????????????????????????????
          !x1=1.1D0;x2=1.999d0
         betta=rtsec(x1,x2,xacc,answer)
       end if! if(u==0.D0)
       !print*,u,betta

RTSEC is a subroutine for the SECANT method.

func1 is

     real(kind=8) function func1(x,ans)
       !use nrtype
       implicit none
       real(kind=8),intent(in)::x
       real(kind=8),intent(in)::ans
       real(kind=8)::PI
       PI=3.141592653589793238462643383279502884197
       func1=(x*sin(PI/x))-(2*PI*ans)
       end function func1

Again note that exending the length of the constant for PI does not
make it more precise WITHOUT affixing a kind indicator such as D0 or
_DP, where DP is appropriately defined.

More important, you are trying to solve an equation

x*sin(PI/x) = 2*PI*A

given A, find x.

If you plot x*sin(PI/x) you will find that it oscillates for small x.
Also multiple solutions are possible for many values of A, and none
for some values of A.

As Acton discusses in his text that I recently referenced in a
different but recent thread, these types of equations are problematic.
Simply printing out "solutions" from test runs of your programs shows
that some of them are returned as NaN.

If you look at the series expansion for

x/(2*pi) * sin(PI/x), you will see that it approaches 1/2 for
relatively small (>1) values of x in a near parabolic fashion.

If you substitute v = PI/x, you get

sin(v) = 2*A * v

if you plot

sin(v) and B*v on the same graph, you will see there are problems
including solutions for v at 0, and multiple or (no, actually 1 at 0)
solution depending on B.

[I understand that your main concern is getting your program to run,
but IMO you also should be concerend with correct output. Again, IMO
breaking this type of problem into smaller components and verifying
these components individually by attaching them to appropriate test
programs is a good strategy.]

Without more investigation of the numerical aspcets of your problem, I
suspect you may in fact be wasting time and effort.
Some exprets say that in the end, you produce more verification and
testing code than actual product!