cmrh_NA/Fortran/helsing.f90

program helsing
! 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
! Author: Sebastien 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

open(unit=20,file='inputfile.dat',status='old')
read(20,*)init
n=init(1)
close(20)
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 program helsing