module randvectors
! This file contains some subroutines : 
! rand0, vecrand
! helsing, rbf
! Author: S. Duminil
implicit none
contains

subroutine rand0(n)
integer n,i,init(4)
real(kind=8),allocatable,dimension(:) ::x
integer*4,dimension(3)  ::timearray
real                     ::rand
call itime(timearray)     ! Get the current time
i = rand ( timearray(1)+timearray(2)+timearray(3) )
close(20)
allocate(x(n))
do i=1,n
x(i)=rand(0);
enddo
open(unit=10,file='rand0.dat',status='unknown')
do i=1,n
write(unit=10,FMT=*) x(i)
enddo
close(10)
deallocate(x)
end subroutine rand0

subroutine vecrand(n)
! vecrand.f90
! This file creates two random vectors and stores them in the file vecrand.dat
!
! 01 February 13
! Created by S.Duminil
integer n,i,j,a,b,init(4)
real(kind=8),allocatable,dimension(:) ::x,y
integer*4,dimension(3)  ::timearray
real                     ::rand
call itime(timearray)     ! Get the current time
i = rand ( timearray(1)+timearray(2)+timearray(3) )

a=1
b=0
do while(a.ge.b)
print*,'Enter the interval limits'
read*,a,b
if (a.ge.b) then
print*,'Choose a<b'
endif
end do
allocate(x(n),y(n))
do i=1,n
x(i)=rand(0)*dble(b-a)-dble(b);
y(i)=rand(0)*dble(b-a)-dble(b);
enddo
open(unit=10,file='vecrand.dat',status='unknown')
do i=1,n
write(unit=10,FMT=*) x(i),y(i)
enddo
close(10)
deallocate(x,y)
end subroutine vecrand


subroutine helsing(n)
! helsing.f90
! This file creates a matrix defined by :
! A(i,j)= -log|z(i)-z(j)| , if i .neq. j
!         -log|r(i)|
! where z(i) are n somehow randomly distributed points in a unit square
! centered at the origin in the complex plane and where each r(i) is a 
! number in (0,d(i)[, d(i) being the distance between the point z(i) 
! and its nearest neighbour.
! For more details, see:
! S. Duminil, A parallel implementation of the CMRH method for dense
! linear systems, Numer. Algor., DOI : 10.1007/s11075-012-9616-4
!
! We store this matrix in matrixfile.dat
! 01 February 2013
! S. Duminil

integer :: m,n,coef,d,i,j,init(4)
real(kind=8) :: l,r
real(kind=8), allocatable, dimension(:) :: xx,yy,value,value2
real(kind=8), allocatable, dimension(:,:) :: A

m=n
allocate(value(2*m),value2(m),xx(m),yy(m),A(n,n))
open(unit=30,file='vecrand.dat',status='old')
read(30,*)value
do i=1,m
   xx(i)=value(2*i-1)
   yy(i)=value(2*i)
end do
close(30)
open(unit=40,file='rand0.dat',status='old')
read(40,*)value2
close(40)
coef=-1
do i=1,n
   d=0
   do j=1,m
      r=dsqrt((xx(i)-xx(j))**2+(yy(j)-yy(i))**2)
      if (r>0) then
         A(i,j)=coef*log(r)
         if (d.eq.0)then
            d=j
            l=dble(0.5*value2(i)*dble(d))
            A(i,i)=-log(l)
         end if
      end if
   end do
end do
deallocate(xx,yy,value,value2)
open(unit=10,file='matrixfile.dat',status='unknown')
do i=1,n
   write(10,*) (A(i,j),j=1,n)
end do
close(10)
deallocate(A)
end subroutine helsing


subroutine rbf(n)
! rbf.f90
! This file creates a matrix defined by :
!
!
! We store this matrix in matrixfile.dat
! 01 February 2013
! S. Duminil

integer :: m,n,d,i,j,init(4)
real(kind=8) :: l,r,coef
real(kind=8), allocatable, dimension(:) :: xx,yy,value,value2
real(kind=8), allocatable, dimension(:,:) :: A
real :: rand

m=n-3
allocate(value(2*m),value2(m),xx(m),yy(m),A(n,n))
open(unit=30,file='vecrand.dat',status='old')
read(30,*)value
do i=1,m
   xx(i)=value(2*i-1)
   yy(i)=value(2*i)
end do
close(30)
do i=1,m
   d=0
   do j=1,m
      r=dsqrt((xx(i)-xx(j))**2+(yy(j)-yy(i))**2)
      if (r>0) then
         A(i,j)=r*r*log(r)
         if (d==0)then
            d=j
            coef=0.5*d*rand(0)
            A(i,i)=coef*coef*log(coef)
         end if
      end if
   end do
end do
do j=1,m
   A(m+1,j)=xx(j)
   A(m+2,j)=yy(j)
   A(m+3,j)=1.d0
   A(j,m+1)=xx(j)
   A(j,m+2)=yy(j)
   A(j,m+3)=1.d0
end do
do i=m+1,n
   do j=m+1,n
      A(i,j)=0.d0
   end do
end do
deallocate(xx,yy,value,value2)
open(unit=10,file='matrixfile.dat',status='unknown')
do i=1,n
   write(10,*) (A(i,j),j=1,n)
end do
close(10)
deallocate(A)
end subroutine rbf

end module randvectors