Fortran 90, problem with reading a file

[fluidistic]

Here's a part of my program in which I'm having an error while compiling:

open(unit=1,file="tabla.dat")
allocate (x(n),y(n),f(n))
do i=0,n
x=(2*i)/n-1
f=1/(1+25*x**2)
write(1,*)x,f
end do

do i=0,n
read(1,*)x(i), y(i)
end do
close(1)

I'm having the error with gfortran: "Fortran runtime error: End of file" for the line "read(1,*)x(i), y(i)".

Basically if I enter n=5, my program opens a file, write in it 6 values for x and 6 values for f. I want my program to read them and to assignate the first value of x to x(0), the first value of f to y(0); the second value of x: x(1), second value of f: y(1), etc.
I've opened my file "tabla.dat" and it looks like this:

-1.0000000      -1.0000000      -1.0000000      3.84615399E-02  3.84615399E-02  3.84615399E-02
  -1.0000000      -1.0000000      -1.0000000      3.84615399E-02  3.84615399E-02  3.84615399E-02
   0.0000000       0.0000000       0.0000000       1.0000000       1.0000000       1.0000000    
   1.0000000       1.0000000       1.0000000      3.84615399E-02  3.84615399E-02  3.84615399E-02

which seems ok since to me (well maybe not numerically but at least its form looks ok to me).
Any idea what's wrong/going on?
Thanks in advance.

 

[Useful nucleus]

(1)You need to use REWIND to let the code go back to the first line in the file and read the values of x and y.

(2)Also, I think the size of your arrays in the allocate statement should be (n+1) rather than (n).

(3)Finally, in the first loop be careful that you keep evaluating all the elements of the arrays x and f and then overwrite them again and again.
I believe that the first loop should be:
do i=0,n
x(i)=(2*i)/n-1
f(i)=1/(1+25*x**2)
write(1,*)x(i),f(i)
end do

 

[fluidistic]

(1)You need to use REWIND to let the code go back to the first line in the file and read the values of x and y.

Then I probably messed up my program but if I can fix my problem with it I'd be glad to use it. I'll try it.

(2)Also, I think the size of your arrays in the allocate statement should be (n+1) rather than (n).

You're right. I accidentally erased my n+1 to n due to another problem in my program.

(3)Finally, in the first loop be careful that you keep evaluating all the elements of the arrays x and f and then overwrite them again and again.
I believe that the first loop should be:
do i=0,n
x(i)=(2*i)/n-1
f(i)=1/(1+25*x**2)
write(1,*)x(i),f(i)
end do

Thank you, I understand. In fact it would be "f(i)=1/(1+25*x(i)**2)"
I could compile the program without the rewind command and I got a "*** glibc detected *** ./a.out: double free or corruption (!prev): 0x08a39e18 ***" error.

I'm done for today and I might come back to this program tomorrow (I hope so).
If you're interested, it's a program to interpolate a set of given points with Lagrange's polynomial.

program test
implicit none
real:: interval,h,x_m,x_n,p,l
integer:: n,j, i,m
real, dimension(:), allocatable :: x,y,f


write(*,*)"Enter the dimension n of the polynomial and the number of divisions of the interval"
read(*,*)n,m

open(unit=1,file="tabla.dat")
allocate (x(n+1),y(n+1),f(n+1))
do i=0,n
x(i)=(2*i)/n-1
f(i)=1/(1+25*x(i)**2)
write(1,*)x(i),f(i)
end do

do i=0,n
read(1,*)x(i), y(i)
end do
close(1)
interval=x(n)-x(0)
h=interval/m
x_m=x(0)
x_n=x(n)

do while (x_m<= x_n)
call Lagrange(x_m,p)
open(unit=2,file="result.dat")
write(2,*)x_m,p
x_m=x_m+h
end do
close(2)


Contains
subroutine Lagrange (x_m,p)
implicit none
real, intent(in):: x_m
real, intent(out)::p
integer:: i, j
p=0.
do i=0,n
  if (j/=i) then
  do j=0, n
  l=l*(x_m-x(j))/(x(i)-x(j))
  end do
  end if
p=p+y(i)*l
end do
end subroutine

end program

My inputs would be n=5, 10 and then 15. For a much greater m, say 100 or 1000. I'd also plot the function $f(x)=\frac{1}{1+25x^2}$ along with the date of "result.dat" to see if my polynomial more or less fits f.
I'll probably post here tomorrow if I need help. Thanks once again.

 

[fluidistic]

