Fortran Solving Function Problem in Code Exchange

[Galizius]

Hello I am trying to exchange following code into two functions which I will need to use in further programs.
The code is as follows:

do i=0,k
  do j=0,k
  if (i/=j) then
    lj(i)=lj(i)*((xx-x(i))/(x(j)-x(i)))
  end if
  enddo
  l(i)=l(i)+y(i)*lj(i)
  write(2,*) x(i),y(i),l(i)
enddo

I would like to swap lj and l arrays for the functions. I have tried something like that:

real(8) function lj(j,xx,x,k)
  real(8) :: x(1000)
  real(8) :: xx
  integer :: i,j,k

  lj=1.d0
do i=0,k
  do j=0,k
  if (i/=j) then
    lj=lj*((xx-x(i))/(x(j)-x(i)))
  end if
  enddo
enddo

real(8) function l(xx,x,y,k)
real(8) :: x(1000),y(1000), xx
integer :: i,j,k
  do i=0,k
    l=l+y(i)*lj(i,xx,x,k)
  enddo

And then I am trying to call it in the main program

do i=0,k
  write(2,*) x(i),y(i),l(x(i),x,y,k)
enddo

But it does not seem to be working as before. What am I doing wrong?

 

[DrClaude]

But it does not seem to be working as before.

This is very vague. You'll get better help when you clearly describe why (you think) it is wrong.

real(8) function l(xx,x,y,k)
real(8) :: x(1000),y(1000), xx
integer :: i,j,k
  do i=0,k
    l=l+y(i)*lj(i,xx,x,k)
  enddo

l is not initialized before being used in the loop.

 

[BvU]

Hello I am trying to exchange following code into two functions which I will need to use in further programs.
What am I doing wrong?

Not sure. In my compiler arrays are 1-based and dividing by 0 is disastrous. From your snippets I can come to something that runs but I have no idea if it's what you intend:

real(8) function lj(j,xx,x,k)
  real(8) :: x(1000)
  real(8) :: xx
  integer :: i,j,k
 
  lj=1.d0
 do i=1,k
  do j=1,k
  if (x(i) /= x(j)) then
  if (i/=j) then
  lj=lj*((xx-x(i))/(x(j)-x(i)))
  end if
  end if
  enddo
 enddo
end
 real(8) function l(xx,x,y,k)
 real(8) :: x(*),y(*), xx,  lj
 external lj
 integer :: i,j,k
  do i=1,k
  l=l+y(i)*lj(i,xx,x,k)
  enddo
 end
 
Program test
real*8 l,x(1000),y(1000)
external l
k=4
do i=1,k
  write(2,*) x(i),y(i),l(x(i),x,y,k)
pause
enddo
end

with output

  0.000000000000000E+000  0.000000000000000E+000  6.790386531088871E-313

suspicious ! We should probably initialize a few things

 

[DrClaude]

I just noticed another problem:

real(8) function lj(j,xx,x,k)
  real(8) :: x(1000)
  real(8) :: xx
  integer :: i,j,k

  lj=1.d0
do i=0,k
  do j=0,k
  if (i/=j) then
    lj=lj*((xx-x(i))/(x(j)-x(i)))
  end if
  enddo
enddo

real(8) function l(xx,x,y,k)
real(8) :: x(1000),y(1000), xx
integer :: i,j,k
  do i=0,k
    l=l+y(i)*lj(i,xx,x,k)
  enddo

The code is accessing x(0) and y(0), which are outside the array boundaries.

 

[DrClaude]

Not sure. In my compiler arrays are 1-based and dividing by 0 is disastrous.

I guess you didn't mean dividing by 0, in which case you beat me to it!

 

[Galizius]

I just noticed another problem:The code is accessing x(0) and y(0), which are outside the array boundaries.

That solved my problem! Thank You so much, I do not know how could I miss that...

 

[DrClaude]

That solved my problem! Thank You so much, I do not know how could I miss that...

Be sure to also implement the fix to the problem I mentioned in post #2. It might work for now, but you never know when it might fail.

 

[Galizius]

Yes, I already did. Thank You again!

 

[Mark44]

real(8) function lj(j,xx,x,k)
  real(8) :: x(1000)
  real(8) :: xx
  integer :: i,j,k
  lj=1.d0
do i=0,k
  do j=0,k
  if (i/=j) then
    lj=lj*((xx-x(i))/(x(j)-x(i)))
  end if
  enddo
enddo

A few comments about your code:

1. The names you have picked for most of your variables are terrible! Single letter names such as i and j are OK for loop control variables, but names like lj, x, xx, and k don't give any indication of what you're going to use them for. There is no reason not to choose a name that suggest what the variable represents. Many years ago there were limits on how long a name you could use for a variable, so having variable names like x, xx, and lj is completely unwarranted, and makes your code harder to read than it should be.
2. Why are you passing j as a parameter of your function? This variable gets set to 0, 1, 2, ..., k. It is strictly a local loop control variable, so there is no point in passing it as a function parameter.
3. There isn't a single comment in your code.
4. There isn't a single whitespace character in your code. The compiler doesn't care, but it's much harder for humans to read this -- lj=lj*((xx-x(i))/(x(j)-x(i))) -- than this --
   lj = lj * ( (xx - x(i) ) / ( x(j) - x(i) ) )