ajout commentaires, authors ( fortran ) , mise en page Readme.txt

Fortran/CMRH.f90

@@ -1,6 +1,9 @@
 program CMRH
 !***************************************************************
 ! CMRH program
+!
+! Created by H. Sadok and M. Heyouni -1999
+! Updated by S. Duminil dec. 2012
 !***************************************************************
 
 external SSWAP,DCOPY,DSCAL,DGEMV,DSWAP

Fortran/PCMRH.f90

@@ -1,6 +1,8 @@
 program PCMRH
 !***************************************************************
 ! Parallel CMRH program with MPI
+!
+! Created by S. Duminil sep. 2012 
 !***************************************************************
 use mpi
 external SSWAP,DCOPY,DSCAL,DGEMV,DSWAP
@@ -36,7 +38,6 @@ call MPI_COMM_RANK(MPI_COMM_WORLD,rang,code)
 call MPI_COMM_SIZE(MPI_COMM_WORLD,Nprocs,code)
 call MPI_BARRIER(MPI_COMM_WORLD,code)
 
-
 call itime(timearray)     ! Get the current time
 i = rand ( timearray(1)+timearray(2)+timearray(3) )
 
@@ -91,8 +92,6 @@ endif
 call MPI_BCAST(solchoice,1,MPI_INTEGER,0,MPI_COMM_WORLD,code)
 call solution(AL,workl,sol,nl,n,solchoice)
 
-
-
 call MPI_ALLGATHER(workl,nl,MPI_DOUBLE_PRECISION,b,nl,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,code)
 
 !********************************************
@@ -132,32 +131,32 @@ D_inv  = 1.d0/beta
 ! First permutation
 !********************************************
 if (indice .ne.1) then
-call dswap(1,b(indice),1,b(1),1)
-call dswap(nl,AL(1,indice),1,AL(1,1),1)
-call position(indice,nl,inl,inlp)
-if (inl==0) then
-if (rang==inl)then
-call dswap(n,AL(inlp,1),nl,AL(1,1),nl)
-end if
-else
-if (rang==inl) then
-call dcopy(n,AL(inlp,1),nl,temp,1)
-end if
-call MPI_BCAST(temp,n,MPI_DOUBLE_PRECISION,inl,MPI_COMM_WORLD,code)
-if (rang==0)then
-call dcopy(n,AL(1,1),nl,temp2,1)
-end if
-call MPI_BCAST(temp2,n,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,code)
-if (rang==inl)then
-call dcopy(n,temp2,1,AL(inlp,1),nl)
-end if
-if (rang==0)then
-call dcopy(n,temp,1,AL(1,1),nl)
-end if
-end if
-if (rang==0)then
-call sswap(1,p(indice),1,p(1),1)
-end if
+   call dswap(1,b(indice),1,b(1),1)
+   call dswap(nl,AL(1,indice),1,AL(1,1),1)
+   call position(indice,nl,inl,inlp)
+   if (inl==0) then
+      if (rang==inl)then
+         call dswap(n,AL(inlp,1),nl,AL(1,1),nl)
+      end if
+   else
+      if (rang==inl) then
+         call dcopy(n,AL(inlp,1),nl,temp,1)
+      end if
+      call MPI_BCAST(temp,n,MPI_DOUBLE_PRECISION,inl,MPI_COMM_WORLD,code)
+      if (rang==0)then
+         call dcopy(n,AL(1,1),nl,temp2,1)
+      end if
+      call MPI_BCAST(temp2,n,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,code)
+      if (rang==inl)then
+         call dcopy(n,temp2,1,AL(inlp,1),nl)
+      end if
+      if (rang==0)then
+         call dcopy(n,temp,1,AL(1,1),nl)
+      end if
+   end if
+   if (rang==0)then
+      call sswap(1,p(indice),1,p(1),1)
+   end if
 end if
 
 !********************************************
@@ -225,14 +224,14 @@ do while((dnorm.gt.tol).and.(k.lt.n))
         end if
     end if
 
-    !********************************************
-    !   j Loop
-    !********************************************
+!********************************************
+!   j Loop
+!********************************************
     do j=1,k