Finally back to this!
I've included a rewind statement, though I'm not really understanding why I need it.

So my code becomes

open(unit=1,file="tabla.dat")
allocate (x(n+1),y(n+1),f(n+1))
do i=0,n
x(i)=(2*i)/n-1
f(i)=1/(1+25*x(i)**2)
write(1,*)x(i),f(i)
rewind (unit=1)
end do

And I still get a "nasty" error. I've tried to put this statement in the "read" line too and in both lines and still get this error.

do i=0,n
read(1,*)x(i), y(i)
rewind (unit=1)
end do

The error is of the form...

 Enter the dimension n of the polynomial and the number of divisions of the interval
5, 100
*** glibc detected *** ./a.out: double free or corruption (!prev): 0x09c57e18 ***
======= Backtrace: =========
/lib/i386-linux-gnu/libc.so.6(+0x6b961)[0x438961]
/lib/i386-linux-gnu/libc.so.6(+0x6d28b)[0x43a28b]
/lib/i386-linux-gnu/libc.so.6(cfree+0x6d)[0x43d41d]
/usr/lib/i386-linux-gnu/libgfortran.so.3(+0x101ed)[0x9c41ed]
/usr/lib/i386-linux-gnu/libgfortran.so.3(+0xa0a92)[0xa54a92]
/usr/lib/i386-linux-gnu/libgfortran.so.3(+0x9a862)[0xa4e862]
./a.out[0x8048f4b]
./a.out[0x8049198]
/lib/i386-linux-gnu/libc.so.6(__libc_start_main+0xe7)[0x3e3e37]
./a.out[0x8048791]
======= Memory map: ========
003cd000-00527000 r-xp 00000000 08:06 658988     /lib/i386-linux-gnu/libc-2.13.so
00527000-00528000 ---p 0015a000 08:06 658988     /lib/i386-linux-gnu/libc-2.13.so
00528000-0052a000 r--p 0015a000 08:06 658988     /lib/i386-linux-gnu/libc-2.13.so
0052a000-0052b000 rw-p 0015c000 08:06 658988     /lib/i386-linux-gnu/libc-2.13.so
0052b000-0052e000 rw-p 00000000 00:00 0 
008f1000-008f2000 r-xp 00000000 00:00 0          [vdso]
009b4000-00a6c000 r-xp 00000000 08:06 656444     /usr/lib/i386-linux-gnu/libgfortran.so.3.0.0
00a6c000-00a6d000 r--p 000b7000 08:06 656444     /usr/lib/i386-linux-gnu/libgfortran.so.3.0.0
00a6d000-00a6e000 rw-p 000b8000 08:06 656444     /usr/lib/i386-linux-gnu/libgfortran.so.3.0.0
00a6e000-00a6f000 rw-p 00000000 00:00 0 
00b1f000-00b43000 r-xp 00000000 08:06 659006     /lib/i386-linux-gnu/libm-2.13.so
00b43000-00b44000 r--p 00023000 08:06 659006     /lib/i386-linux-gnu/libm-2.13.so
00b44000-00b45000 rw-p 00024000 08:06 659006     /lib/i386-linux-gnu/libm-2.13.so
00b98000-00bb4000 r-xp 00000000 08:06 659007     /lib/i386-linux-gnu/ld-2.13.so
00bb4000-00bb5000 r--p 0001b000 08:06 659007     /lib/i386-linux-gnu/ld-2.13.so
00bb5000-00bb6000 rw-p 0001c000 08:06 659007     /lib/i386-linux-gnu/ld-2.13.so
00d6b000-00d85000 r-xp 00000000 08:06 682158     /lib/i386-linux-gnu/libgcc_s.so.1
00d85000-00d86000 r--p 00019000 08:06 682158     /lib/i386-linux-gnu/libgcc_s.so.1
00d86000-00d87000 rw-p 0001a000 08:06 682158     /lib/i386-linux-gnu/libgcc_s.so.1
08048000-0804a000 r-xp 00000000 08:06 524604     /home/fluidistic/Escritorio/fortran/a.out
0804a000-0804b000 r--p 00001000 08:06 524604     /home/fluidistic/Escritorio/fortran/a.out
0804b000-0804c000 rw-p 00002000 08:06 524604     /home/fluidistic/Escritorio/fortran/a.out
09c55000-09c76000 rw-p 00000000 00:00 0          [heap]
b7600000-b7621000 rw-p 00000000 00:00 0 
b7621000-b7700000 ---p 00000000 00:00 0 
b77c1000-b77c3000 rw-p 00000000 00:00 0 
b77e0000-b77e2000 rw-p 00000000 00:00 0 
bf949000-bf96a000 rw-p 00000000 00:00 0          [stack]
Aborted

