1.

program autoc
! autocorelation function for random numbers

  use plot
  implicit none

  integer, parameter :: n=100
  real, dimension(maxpoint) :: x, t, R

  integer :: i, k

  call random_number(x)

  do k=1,n/2
     t(k)=k
     R(k)= ac(x, n, k)
     write(*,*) k,R(k)
  end do

  ! plot R(k)
  call init("yy.ps", 3.0, 10.0, cm, cm)
  call frame(0.0, 0.0, 10.0, 5.0)
  call plotcurve(0.1*t, 100*R, n/2, 0.5, 0)
  call finish()

contains

real function ac(x, n, k)
! autocorrelation for delay k
! x(1:n)=observed values
  integer, intent(in) :: n, k
  real, dimension(n), intent(in) :: x
  real :: s
  integer :: i

  s = 0.0
  do i=1,n/2
     s = x(i)*x(i+k)
  end do

  ac = s/(n/2-1)
end function

end program

           1  1.05944360E-02
           2  4.67108050E-03
           3  1.59435032E-03
           4  8.45734589E-03
           5  1.12965358E-02
           6  1.30376127E-02
           7  1.00646457E-02
           8  6.84996694E-03
           9  5.15948376E-03
          10  5.80240833E-03
          11  7.61120906E-03
          12  1.37372604E-02
          13  1.36336777E-02
          14  1.02736261E-02
          15  1.31293545E-02
          16  1.28400978E-03
          17  1.01044979E-02
          18  1.03486776E-02
          19  1.30342245E-02
          20  9.72115528E-03
          21  1.12028243E-02
          22  7.68955750E-03
          23  8.49432312E-04
          24  6.61287503E-03
          25  8.22773948E-03
          26  1.89325050E-03
          27  8.08602478E-03
          28  7.15783192E-03
          29  1.21949157E-02
          30  4.18222090E-03
          31  9.21843760E-03
          32  9.15350020E-03
          33  6.93356758E-03
          34  3.60081764E-03
          35  1.05391047E-03
          36  1.39379269E-03
          37  7.56113790E-03
          38  5.17026568E-03
          39  2.08546204E-04
          40  1.09151965E-02
          41  8.54693912E-03
          42  1.06493514E-02
          43  1.31269004E-02
          44  1.57267659E-03
          45  4.38392535E-03
          46  8.21576361E-03
          47  6.62867678E-04
          48  1.57214561E-03
          49  2.97295186E-03
          50  1.38448318E-03

If you plot this you'll see that there is no obvious correlation


2.

program struc
! structure function for random numbers

  use plot
  implicit none

  integer, parameter :: n=100
  real, dimension(maxpoint) :: t, x, R
  integer :: i, k

  call random_number(x)

  do k=1,n/2
     t(k)= k
     R(k)= strf(x, n, k)
     write(*,*) k,R(k)
  end do

  ! plot R(k)
  call init("yy.ps", 3.0, 10.0, cm, cm)
  call frame(0.0, 0.0, 10.0, 5.0)
  call plotcurve(0.1*t, 10*R, n/2, 0.5, 0)
  call finish()

contains

real function strf(x, n, k)
! structure function for delay k
! x(1:n)=observed values
  integer, intent(in) :: n, k
  real, dimension(n), intent(in) ::  x

  real :: s
  integer :: i

  s=0
  do i=1,n-k
     s = s + (x(i+k)-x(i))**2
  end do
  strf = s/(n-k)

end function
end program

           1  0.13625167    
           2  0.15573707    
           3  0.14728901    
           4  0.18636374    
           5  0.16529454    
           6  0.14978512    
           7  0.17042172    
           8  0.15879260    
           9  0.17032824    
          10  0.16005220    
          11  0.15402590    
          12  0.12056794    
          13  0.13815081    
          14  0.14472766    
          15  0.17639984    
          16  0.19760177    
          17  0.16122319    
          18  0.15341517    
          19  0.15879562    
          20  0.16402225    
          21  0.16452965    
          22  0.17126027    
          23  0.16306941    
          24  0.13441356    
          25  0.15927915    
          26  0.19143927    
          27  0.16602427    
          28  0.19553414    
          29  0.19642645    
          30  0.18198365    
          31  0.15754066    
          32  0.18063410    
          33  0.19244763    
          34  0.19912025    
          35  0.18675356    
          36  0.19061990    
          37  0.15802692    
          38  0.17146778    
          39  0.17659733    
          40  0.19796273    
          41  0.16554521    
          42  0.17762996    
          43  0.19682176    
          44  0.16874719    
          45  0.15308674    
          46  0.18712041    
          47  0.16846226    
          48  0.20016560    
          49  0.17400251    
          50  0.17879173    