-        !********************************************
-        !   A(j,k)=workl(j)
-        !   workl(j)=0
-        !********************************************
+!********************************************
+!   A(j,k)=workl(j)
+!   workl(j)=0
+!********************************************
         call position(j,nl,knl,knlp)
         if (rang==knl) then
             AL(knlp,k)=workl(knlp)
@@ -240,9 +239,9 @@ do while((dnorm.gt.tol).and.(k.lt.n))
             alpha=-AL(knlp,k)
         end if
         call MPI_BCAST(alpha,1,MPI_DOUBLE_PRECISION,knl,MPI_COMM_WORLD,code)
-        !********************************************
-        !   workl=workl - alpha A(:,j)
-        !********************************************
+!********************************************
+!   workl=workl - alpha A(:,j)
+!********************************************
         call position(j+1,nl,knl,knlp)
         if (rang >knl)then
             call daxpy(nl,alpha,AL(1,j),1,workl,1)
@@ -251,19 +250,19 @@ do while((dnorm.gt.tol).and.(k.lt.n))
             call daxpy(nl-knlp+1,alpha,AL(knlp,j),1,workl(knlp),1)
         end if
     end do
-    !********************************************
-    ! End j Loop
-    !********************************************
+!********************************************
+! End j Loop
+!********************************************
     if (k.lt.n)then
         call MPI_ALLGATHER(workl,nl,MPI_DOUBLE_PRECISION,b,nl,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,code)
-        !********************************************
-        !   Search indice i_0
-        !********************************************
+!********************************************
+!   Search indice i_0
+!********************************************
         indice=k+idamax(n_k,b(k1),1)
         H_k1=b(indice)
-        !********************************************
-        !   Permutation
-        !********************************************
+!********************************************
+!   Permutation
+!********************************************
         if (indice.ne.k1) then
             if (rang==0) then
                 call sswap(1,p(indice),1,p(k1),1)
@@ -294,13 +293,13 @@ do while((dnorm.gt.tol).and.(k.lt.n))
                 end if
             end if
         end if
-        !********************************************
-        !   End Permutation
-        !********************************************
+!********************************************
+!   End Permutation
+!********************************************
         if (H_k1.ne.0.d0)then
-            !********************************************
-            !   b=b/(b)_(i_0)
-            !********************************************
+!********************************************
+!   b=b/(b)_(i_0)
+!********************************************
             D_inv=1.d0/H_k1
             ! call dcopy(n,work,1,b,1)
             call dscal(n,D_inv,b,1)
@@ -309,9 +308,9 @@ do while((dnorm.gt.tol).and.(k.lt.n))
         H_k1=0.d0
         call dscal(n,0.d0,b,1)
     end if
-    !********************************************
-    !   Application of Givens rotations
-    !********************************************
+!********************************************
+!   Application of Givens rotations
+!********************************************
     if (k.gt.1)then
         do i=2,k
             call position(i-1,nl,knl,knlp)
@@ -357,9 +356,9 @@ do while((dnorm.gt.tol).and.(k.lt.n))
             AL(knlp,k)=c(k)*btemp+s(k)*H_k1
         end if
     end if
-    !********************************************
-    ! End Givens rotations
-    !********************************************
+!********************************************
+! End Givens rotations
+!********************************************
 end do
 
 !********************************************
@@ -534,23 +533,23 @@ end subroutine matrice
 ! Solution Choice
 !********************************************
 subroutine solution(A,b,sol,nl,n,solchoice)
-integer    ::n,i,nl,solchoice
-real(kind=8),dimension(n,n) ::A
-real(kind=8),dimension(n)   ::sol
-real(kind=8),dimension(nl)  ::b
-real   ::rand
-do i=1,n
-if (solchoice==1) then
-sol(i)=1.d0
-end if
-if (solchoice==2) then
-sol(i)=dble(n-i+1)
-end if
-if (solchoice==3)then
-sol(i)=rand(0)
-end if
-end do
-call dgemv('N',nl,n,1.d0,A,nl,sol,1,0.d0,b,1)
+  integer    ::n,i,nl,solchoice
+  real(kind=8),dimension(n,n) ::A
+  real(kind=8),dimension(n)   ::sol
+  real(kind=8),dimension(nl)  ::b
+  real   ::rand
+  do i=1,n
+     if (solchoice==1) then
+        sol(i)=1.d0
+     end if
+     if (solchoice==2) then
+        sol(i)=dble(n-i+1)
+     end if
+     if (solchoice==3)then
+        sol(i)=rand(0)
+     end if
+  end do
+  call dgemv('N',nl,n,1.d0,A,nl,sol,1,0.d0,b,1)
 end subroutine solution
 
 