Hmm what the heck is happening?!

Edit: I think I found a way not to use the rewind statement. It's when I open a file. I must use the command status="old" so that I'm not rewriting over and over the same file. But I still get the memory error.
However my tabla.dat looks like

 -1.0000000      3.84615399E-02
  -1.0000000      3.84615399E-02
  -1.0000000      3.84615399E-02
  -1.0000000      3.84615399E-02
  -1.0000000      3.84615399E-02
   0.0000000       1.0000000    
   0.0000000       1.0000000    
   0.0000000       1.0000000    
   0.0000000       1.0000000    
   0.0000000       1.0000000    
   1.0000000      3.84615399E-02

if I enter n=10 and m=100.
So basically I'm trying to use a polynomial of degree 10 to interpolate 11 points.
As you can see, my x values are -1, -1, -1, -1, 0, etc. which is very weird to me. I'm going to investigate this. Meanwhile if you have any help to provide, feel free to do it :) Thanks again.
Wow I fixed this problem! I had to put "." after numbers.
So that my code becomes

open(unit=1,file="tabla.dat", status='old')
allocate (x(n+1),y(n+1),f(n+1))
do i=0,n
x(i)=(2.*i)/n-1.
f(i)=1./(1.+25.*x(i)**2.)
write(1,*)x(i),f(i)
end do