Again, no correlation

3.

program forktest
! find the minimum of a 2D function by forking
! this is not a very decent solution!
  implicit none
  real, parameter :: limit=0.0001
  integer, parameter :: maxcount=1000
  real :: x0=1.0, y0=2.0, x1=100, y1=100, d=5.0, diff
  integer :: i=0

  x0 = forkx(x0, y0)
  write(*,*) i, x0, y0, f(x0, y0)

  do 
    y0 = forky(x0, y0)
    x0 = forkx(x0, y0)

    write(*,*) i, x0, y0, f(x0, y0)
    diff = sqrt((x1-x0)**2 + (y1-y0)**2)
    if (diff < limit) exit
    i = i+1
    if (i > maxcount) exit
    x1=x0; y1=y0
  end do

contains

real function f(x, y)
! Rosenbrock banana function
  real, intent(in) :: x, y
  f = 100*(y-x**2)**2 + (1-x)**2
end function

real function dfdx(x, y)
  real, intent(in) :: x, y
  dfdx = -400*x*(y-x**2) - 2*(1-x)
end function

real function dfdy(x, y)
  real, intent(in) :: x, y
  dfdy = 200*(y-x**2)
end function

real function forkx(x0, yy)
! minimize f(x,y) for a fixed value y=yy
! proceed in the direction of the x axis
  real, intent(in) :: x0, yy
  real :: a, b, c, fa, fb, fc, df

  a=x0; df = dfdx(x0, yy)
  if (df > 0) then
     b = a - d
  else
     b = a + d
  end if
  c = (a+b)/2

  do while (abs(a-c)>limit)
    df = dfdx(c, yy)
    if (df > 0) then
       b = c
    else
       a = c
    end if
    c = (a+b)/2

!  write(*,*) a,b,c, f(c, yy)
  end do
  forkx=c
end function


real function forky(xx, y0)
! minimize f(x,y) for a fixed value x=xx
  real, intent(in) :: y0, xx
  real :: a, b, c, fa, fb, fc, df

  a=y0; df = dfdy(xx, y0)
  if (df > 0) then
     b = a - d
  else
     b = a + d
  end if
  c = (a+b)/2

  do while (abs(a-c)>limit)
    df = dfdy(xx, c)
    if (df > 0) then
       b = c
    else
       a = c
    end if
    c = (a+b)/2
!  write(*,*) a,b,c, f(xx, c)
  end do
  forky=c

end function

end program


           0   1.4137421       2.0000000      0.17136028    
           0  -3.5861816      -2.9999237       25176.967    
           1  -1.4111176       2.0000000       5.8211393    
           2  1.61743164E-03  -2.9999237       900.95264    
           3  0.16130066      3.05175781E-05  0.77095103    
           4  0.21142578      2.60467529E-02  0.65664685    
           5  0.24522400      4.47387695E-02  0.59339064    
           6  0.27108765      6.00738525E-02  0.54930854    
           7  0.29237366      7.34252930E-02  0.51527232    
           8  0.31091309      8.55560303E-02  0.48718601    
           9  0.32701111      9.66186523E-02  0.46355939    
          10  0.34143066      0.10691833      0.44303846    
     ...
         995  0.97978210      0.95980835      4.11473535E-04
         996  0.97985840      0.96003723      4.06410603E-04
         997  0.97993469      0.96011353      4.05128434E-04
         998  0.98001099      0.96034241      4.00187215E-04
         999  0.98008728      0.96041870      3.98837437E-04
        1000  0.98016357      0.96049500      3.98576172E-04

Getting close, but the convergence is really slow



4.
v = can be vectorised
p = can be parallelised

module linear
...
implicit none

contains