@@ -562,15 +561,3 @@ integer :: k,nl,knl,knlp
 knl=(k-1)/nl
 knlp=k-knl*nl
 end subroutine position
-
-
-
-
-
-
-
-
-
-
-
-

Fortran/helsing.f90

@@ -3,13 +3,16 @@ program helsing
 ! 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
+! 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

Fortran/initialisation.f90

@@ -1,6 +1,12 @@
 program initialisation
-! This program creates a file inputfile.dat
-! 
+! This program creates the file inputfile.dat which is used in 
+! CMRH (sequential) or PCMRH (parallel) test suite.
+!
+! - output choice (no prints/ all iterations & results / only final results) 
+! - size of matrix A.
+! - choice between 5 types of matrices A
+! - choice between 3 types of solution vectors
+! - tolerance input.
 
 use randvectors
 integer :: iprint,n,choix,choixsol
@@ -15,8 +21,8 @@ print*,'          2i-1 / n-i+j for j=i+1,n '
 print*,'2- A(i,j)=i for j=1,i '
 print*,'          j for j=i+1,n '
 print*,'3-helsing.f90'
-print*,'4-rbf.f90'
-print*,'5- Matrix File (copy in matrixfile.dat) '
+print*,'4-rbf.f90 (Radial basis function matrix - see rbf.f90 - )'
+print*,'5- Matrix File (put your matrixfile.dat in current directory) '
 read*,choix
 if (choix==3) then
     call rand0(n)
@@ -34,7 +40,7 @@ print*,'Enter solution choice : 1- [1,...,1]^T  / 2- [n,...,1]^T '
 print*,'3-rand vector'
 read*,choixsol
 
-print*,'Enter tol :'
+print*,'Enter tolerance:'
 read*,choixtol
 
 open(unit=10,file='inputfile.dat',status='unknown')

Fortran/rand0.f90

@@ -1,7 +1,7 @@
 program rand0
 ! rand0.f90
 ! This file creates a random vector of size n and stores it in rand0.dat
-
+! Author: Sebastien Duminil
 
   implicit none
   integer n,i,init(4)

Fortran/randvectors.f90

@@ -2,7 +2,7 @@ module randvectors
 ! This file contains some subroutines : 
 ! rand0, vecrand
 ! helsing, rbf
-
+! Author: S. Duminil
 implicit none
 contains
 
@@ -67,9 +67,13 @@ subroutine helsing(n)
 ! 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
+! 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
@@ -85,8 +89,8 @@ 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)
+   xx(i)=value(2*i-1)
+   yy(i)=value(2*i)
 end do
 close(30)
 open(unit=40,file='rand0.dat',status='old')
@@ -94,23 +98,23 @@ 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
+   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)
+   write(10,*) (A(i,j),j=1,n)
 end do
 close(10)
 deallocate(A)
@@ -137,51 +141,45 @@ 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)
+   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
+   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
+   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
+   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)
+   write(10,*) (A(i,j),j=1,n)
 end do
 close(10)
 deallocate(A)
 end subroutine rbf
 
-
-
-
-
-
-
 end module randvectors
 

Fortran/rbf.f90

@@ -4,7 +4,7 @@ program rbf
 ! 
 !
 ! We store this matrix in matrixfile.dat
-! 
+! Author: Sebastien Duminil
 
 integer :: m,n,d,i,j,init(4)
 real(kind=8) :: l,r,coef

Readme.txt

@@ -37,10 +37,10 @@
 	- A parallel implementation of the CMRH method for dense linear systems, 
 	  Numer. Algorithms, 63 (2013), 127-142
 
-	SOFWARE REVISION DATE:
+      SOFWARE REVISION DATE:
 	  v1.0, January 2013
 
-	SOFTWARE LANGUAGE:
+      SOFTWARE LANGUAGE:
 	  Fortran 90 and MATLAB 8.0 (R2012b)
 
 =============================================================================