Now my tabla.dat is perfect (I've plotted it along with f(x) and the generated points are just perfect).
However I still have this memory error... I really have no idea how to fix it. I really need help on this now.

 

[Useful nucleus]

Let's say you have a fortran code and a file X.dat . Let's say that your code wrote ten lines in the file X.dat. Now let's assume that you ask your code to read. The code will automoatically go to line number 11 to read something but , then it will find it empty and starts complaining! A REWIND statement tells the code to go to the very first line to start reading from there then it will find something to read!
Your REWIND statement shoudl be between the two loops!

NOTE: I do not see why you need to write in a file and then read from it! Why not directly use f(i) as if it is y(i).

 

[fluidistic]

Let's say you have a fortran code and a file X.dat . Let's say that your code wrote ten lines in the file X.dat. Now let's assume that you ask your code to read. The code will automoatically go to line number 11 to read something but , then it will find it empty and starts complaining! A REWIND statement tells the code to go to the very first line to start reading from there then it will find something to read!
Your REWIND statement shoudl be between the two loops!

I just tried the "rewind (unit=1)" between the 2 loops and I still get the memory error message.
Anyway since I write

do i=0, 9
read(1,*) some stuff
end do

I really don't understand why my program will try to read the line 11 if only 10 lines are written and I explicitly tell my program to assign/read values for the 10 first lines.

NOTE: I do not see why you need to write in a file and then read from it! Why not directly use f(i) as if it is y(i).

Because I wanted to plot the data from a separated file in GNUplot. So that I run GNUplot and type "plot tabla.dat" and I have the points plotted.

Note: Could you run the program? I'm using gfortran and I've heard ifortran is "better" and sometimes you can get errors from 1 compiler while the other runs the program fine.

 

[fluidistic]

I really don't understand what's wrong with my read line. I have like 4 other programs with exactly the same read loop line and they work great. I can post one here if you want to look at it. I'm really at a lost about this memory error. :/
EDIT:
I've searched on the internet about this error and typed in a terminal "export MALLOC_CHECK_=0".
Now I don't have the memory error. I do get the error "At line 20 of file test.f90 (unit = 1, file = 'tabla.dat')
Fortran runtime error: End of file" though. I'm trying to fix it but failing as of now.


EDIT2: Some news! I've disabled the part of the program where it writes x(i) and f(i). Instead I've modified the file "tabla.dat" and after executing the program the error of reading disappeared. In other words, Fortran cannot read what it writes; or at least not without telling it explicitly to understand its own notation.
Here is what fortran writes:

-1.0000000      3.84615399E-02
 -0.80000001      5.88235259E-02
 -0.60000002      9.99999940E-02
 -0.40000001      0.19999999    
 -0.20000000      0.50000000    
   0.0000000       1.0000000    
  0.20000000      0.50000000    
  0.40000001      0.19999999    
  0.60000002      9.99999940E-02
  0.80000001      5.88235259E-02
   1.0000000      3.84615399E-02

If I remove the parts "E-02" the program will run without the end of file error.
However I get a Spanish message "Violación de segmento" message. I think it would translate as "Segment fault". Since it's in Spanish I'm guessing linux (ubuntu) has a problem with my program.
If you look at my code/program, I should get some data inside the "result.dat" file. Mine is empty. But at least it was created, a sign my program ran.

Is there any hope to solve my problem with this program?!

EDIT3: I think the problem resides in my subroutine, when I write

if (j/=i) then
  do j=0, n
  l=l*(x_m-x(j))/(x(i)-x(j))
  end do
  end if

since I never defined a value for j before. I don't really know how to proceed.

EDIT4: Nevermind, this isn't the only problem. I've changed the order of the do and if statements... I still get the Segment default error... and my file "result.dat" is still empty.
If someone can help me... here is my code as of now:

program test
implicit none
real:: interval,h,x_m,x_n,p,l
integer:: n,j, i,m
real, dimension(:), allocatable :: x,y,f


write(*,*)"Enter the dimension n of the polynomial and the number of divisions of the interval"
read(*,*)n,m

open(unit=1,file="tabla.dat", status='old')
allocate (x(n+1),y(n+1),f(n+1))
!do i=0,n
!x(i)=(2.*i)/n-1.
!f(i)=1./(1.+25.*x(i)**2.)
!write(1,*)x(i),f(i)
!end do

do i=0,n
read(1,*)x(i),y(i)
end do
close(1)

interval=x(n)-x(0)
h=interval/m
x_m=x(0)
x_n=x(n)

do while (x_m<= x_n)
call Lagrange(x_m,p)
open(unit=2,file="result.dat")
write(2,*)x_m,p
x_m=x_m+h
end do
close(2)


Contains
subroutine Lagrange (x_m,p)
implicit none
real, intent(inout)::p, x_m
integer:: i, j
p=0.

do i=0,n
  do j=0, n
  if (j/=i) then
  l=l*(x_m-x(j))/(x(i)-x(j))
  end if
  end do
p=p+y(i)*l
end do
end subroutine

end program

My file "tabla.dat" is

 -1.0000000       3
 -0.80000001      5
 -0.60000002      9
 -0.40000001      0    
 -0.20000000      0   
   0.0000000       1    
  0.20000000      0 
  0.40000001      0  
  0.60000002      9
  0.80000001      5
   1.0000000      3

And I'm having a Segment fault message when executing the program.
EDIT5 I got it! I was closing (1) too fast. Changed it to under the "close(2)" line. Now I have a result.dat file!

So my only remaining question is, how do I get fortran to understand its own notation?
In other words, I want my program to read and understand the notation of

  -1.0000000      3.84615399E-02
  -1.0000000      3.84615399E-02
  -1.0000000      3.84615399E-02
  -1.0000000      3.84615399E-02
  -1.0000000      3.84615399E-02
   0.0000000       1.0000000    
   0.0000000       1.0000000    
   0.0000000       1.0000000    
   0.0000000       1.0000000    
   0.0000000       1.0000000    
   1.0000000      3.84615399E-02

How can I do that?! I'm so excited to having found out the problems!

 

[fluidistic]

I solved the problem. For some reasons I wasn't able to make the rewind command work properly.
What I did was close(1) between the 2 loops and then I opened (1) again and closed it after the second loop.
I had made a small mistake in the algorithm too, now all is fixed and thus the thread is solved.
Thank you for having helped me!

 

[Pinčiukas]

fortran arrays begin with x(1), not x(0), so i guess there could be more errors.

in this section you're accessing wrong array elements:
allocate (x(n+1),y(n+1),f(n+1))
<...>
do i=0,n
read(1,*)x(i),y(i)
end do

it can be corrected:
allocate (x(n+1)) -> allocate (x(0:n+1))

well, I could be wrong in case everything works fine already:)

 

[fluidistic]

fortran arrays begin with x(1), not x(0), so i guess there could be more errors.

in this section you're accessing wrong array elements:
allocate (x(n+1),y(n+1),f(n+1))
<...>
do i=0,n
read(1,*)x(i),y(i)
end do

it can be corrected:
allocate (x(n+1)) -> allocate (x(0:n+1))

well, I could be wrong in case everything works fine already:)

You are right, I forgot to mention I changed this part also. I tried with x(n+1) and it didn't work. Instead, x(0:n) worked.