subroutine lu (A, order, n)
...

  integer, intent(in) :: n
  real, dimension(n,n), intent(inout) :: A
  integer, dimension(n), intent(out) :: order

  real, dimension(n) :: tmp
  real :: s, maxv
  integer i,j,k,l, m, maxind, t

  ! original order of the lines
  do i=1,n              ! v p
     order(i) = i
  end do

  ! first line
  ! pivoting, find the largest element on the 1. column
  maxind = 1
  maxv = A(1,1)
  do l=2,n
     if (A(l, 1) > maxv) then
        maxv = A(l, 1) ; maxind = l
     end if
  end do
  ! exchange the line found with the first line
  tmp=A(maxind,:)
  A(maxind,:) = A(1,:)
  A(1,:)=tmp
  ! keep track of the order of the lines
  t = order(maxind)
  order(maxind) = 1
  order(1) = t

  do j=2,n                  ! v p
     A(1,j)=A(1,j)/A(1,1)
  end do

  ! other lines
  do m=2,n
     ! first the column
     ! pivoting, find the largest element below A(m,m)
     maxind = m
     maxv = A(m,m)
     do l=m+1,n
       if (A(l, m) > maxv) then
          maxv = A(l, m) ; maxind = l
       end if
     end do
     ! exchange lines
     tmp=A(maxind,:)
     A(maxind,:) = A(m,:)
     A(m,:)=tmp
     t = order(maxind)
     order(maxind) = m
     order(m) = t

     ! handle the column below A(m,m)
     do i=m,n
        s=0.0
        do k=1,m-1            ! v
           s=s+A(i,k)*A(k,m)
        end do
        A(i,m)=A(i,m)-s
     end do

     ! then the line to the right of A(m,m)
     do j=m+1,n
        s=0.0
        do k=1,m-1            ! v
           s=s+A(m,k)*A(k,j)
        end do
        A(m,j)=(A(m,j)-s)/A(m,m)
     end do
   end do
end subroutine

subroutine solve (A, b, order, n, x)
! Solve the equation Ax=b.
! A(n, n) contains the LU decomposition of the coefficient matrix,
! b is the right-hand side of the equation.
! order is the vector returned by the subroutine lu, giving the
! permuted order of the lines of the matrix.
! The solution will be returned in the vector x(n), 
! the order of which will correspond to the original matrix

  integer, intent(in) :: n
  real, dimension(n,n), intent(in) :: A
  integer, dimension(n), intent(in) :: order
  real, dimension(n), intent(in) :: b
  real, dimension(n) :: x
  integer i,k
  real s    

  ! permute the right-hand side to correspond to the coefficien matrix
  do i=1,n               ! v p
     x(order(i)) = b(i)
  end do

  x(1) = x(1)/A(1,1)

  do i=2,n
     s=0.0
     do k=1,i-1           ! v
        s=s+A(i,k)*x(k)
     end do
     x(i) = (x(i)-s) / A(i,i)
  end do  

  do i=n-1,1,-1
     s=x(i)
     do k=i+1,n         ! v
        s=s-A(i,k)*x(k)
     end do
     x(i)=s
  end do

end subroutine

function inverse (A, order, n)
! Find the inverse of the matrix A(n,n).
! Before calling the function A must be replaced by its LU decomposition.
! order is the vector returned by the subroutine lu, giving the
! permuted order of the lines of the matrix.

  integer, intent(in) :: n
  real, dimension(n,n), intent(in) :: A
  integer, dimension(n), intent(in) :: order
  real, dimension(n,n) :: inverse, C, S
  integer i

  C=0.0
  do i=1,n         ! v p
     C(i,i)=1.0
  end do

  do i=1,n
     call solve(A, C(:,i), order, n, S(:,i))
  end do

  inverse = S

end function

function reorder (A, order, n)
! permute the lines of the matrix A(n,n) in the order 
! given in the vector order

  real, dimension(n,n) :: reorder
  integer, intent(in) :: n
  real, dimension(n,n), intent(in) :: A
  integer, dimension(n), intent(in) :: order
  real, dimension(n,n) :: C
  integer i

  do i=1,n                ! v p
     C(order(i),:)=A(i,:)
  end do

  reorder = C
end function

subroutine dump(A, n, s)
! Write the matrix A(n,n)
! The string s is the title text

  integer, intent(in) :: n
  real, intent(in), dimension(n,n) :: A
  character(len=*) :: s
  integer i

  write(*,'(/A)') s
  do i=1,n
     write(*,'(10F6.2)') A(i,:)
  end do

end subroutine

end module