Initial commit

Fortran/CMRH.f90

@@ -0,0 +1,429 @@
+program CMRH
+!***************************************************************
+! CMRH program
+!***************************************************************
+
+external SSWAP,DCOPY,DSCAL,DGEMV,DSWAP
+intrinsic dabs,dsqrt
+!********************************************
+! VARIABLES
+!********************************************
+integer             ::n,code,i,j,k,choix,solchoice,iprint
+real(kind=8) ::init(5)
+real(kind=8),allocatable,dimension(:,:)  ::A
+real(kind=8)                           ::btemp,bmax,btemp2,err,res
+real(kind=8),allocatable,dimension(:)   ::b,sol,rs,temp,x,work,tempn,c,s,temp2,xx,yy,valeurs
+integer,allocatable,dimension(:) ::p
+real*8                          ::dnorm2,D_inv,beta,dnorm,dnrm2,tol,alpha,H_k1,dgam,r
+integer                         ::indice,idamax,k1,n_k
+real dtime, delta, delta0, da(2), dt1, dt2, dt3
+real etime, eelta, eelta0, ea(2), et1, et2, et3
+real t1, t2, t3, tm1, tm2, tm3
+real dtime1, delta1, delta01, da1(2)
+real etime1, eelta1, eelta01, ea1(2),l,p1,p2
+integer nb_ligne
+real                     ::rand
+da(1)=0.d0
+da(2)=1.d0
+
+
+!call itime(timearray)     ! Get the current time
+!i = rand ( timearray(1)+timearray(2)+timearray(3) )
+
+!************************************
+! Matrix Initialisation
+!************************************
+open(unit=10,file='inputfile.dat',status='old')
+read(10,*)init
+close(10)
+iprint=int(init(1))
+n=int(init(2))
+choix=int(init(3))
+solchoice=int(init(4))
+tol=init(5)
+
+!call init(iprint,n,choix,solchoice)
+
+allocate(b(n),rs(n),x(n),A(n,n),work(n),sol(n))
+if (choix<3)then
+call matrice(A,n,n,2.d0,1.d0,choix)
+else
+allocate(valeurs(n))
+open(unit=30,file='matrixfile.dat',status='old')
+do i=1,n
+read(30,*)valeurs
+A(i,:)=valeurs(:)
+enddo
+close(30)
+endif
+!********************************************
+! End Matrix initialisation
+!********************************************
+
+
+!********************************************
+!   Solution initialisation
+!********************************************
+
+call solution(A,b,sol,n,n,solchoice)
+
+!********************************************
+! End solution initialisation
+!********************************************
+
+
+!********************************************
+!   Permutation Vector initialisation
+!   and computation of initial residual
+!********************************************
+allocate(p(n))
+do i=1,n
+    p(i)=i
+end do
+dnorm2=dnrm2(n,b,1)
+if (iprint==2)then
+    print*,''
+    print*,'||r0|| =',dnorm2
+    print*,''
+end if
+
+!********************************************
+!   Computation of max b(i)
+!********************************************
+indice = idamax(n,b,1)
+
+
+beta   = b(indice)
+D_inv  = 1.d0/beta
+
+!********************************************
+! First permutation
+!********************************************
+if (indice .ne.1) then
+    call dswap(1,b(indice),1,b(1),1)
+    call dswap(n,A(1,indice),1,A(1,1),1)
+    call dswap(n,A(indice,1),n,A(1,1),n)
+    call sswap(1,p(indice),1,p(1),1)
+end if
+
+
+!********************************************
+! End permutation
+!********************************************
+
+
+! call dscal(n,0.d0,x,1)
+call dscal(n,D_inv,b,1)
+rs(1)=beta
+dnorm=dabs(beta)
+allocate(c(n),s(n))
+call dscal(n,0.d0,x,1)
+! Fin initialisation de A
+!print*,'Fin initialisation des matrices'
+k=0
+
+if (iprint==2)then
+    write(*,*)'iter., pseudo-residual norm'
+end if
+!tol=1.d-7
+delta0 = dtime(da)
+eelta0 = etime(ea)
+
+open(unit=30,file='residual.dat',status='unknown')
+write(30,*)'Iteration,  Pseudo-Residual Norm' 
+!********************************************
+! Loop
+!********************************************
+do while((dnorm.gt.tol).and.(k.lt.n))
+k=k+1;
+n_k=n-k;
+k1=k+1;
+
+if (iprint==2)then
+    write(*,10)k-1,dnorm
+end if
+write(30,10)k-1,dnorm
+10 format(i5,2x,e12.4)
+
+!********************************************
+!   workl = A*b
+!********************************************
+call dcopy(n,A(1,k),1,work,1)
+call dgemv('N',n,n_k,1.d0,A(1,k1),n,b(k1),1,1.d0,work,1)
+
+!    ! call dscal(nl,0.d0,workl,1)
+!    ! do j=k1,n
+!    !    do i=1,nl
+!    !       workl(i)=workl(i)+AL(i,j)*b(j)
+!    !    end do
+!    ! end do
+
+!********************************************
+!   A(:,k)=b
+!********************************************
+call dcopy(n_k,b(k1),1,A(k1,k),1)
+
+
+!********************************************
+!   j Loop
+!********************************************
+do j=1,k
+!********************************************
+!   A(j,k)=workl(j)
+!   workl(j)=0
+!********************************************
+    A(j,k)=work(j)
+    work(j)=0.d0
+    alpha=-A(j,k)
+!********************************************
+!   workl=workl - alpha A(:,j)
+!********************************************
+    call daxpy(n-j,alpha,A(j+1,j),1,work(j+1),1)
+end do
+!********************************************
+! End j Loop
+!********************************************
+if (k.lt.n)then
+!********************************************
+!   Search indice i_0
+!********************************************
+    indice=k+idamax(n_k,work(k1),1)
+    H_k1=work(indice)
+!********************************************
+!   Permutation
+!********************************************
+    if (indice.ne.k1) then
+        call sswap(1,p(indice),1,p(k1),1)
+        call dswap(1,work(indice),1,work(k1),1)
+        call dswap(n,A(1,indice),1,A(1,k1),1)
+        call dswap(n,A(indice,1),n,A(k1,1),n)
+    end if
+!********************************************
+!   End Permutation
+!********************************************
+    if (H_k1.ne.0.d0)then
+!********************************************
+!   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)
+    end if
+else    
+    H_k1=0.d0
+    call dscal(n,0.d0,b,1)
+end if
+!********************************************
+!   Application of Givens rotations
+!********************************************
+if (k.gt.1)then
+    do i=2,k
+        btemp=A(i-1,k)
+        A(i-1,k)=c(i-1)*btemp+s(i-1)*A(i,k)
+        A(i,k)=-s(i-1)*btemp+c(i-1)*A(i,k)
+    end do
+end if
+
+if (k.le.n)then
+    btemp=A(k,k)
+    dgam=dsqrt(btemp**2+H_k1**2)
+    c(k)=dabs(btemp)/dgam
+    s(k)=(H_k1*dabs(btemp))/(btemp*dgam)
+    rs(k1)=-s(k)*rs(k)
+    rs(k)=c(k)*rs(k)
+    dnorm=dabs(rs(k1))
+    A(k,k)=c(k)*btemp+s(k)*H_k1
+end if
+!********************************************
+! End Givens rotations
+!********************************************
+end do
+!********************************************
+!   Solve H_k d_k = alpha e_1
+!********************************************
+call dtrsv('U','N','N',k,A,n,rs,1)
+call dcopy(k,rs,1,work,1)
+call dtrmv('L','N','U',k,A,n,rs,1)
+call dgemv('N',n-k,k,1.d0,A(k1,1),n,work,1,0.d0,c,1)
+call daxpy(k,1.d0,rs,1,x,1)
+call daxpy(n-k,1.d0,c,1,x(k1),1)
+
+
+
+
+
+do j=1,n
+b(p(j))=x(j)
+end do
+
+
+!!***********************************************
+!!allocate(b(n),rs(n),x(n),A(n,n),work(n),sol(n))
+!if (choix<3)then
+!call matrice(A,n,n,2.d0,1.d0,choix)
+!else
+!allocate(valeurs(n))
+!open(unit=30,file='matrixfile.dat',status='old')
+!do i=1,n
+!read(30,*)valeurs
+!A(i,:)=valeurs(:)
+!enddo
+!close(30)
+!endif
+!!********************************************
+!! End Matrix initialisation
+!!********************************************
+
+
+!!********************************************
+!!   Solution initialisation
+!!********************************************
+
+!call solution(A,work,sol,n,n,solchoice)
+
+!!*********************************************
+!!   Residual and Error Norm
+!!*********************************************
+!call dgemv('N',n,n,-1.d0,A(k1,1),n,b,1,1.d0,work,1)
+!dnorm=dnorm2(n,work,1)
+!call dcopy(n,sol,1,rs,1)
+!call daxpy(n,-1.d0,b,1,rs,1)
+!dnorm2=dnorm2(n,rs,1)
+!!***********************************************
+
+
+delta = dtime(da)
+eelta = etime(ea)
+dt2   = delta-delta0
+et2   = eelta-eelta0
+t2    = dt2+et2
+tm2   = 0.5*(t2)
+!tm2=0.d0
+if (iprint >= 2) then
+print*,'Results : iter, pseudo-residual norm, time'
+write(*,40)k,dnorm, tm2
+end if
+write(30,20)k,dnorm
+close(30)
+open(unit=100,file='tempsseq.dat',status='unknown',access='append')
+write(unit=100,FMT=*)'Iteration Number :',k,' temps :',tm2
+close(100)
+!20 format(i5,2x,e12.4,e12.4)
+
+
+if (iprint==2) then
+    write(*,*)
+    write(*,*)' Swap, CMRH Solution, Exact Solution : '
+    write(*,*)' ______________________________________________ '
+    if (n.lt.10)then
+        do i=1,n
+            write(*,20)p(i), b(i), sol(i)
+        end do
+    else
+        do i=1,10
+            write(*,20)p(i), b(i), sol(i)
+        end do
+    end if
+end if
+
+30 format(A37,4(f14.7,1x))
+20 format(1x,i5,2(1x,f14.8))
+40 format(1x,i5,2(1x,e12.4))
+
+
+
+deallocate(A,b,sol,work,c,s,rs,x)
+deallocate(p)
+
+
+end program CMRH
+!********************************************
+!    End CMRH
+!********************************************
+
+
+
+!********************************************
+!   Matrix Choice
+!********************************************
+subroutine matrice(A,nl,n,alpha,beta,choix)
+integer            ::n,i,j,choix,nl
+real(kind=8),dimension(nl,n)::A
+real*8               ::alpha,beta
+if (choix==1) then
+do i = 1, nl
+    do j = 1, i
+        A(i,j) = (2.d0*dble(j)-1.d0)/dble(n-i+j)
+    end do
+    do j = i+1, n
+        A(i,j) = (2.d0*dble(i)-1.d0)/dble(n-i+j)
+    end do
+end do
+end if
+if (choix==2)then
+    do i=1,nl
+        do j=1,i
+            A(i,j)=dble(i)
+        end do
+        do j=i+1,n
+            A(i,j)=dble(j)
+        end do
+    end do
+end if
+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
+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
+
+
+
+
+!subroutine init(iprint,n,choix,choixsol)
+!integer :: iprint,n,choix,choixsol
+!print*,'Would you print results ? : 1-No / 2-All / 3-Time and final residual norm'
+!read*,iprint
+!print*,'Enter matrix size : '
+!read*,n
+!print*,'Enter matrix choice : (more details in README.txt) '
+!print*,'1- A(i,j)=2j-1 / n-i+j for j=1,k   '
+!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- Matrix File (copy in matrixfile.dat) '
+!read*,choix
+!print*,'Enter solution choice : 1- [1,...,1]^T  / 2- [n,...,1]^T '
+!read*,choixsol
+!endsubroutine init
+
+
+
+
+
+
+
+
+
+

Fortran/PCMRH.f90

@@ -0,0 +1,576 @@
+program PCMRH
+!***************************************************************
+! Parallel CMRH program with MPI
+!***************************************************************
+use mpi
+external SSWAP,DCOPY,DSCAL,DGEMV,DSWAP
+intrinsic dabs,dsqrt
+!********************************************
+! VARIABLES
+!********************************************
+integer,parameter   ::Nprocs_max=16
+integer             ::rang,Nprocs,n,nl,code,i,j,k,choix,solchoice,iprint
+real(kind=8) ::init(5)
+integer,dimension(MPI_STATUS_SIZE) ::statut
+real(kind=8),allocatable,dimension(:,:)  ::AL
+real(kind=8)                           ::btemp,bmax,btemp2
+real(kind=8),allocatable,dimension(:)   ::b,sol,rs,temp,x,workl,work,tempn,c,s,xl,cl,temp2,xx,yy,valeurs
+integer,allocatable,dimension(:) ::p
+real*8                          ::dnorm2,D_inv,beta,dnorm,dnrm2,tol,alpha,H_k1,dgam,r
+integer                         ::indice,idamax,k1,n_k,knl,knlp,inl,inlp
+real*8 dtime, delta, delta0, dt1, dt2, dt3
+real*8 etime, eelta, eelta0, et1, et2, et3
+real*8 t1, t2, t3, tm1, tm2, tm3
+real*8 dtime1, delta1, delta01
+real*8 etime1, eelta1, eelta01,l,p1,p2
+integer nb_ligne
+real(kind=8)  ::temps,tempsd,tempsf,temps2
+integer*4,dimension(3)  ::timearray
+real                     ::rand
+real*4 da(2),ea(2),da1(2),ea1(2)
+!********************************************
+! MPI CONFIG
+!********************************************
+call MPI_INIT(code)
+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) )
+
+!************************************
+! Matrix Initialisation
+!************************************
+if (rang==0) then
+    open(unit=10,file='inputfile.dat',status='old')
+    read(10,*)init
+    close(10)
+    iprint=int(init(1))
+    n=int(init(2))
+    choix=int(init(3))
+    solchoice=int(init(4))
+    tol=init(5)
+end if
+tempsd=MPI_WTIME()
+call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,code)
+call MPI_BCAST(choix,1,MPI_INTEGER,0,MPI_COMM_WORLD,code)
+call MPI_BCAST(tol,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,code)
+
+if (mod(n,Nprocs)==0)then
+    nl=n/Nprocs
+else
+    nl=(n/Nprocs)+1
+end if
+
+allocate(b(n),rs(n),x(n),AL(nl,n),workl(nl),sol(n),temp(n),temp2(n))
+if (choix<3)then
+    call matrice(AL,nl,n,2.d0,1.d0,choix,rang)
+else
+    allocate(valeurs(n))
+    open(unit=30,file='matrixfile.dat',status='old')
+    do i=1,n
+        call position(i,nl,inl,inlp)
+        read(30,*)valeurs
+        if (rang.eq.inl)then
+            AL(inlp,:)=valeurs(:)
+        endif
+    enddo
+    close(30)
+endif
+!********************************************
+! End Matrix initialisation
+!********************************************
+
+
+!********************************************
+!   Solution initialisation
+!********************************************
+
+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)
+
+!********************************************
+! End solution initialisation
+!********************************************
+
+
+!********************************************
+!   Permutation Vector initialisation
+!   and computation of initial residual
+!********************************************
+if (rang==0) then
+    allocate(p(n))
+    do i=1,n
+        p(i)=i
+    end do
+end if
+dnorm2=dnrm2(n,b,1)
+if (rang==0)then
+    if (iprint==2)then
+        print*,''
+        print*,'la norme 2 de r0=',dnorm2
+        print*,''
+    end if
+end if
+
+!********************************************
+!   Computation of max b(i) 
+!********************************************
+indice = idamax(n,b,1)
+
+
+beta   = b(indice)
+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
+end if
+
+!********************************************
+! End permutation
+!********************************************
+
+
+! call dscal(n,0.d0,x,1)
+call dscal(n,D_inv,b,1)
+rs(1)=beta
+dnorm=dabs(beta)
+allocate(tempn(nl),work(n),c(n),s(n),xl(nl),cl(nl))
+call dscal(nl,0.d0,xl,1)
+
+call MPI_BARRIER(MPI_COMM_WORLD,code)
+k=0
+if (rang==0) then
+    if (iprint==2)then
+        write(*,*)'iter., pseudo residu '
+    end if
+    open(unit=30,file='presidual.dat',status='unknown')
+    write(30,*)'Iter., Pseudo Residual Norm'
+end if
+!tol=1.d-7
+delta0 = dtime(da)
+eelta0 = etime(ea)
+tempsf=MPI_WTIME()
+
+!********************************************
+! Loop
+!********************************************
+do while((dnorm.gt.tol).and.(k.lt.n))
+    k=k+1;
+    n_k=n-k;
+    k1=k+1;
+    if (rang==0) then
+        if (iprint==2)then
+            write(*,10)k-1,dnorm
+        end if
+        write(30,10)k-1,dnorm
+    end if
+10 format(i5,2x,e12.4)
+
+!********************************************
+!   workl = A*b
+!********************************************
+    call dgemv('N',nl,n_k+1,1.d0,AL(1,k),nl,b(k),1,0.d0,workl,1)
+
+!    ! call dscal(nl,0.d0,workl,1)
+!    ! do j=k1,n
+!    !    do i=1,nl
+!    !       workl(i)=workl(i)+AL(i,j)*b(j)
+!    !    end do
+!    ! end do
+
+!********************************************
+!   A(:,k)=b
+!********************************************
+    call position(k1,nl,knl,knlp)
+    if (rang>knl)then
+        call dcopy(nl,b(rang*nl+1),1,AL(1,k),1)
+    else
+        if (rang==knl) then
+            call dcopy(nl-knlp+1,b(k1),1,AL(knlp,k),1)
+        end if
+    end if
+
+    !********************************************
+    !   j Loop
+    !********************************************
+    do j=1,k
+        !********************************************
+        !   A(j,k)=workl(j)
+        !   workl(j)=0
+        !********************************************
+        call position(j,nl,knl,knlp)
+        if (rang==knl) then
+            AL(knlp,k)=workl(knlp)
+            workl(knlp)=0.d0
+            alpha=-AL(knlp,k)
+        end if
+        call MPI_BCAST(alpha,1,MPI_DOUBLE_PRECISION,knl,MPI_COMM_WORLD,code)
+        !********************************************
+        !   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)
+        end if
+        if (rang==knl) then
+            call daxpy(nl-knlp+1,alpha,AL(knlp,j),1,workl(knlp),1)
+        end if
+    end do
+    !********************************************
+    ! 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
+        !********************************************
+        indice=k+idamax(n_k,b(k1),1)
+        H_k1=b(indice)
+        !********************************************
+        !   Permutation
+        !********************************************
+        if (indice.ne.k1) then
+            if (rang==0) then
+                call sswap(1,p(indice),1,p(k1),1)
+            end if
+            call dswap(1,b(indice),1,b(k1),1)
+            call dswap(nl,AL(1,indice),1,AL(1,k1),1)
+            call MPI_BARRIER(MPI_COMM_WORLD,code)
+            call position(k1,nl,knl,knlp)
+            call position(indice,nl,inl,inlp)
+            if (inl==knl) then
+                if (rang==inl)then
+                    call dswap(n,AL(inlp,1),nl,AL(knlp,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==knl)then
+                    call dcopy(n,AL(knlp,1),nl,temp2,1)
+                end if
+                call MPI_BCAST(temp2,n,MPI_DOUBLE_PRECISION,knl,MPI_COMM_WORLD,code)
+                if (rang==inl)then
+                    call dcopy(n,temp2,1,AL(inlp,1),nl)
+                end if
+                if (rang==knl)then
+                    call dcopy(n,temp,1,AL(knlp,1),nl)
+                end if
+            end if
+        end if
+        !********************************************
+        !   End Permutation
+        !********************************************
+        if (H_k1.ne.0.d0)then
+            !********************************************
+            !   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)
+        end if
+    else
+        H_k1=0.d0
+        call dscal(n,0.d0,b,1)
+    end if
+    !********************************************
+    !   Application of Givens rotations
+    !********************************************
+    if (k.gt.1)then
+        do i=2,k
+            call position(i-1,nl,knl,knlp)
+            call position(i,nl,inl,inlp)
+            if (inl==knl)then
+                if (rang==inl)then
+                    btemp=AL(knlp,k)
+                    AL(knlp,k)=c(i-1)*btemp+s(i-1)*AL(inlp,k)
+                    AL(inlp,k)=-s(i-1)*btemp+c(i-1)*AL(inlp,k)
+                end if
+            else
+                if (rang==inl)then
+                    btemp=AL(inlp,k)
+                end if
+                if (rang==knl)then
+                    btemp2=AL(knlp,k)
+                end if
+                call MPI_BCAST(btemp,1,MPI_DOUBLE_PRECISION,inl,MPI_COMM_WORLD,code)
+                call MPI_BCAST(btemp2,1,MPI_DOUBLE_PRECISION,knl,MPI_COMM_WORLD,code)
+                net=net+2
+                if (rang==inl)then
+                    AL(inlp,k)=-s(i-1)*btemp2+c(i-1)*btemp
+                end if
+                if (rang==knl)then
+                    AL(knlp,k)=c(i-1)*btemp2+s(i-1)*btemp
+                end if
+            end if
+        end do
+    end if
+    if (k.le.n)then
+        call position(k,nl,knl,knlp)
+        if (rang==knl)then
+            btemp=AL(knlp,k)
+        end if
+        call MPI_BCAST(btemp,1,MPI_DOUBLE_PRECISION,knl,MPI_COMM_WORLD,code)
+        dgam=dsqrt(btemp**2+H_k1**2)
+        c(k)=dabs(btemp)/dgam
+        s(k)=(H_k1*dabs(btemp))/(btemp*dgam)
+        rs(k1)=-s(k)*rs(k)
+        rs(k)=c(k)*rs(k)
+        dnorm=dabs(rs(k1))
+        if (rang==knl)then
+            AL(knlp,k)=c(k)*btemp+s(k)*H_k1
+        end if
+    end if
+    !********************************************
+    ! End Givens rotations
+    !********************************************
+end do
+
+!********************************************
+!   Solve H_k d_k = alpha e_1
+!********************************************
+call position(k,nl,knl,knlp)
+do j=k,1,-1
+    call position(j,nl,inl,inlp)
+    if (rang==inl) then
+        rs(j)=rs(j)/AL(inlp,j)
+        btemp=rs(j)
+    end if
+    call MPI_BCAST(btemp,1,MPI_DOUBLE_PRECISION,inl,MPI_COMM_WORLD,code)
+    rs(j)=btemp
+    work(j)=btemp
+    call position(j-1,nl,inl,inlp)
+    if (rang==inl)then
+        call daxpy(inlp,-btemp,AL(1,j),1,rs(rang*nl+1),1)
+    end if
+    if (rang.lt.inl)then
+        call daxpy(nl,-btemp,AL(1,j),1,rs(rang*nl+1),1)
+    end if
+end do
+
+do j=k,1,-1
+    btemp=rs(j)
+    call position(k,nl,knl,knlp)
+    call position(j+1,nl,inl,inlp)
+    if (rang==inl)then
+        if (rang==knl)then
+            call daxpy((knlp-inlp+1),btemp,AL(inlp,j),1,rs(rang*nl+inlp),1)
+        end if
+        if (rang.lt.knl)then
+            call daxpy(nl-inlp+1,btemp,AL(inlp,j),1,rs(rang*nl+inlp),1)
+        end if
+    end if
+    if (rang.gt.inl)then
+        if (rang==knl)then
+            call daxpy(knlp,btemp,AL(1,j),1,rs(rang*nl+1),1)
+        end if
+        if (rang.lt.knl)then
+            call daxpy(nl,btemp,AL(1,j),1,rs(rang*nl+1),1)
+        end if
+    end if
+end do
+!********************************************
+!   End Loop
+!********************************************
+
+!********************************************
+!   Update x_k=x_0 + L_k d_k
+!********************************************
+call position(k1,nl,knl,knlp)
+if (rang.gt.knl)then
+    call dgemv('N',nl,k,1.d0,AL(1,1),nl,work,1,0.d0,c(rang*nl+1),1)
+end if
+if (rang==knl)then
+    call dgemv('N',nl-knlp+1,k,1.d0,AL(knlp,1),nl,work,1,0.d0,c(rang*nl+knlp),1)
+end if
+call position(k,nl,knl,knlp)
+if (rang.lt.knl)then
+    call daxpy(nl,1.d0,rs(rang*nl+1),1,xl,1)
+end if
+if (rang==knl)then
+    call daxpy(knlp,1.d0,rs(rang*nl+1),1,xl,1)
+    call daxpy(nl-knlp,1.d0,c(rang*nl+knlp+1),1,xl(knlp+1),1)
+end if
+if (rang.gt.knl)then
+    call daxpy(nl,1.d0,c(rang*nl+1),1,xl,1)
+end if
+
+!********************************************
+!   Reorder the components of x_k
+!********************************************
+call MPI_GATHER(xl,nl,MPI_DOUBLE_PRECISION,x,nl,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,code)
+call MPI_BARRIER(MPI_COMM_WORLD,code)
+if (rang==0)then
+    do j=1,n
+        b(p(j))=x(j)
+    end do
+end if
+
+
+tempsf=(MPI_WTIME()-tempsf)
+tempsd=(MPI_WTIME()-tempsd)
+!write(*,*)'rang et temps',rang,tempsf,tempsd
+call MPI_REDUCE(tempsf,temps,1,MPI_DOUBLE_PRECISION,MPI_MAX,0,MPI_COMM_WORLD,code)
+call MPI_REDUCE(tempsd,temps2,1,MPI_DOUBLE_PRECISION,MPI_MAX,0,MPI_COMM_WORLD,code)
+if (rang==0)then
+    if (iprint >= 2) then
+        print*,'Results : iter, norm residual, time'
+        write(*,40)k,dnorm,temps,temps2
+    end if
+    write(30,10),k,dnorm
+    close(30)
+    open(unit=10,file='temps.txt',status='unknown',access='append')
+    write(unit=10,FMT=*)'Proc. Number :',Nprocs,'Iter. Number:',k,'Time:',temps,temps2
+    close(10)
+end if
+
+if (rang==0)then
+    if (iprint==2) then
+        write(*,*)
+        write(*,*)' la perm. et les solutions CMRH et EXACTE sont : '
+        write(*,*)' ______________________________________________ '
+        if (n.lt.10)then
+            do i=1,n
+                write(*,20)p(i), b(i), sol(i)
+            end do
+        else
+            do i=1,10
+                write(*,20)p(i), b(i), sol(i)
+            end do
+        end if
+    end if
+end if
+
+30 format(A37,4(f14.7,1x))
+20 format(1x,i5,2(1x,f14.8))
+40 format(1x,i5,4(1x,e12.4))
+
+
+
+!  deallocate(AL,b,sol,workl,temp,tempn,work,c,s,rs,xl,x,cl)
+if (rang==0) then
+    deallocate(p)
+end if
+call MPI_FINALIZE(code)
+
+
+end program PCMRH
+!********************************************
+!    End Parallel CMRH
+!********************************************
+
+
+
+!********************************************
+!   Matrix Choice
+!********************************************
+subroutine matrice(A,nl,n,alpha,beta,choix,rang)
+integer            ::n,i,j,choix,rang,nl,k
+real(kind=8),dimension(nl,n)::A
+real*8               ::alpha,beta
+if (choix==1) then
+    do i = 1, nl
+        k=rang*nl+i
+        do j = 1, k
+            A(i,j) = (2.d0*dble(j)-1.d0)/dble(n-k+j)
+        end do
+        do j = k+1, n
+            A(i,j) = (2.d0*dble(k)-1.d0)/dble(n-k+j)
+        end do
+    end do
+end if
+if (choix==2)then
+    do i=1,nl
+        k=rang*nl+i
+        do j=1,k
+            A(i,j)=dble(k)
+        end do
+        do j=k+1,n
+            A(i,j)=dble(j)
+        end do
+    end do
+end if
+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)
+end subroutine solution
+
+
+!********************************************
+!   Where are we ?
+!********************************************
+subroutine position(k,nl,knl,knlp)
+integer :: k,nl,knl,knlp
+knl=(k-1)/nl
+knlp=k-knl*nl
+end subroutine position
+
+
+
+
+
+
+
+
+
+
+
+

Fortran/helsing.f90

@@ -0,0 +1,59 @@
+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
+
+
+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
+
+
+

Fortran/initialisation.f90

@@ -0,0 +1,43 @@
+program initialisation
+! This program creates a file inputfile.dat
+! 
+
+use randvectors
+integer :: iprint,n,choix,choixsol
+real*8 ::choixtol
+print*,'Would you print results ? : 1-No / 2-All / 3-Time and final residual norm'
+read*,iprint
+print*,'Enter matrix size : '
+read*,n
+print*,'Enter matrix choice : (more details in README.txt) '
+print*,'1- A(i,j)=2j-1 / n-i+j for j=1,k   '
+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) '
+read*,choix
+if (choix==3) then
+    call rand0(n)
+    call vecrand(n)
+    call helsing(n)
+end if
+if (choix==4) then
+    call rand0(n)
+    call vecrand(n)
+    call rbf(n)
+end if
+
+
+print*,'Enter solution choice : 1- [1,...,1]^T  / 2- [n,...,1]^T '
+print*,'3-rand vector'
+read*,choixsol
+
+print*,'Enter tol :'
+read*,choixtol
+
+open(unit=10,file='inputfile.dat',status='unknown')
+write(10,*) iprint,n,choix,choixsol,choixtol
+close(10)
+end  program initialisation

Fortran/makefile

@@ -0,0 +1,73 @@
+# OSX fortran compilers
+F=gfortran
+# Linux fortran compilers
+#F=ifort
+M=mpif90
+
+OPTS= -o
+# BLAS and LAPACK DIRECTORY for intel
+#OPTS2=-lm -L /opt/intel/mkl721/lib/64 -lmkl_lapack -lmkl_ipf -lguide -lpthread
+# BLAS and LAPACK DIRECTORY for OSX
+#OPTS2=-DYA_BLAS -DYA_LAPACK -DYA_BLASMULT -framework vecLib
+# BLAS and LAPACK DIRECTORY on baru
+OPTS2= -lm -L/usr/lib -llapack -lblas 
+
+
+FILE= randvectors.f90
+
+
+all:
+	$(F) $(OPTS) rand0 rand0.f90
+	$(F) $(OPTS) vecrand vecrand.f90
+	$(F) $(OPTS) helsing helsing.f90
+	$(F) $(OPTS) rbf rbf.f90
+	$(F) -c randvectors.f90
+	$(F) $(OPTS) initialisation initialisation.f90 $(FILE) 
+	$(F) $(OPTS) CMRH CMRH.f90 $(OPTS2)
+	$(M) $(OPTS) PCMRH PCMRH.f90 $(OPTS2)
+
+allseq: 
+	$(F) $(OPTS) rand0 rand0.f90
+	$(F) $(OPTS) vecrand vecrand.f90
+	$(F) $(OPTS) helsing helsing.f90
+	$(F) $(OPTS) rbf rbf.f90
+	$(F) -c randvectors.f90
+	$(F)  initialisation.f90 $(FILE) $(OPTS) initialisation
+	$(F) $(OPTS) CMRH CMRH.f90 $(OPTS2)
+
+allpar:
+	$(F) $(OPTS) rand0 rand0.f90
+	$(F) $(OPTS) vecrand vecrand.f90
+	$(F) $(OPTS) helsing helsing.f90
+	$(F) $(OPTS) PCMRH PCMRH.f90 $(OPTS2)	
+rand0:
+	$(F) $(OPTS) rand0 rand0.f90
+vecrand:
+	$(F) $(OPTS) vecrand vecrand.f90
+helsing:
+	$(F) $(OPTS) helsing helsing.f90
+PCMRH:
+	$(M) $(OPTS) PCMRH PCMRH.f90 $(OPTS2)
+CMRH:
+	$(F) $(OPTS) CMRH CMRH.f90 $(OPTS2)    
+clean:
+	rm *.dat
+	rm *.o
+	rm *.mod    
+	rm rand0
+	rm vecrand
+	rm helsing
+	rm rbf
+	rm initialisation
+	rm CMRH
+	rm PCMRH    
+
+testCMRH:
+	./initialisation
+	./CMRH
+
+testPCMRHbaru:
+	touch ~/.mpd.conf
+	chmod 600 ~/.mpd.conf
+	mpd &
+	mpirun -np 4 PCMRH

Fortran/rand0.f90

@@ -0,0 +1,27 @@
+program rand0
+! rand0.f90
+! This file creates a random vector of size n and stores it in rand0.dat
+
+
+  implicit none
+  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) )
+  open(unit=20,file='inputfile.dat',status='old')
+  read(20,*)init
+  n=init(1)
+  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 program rand0

Fortran/randvectors.f90

@@ -0,0 +1,187 @@
+module randvectors
+! This file contains some subroutines : 
+! rand0, vecrand
+! helsing, rbf
+
+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
+

Fortran/rbf.f90

@@ -0,0 +1,65 @@
+program rbf
+! rbf.f90
+! This file creates a Radial basis function matrix.
+! 
+!
+! We store this matrix in matrixfile.dat
+! 
+
+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
+
+open(unit=20,file='inputfile.dat',status='old')
+read(20,*)init
+n=init(1)
+close(20)
+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 program rbf
+
+
+

Fortran/vecrand.f90

@@ -0,0 +1,38 @@
+program vecrand
+! vecrand.f90 
+! This file creates two random vectors and stores them in the file vecrand.dat
+!
+! 
+ 
+  implicit none
+  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) )
+
+open(unit=20,file='inputfile.dat',status='old')
+read(20,*)init
+n=init(1)
+  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 program vecrand

Matlab/CMRH.m

@@ -0,0 +1,136 @@
+function [x,tnorm,flag]=CMRH(A,b,tol,maxiter,x0,iprint)
+%
+% [x,tnorm,flag] = (A,b,tol,maxiter,x0,iprint)
+%
+% This function solves the linear system Ax=b using the CMRH method
+%
+% Input :   
+%           A       : the n times n matrix
+%           b       : the right hand side vector
+%           tol     : the stopping criterion 
+%           maxiter : the maximum number of iterations
+%           x0      : the initial approximation
+%           iprint  : iprint = 1 to print some components of the 
+%                     approximated solution and the norm of the residual
+%                     in the Matlab Command Window.
+%                   : iprint = 0 if we do not want to print the components
+%                     and the norm of the residuals in the matlab Command 
+%                     Window. 
+%
+% Output :  
+%           x       : the solution 
+%           tnorm   : the pseudo residual norm array
+%           flag    : the stopping criterion uses
+%                   if flag = 0  : we use tol
+%                   if flag = 1  : we use maxiter 
+%
+
+n=length(A(:,1));
+
+if (nargin < 5)
+    x = x0;
+else
+    x = zeros(n,1);
+end
+if (nargin < 4)
+    maxiter = n;
+end
+if (nargin < 3)
+    tol = 1.d-7;
+end
+
+rs = zeros(n,1); s = zeros(n,1); c = zeros(n,1);
+p = 1:n; b = b-A*x;
+[~,indice] = max(abs(b(1:n)));
+beta = b(indice); D_inv = 1/beta; rs(1) = beta;
+dnorm = abs(beta); tnorm = dnorm;
+if ((indice < 1)||(indice > 1))
+    temp = p(1); p(1) = p(indice); p(indice) = temp;
+    temp = b(1); b(1) = b(indice); b(indice) = temp;
+    Ai = A(:,indice);
+    A(:,indice) = A(:,1); A(:,1) = Ai; 
+    Aip = A(indice,:); A(indice,:)=A(1,:); A(1,:) = Aip;
+end
+b = D_inv.*b;
+
+k = 0;
+while((dnorm > tol)&&(k < maxiter))
+    k = k+1;
+    n_k = n-k;
+    k1 = k+1;
+    if (iprint==1)
+       fprintf(' %3.0f %6.4e \n',k-1,dnorm);
+    end
+    work = A(:,k);
+    work = A(:,k1:n)*b(k1:n)+work;
+    A(k1:n,k) = b(k1:n);
+    for j = 1:k
+        A(j,k) = work(j);
+        work(j) = 0.d0;
+        j1 = j+1;
+        work(j1:n) = -A(j,k).*A(j1:n,j)+work(j1:n);
+    end
+    if (k < n)
+        [~,indice] = max(abs(work(k1:n)));
+        indice = k+indice;
+        H_k1 = work(indice);
+        if ((indice < k1)||(indice > k1))
+            temp = p(k1); p(k1) = p(indice); p(indice) = temp;
+            temp = work(k1); work(k1) = work(indice); work(indice) = temp;
+            Ai = A(:,indice); A(:,indice) = A(:,k1);
+            A(:,k1) = Ai; Aip = A(indice,:); A(indice,:) = A(k1,:);
+            A(k1,:) = Aip;
+        end
+        if ((H_k1 < 0)||(H_k1 > 0))
+            D_inv = 1/H_k1; b = work; b = D_inv.*b;
+        end
+    else
+        H_k1 = 0;
+        b = zeros(n,1);
+    end
+    
+    % solving the least squares problem    
+    if (k > 1)
+        for i = 2:k
+            i_1 = i-1;
+            dtemp = A(i_1,k);
+            A(i_1,k) = c(i_1)*dtemp+s(i_1)*A(i,k);
+            A(i,k) = -s(i_1)*dtemp+c(i_1)*A(i,k);
+        end
+    end
+    if (k <= n)
+        dgam = sqrt(A(k,k)^2+H_k1^2);
+        c(k) = abs(A(k,k))/dgam;
+        s(k) = (H_k1*abs(A(k,k)))/(A(k,k)*dgam);
+        rs(k1) = -s(k)*rs(k);
+        rs(k) = c(k)*rs(k);
+        A(k,k) = c(k)*A(k,k)+s(k)*H_k1;
+        dnorm = abs(rs(k1));
+    end
+    tnorm = [tnorm dnorm];
+end
+if (k==maxiter)
+    flag=1;
+else
+    flag=0;
+end
+if (iprint==1)
+    fprintf(' %3.0f %6.4e \n',k,dnorm);
+end
+
+% computing d the solution of the least square problem
+for j = k:-1:1
+    rs(j) = rs(j)/A(j,j);
+    rs(1:j-1) = rs(1:j-1)-rs(j).*A(1:j-1,j);
+end
+work(1:k) = rs(1:k);
+for j = k:-1:1
+    rs(j+1:k) = rs(j+1:k)+rs(j).*A(j+1:k,j);
+end
+c = A(k1:n,1:k)*work(1:k);
+x(1:k) = rs(1:k)+x(1:k);
+x(k1:n) = c(1:n_k)+x(k1:n);
+
+% Reorder the components of x to get the solution
+xx(p(1:n)) = x(1:n);
+x = xx';

Matlab/helsing.m

@@ -0,0 +1,39 @@
+function A = helsing(n)
+%
+% A = helsing(n)
+% 
+% This file creates a matrix A whose entries are 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. 
+
+%
+% Input : 
+%           n       : size of the matrix A
+% Output : 
+%           A       : square matrix of size n as described above 
+%
+% 07 February 13
+%
+xi=rand(n,1)*2-1;
+yi=rand(n,1)*2-1;
+[xj,yj]=meshgrid(xi,yi);
+xj=xj'-xj;
+yj=yj'-yj;
+coef=-1;
+for i=1:max(size(xj))
+    d0=0;
+    for j=1:n
+        temp=sqrt((xi(i)-xi(j))^2+(yi(j)-yi(i))^2);
+        if (temp>0)
+           A(i,j)=coef*log(temp);
+           if (d0==0)
+               d0=j;
+               A(i,i)=-log(0.5*d0*rand());
+           end
+        end
+    end
+end

Matlab/main.m

@@ -0,0 +1,79 @@
+% main.m
+%
+% This script allows to perform some numerical tests by giving the results
+% obtained when solving a linear system A*x=b by using the CMRH algorithm
+%
+% You have to choose :
+%    - the size of the system, 
+%    - the coefficient matrix, 
+%    - the exact solution (so that we can compute the error), 
+%    - the tolerance 
+%
+
+
+fprintf('Hello world \n')
+fprintf('Enter the matrix choice :\n')
+fprintf('1- A(i,j)= (alpha*j-beta) / (n-i+j)  if j <= i    \n')
+fprintf('           (alpha*i-beta) / (n-i+j)  if j > i   \n') 
+fprintf('2- A(i,j)= i if j <= i    \n')
+fprintf('           j if j > i  \n')
+fprintf('3- A(i,j)= abs(i-j)+1/(i-j) if j < i  and j > i  \n')
+fprintf('           0  if i = j   \n')
+fprintf('4- A(i,j)= i+j    \n')
+fprintf('5- helsing matrix \n')
+fprintf('6- Matrix Market matrix \n')
+matrixchoice=input('Enter your choice : ');
+if (matrixchoice <6)
+    n=input('Enter the size of the matrix : ');
+    if (matrixchoice <5)
+         A = matrice(n,matrixchoice);
+    end
+    if (matrixchoice==5)
+        A = helsing(n);
+    end
+else
+    filename = input('Enter the matrix file name :','s');
+    [A,rows,cols,entries,rep,field,symm] = mmread(filename);
+    n = rows;
+end
+
+fprintf('Enter the solution choice :\n')
+fprintf('1-sol=random vector \n')
+fprintf('2-sol=[1,1,...,1]^T \n')
+fprintf('3-sol=[1,2,...,n]^T \n')
+solchoice = input('Enter your choice : ');
+[b,sol] = solution(A,solchoice);
+
+tol = input('Enter the tolerance : ');
+iprint = input('Print results : 1-Yes, 2-No : ');
+
+
+V = A; bc = b;
+[x,pseudores,flag] = CMRH(A,b,tol,n,zeros(n,1),iprint);
+
+fid = fopen('pseudores.dat','w');
+for i=1:length(pseudores)
+    fprintf(fid,'%i \n',pseudores(i));
+end
+fclose(fid); 
+
+b = x;
+if (iprint==1)
+    fprintf('\n approx. and exact solutions :\n');
+    if (n < 10)
+       for i = 1:n
+         fprintf(' %6.4e %6.4e\n',b(i),sol(i));
+       end
+    else
+       for i = 1:10
+         fprintf(' %6.4e %6.4e\n',b(i),sol(i));
+       end
+    end
+    % we compute the norm of the residu
+    bc = -V*b+bc;
+    re_cmrh = norm(bc,2);
+    % we compute the norm of the error
+    sol = -b+sol;
+    er_cmrh = norm(sol,2);
+    fprintf('residual norm: %6.4e  error norm: %6.4e \n',re_cmrh,er_cmrh);
+end

Matlab/matrice.m

@@ -0,0 +1,63 @@
+function A=matrice(n,i_mat)
+%
+% A = matrice(n,i_mat)
+
+% This function allows to choose the coefficient matrix of
+% the linear system A*x=b.  
+%
+% Input :   
+%           n       : the size of the coefficient matrix of 
+%                     the linear system A*x=b
+%           i_mat   : integer (i_rhs = 1 or 2 or 3 or 4)
+% if i_mat = 1 then A(i,j) = (alpha*j-beta) / (n-i+j)  if j <= i 
+%              and  A(i,j) = (alpha*i-beta) / (n-i+j)  if j > i   
+% if i_mat = 2 then A(i,j) = i if j <= i   
+%              and  A(i,j) = j if j > i 
+% if i_mat = 3 then A(i,j) = abs(i-j)+1/(i-j) if (j < i  and j > i)
+%              and  A(i,j) = 0  if i=j   \n')
+% if i_mat = 4 then A(i,j)= i+j
+%
+% Output :  
+%           A       : the coefficient matrix of the linear system A*x=b
+%
+%
+if (i_mat==1)
+    a=input('Enter alpha : ');
+    b=input('Enter beta : ');
+    for i=1:n
+        for j=1:i
+            A(i,j)=(a*j-b)/(n-i+j);
+        end
+        for j=i+1:n
+            A(i,j)=(a*i-b)/(n-i+j);
+        end
+    end
+end
+if (i_mat==2)
+   for i=1:n
+       for j=1:i
+           A(i,j)=i;
+       end
+       for j=i+1:n
+           A(i,j)=j;
+       end
+   end
+end
+if (i_mat==3)
+   for i=1:n
+       for j=1:i-1
+           A(i,j)=abs(i-j)+1/(i-j);
+       end
+       for j=i+1:n
+           A(i,j)=abs(i-j)+1/(i-j);
+       end
+       A(i,i)=0;
+   end
+end
+if (i_mat==4)
+   for i=1:n
+       for j=1:n
+           A(i,j)=i+j;
+       end
+   end
+end

Matlab/mmread.m

@@ -0,0 +1,221 @@
+function  [A,rows,cols,entries,rep,field,symm] = mmread(filename)
+%
+% function  [A] = mmread(filename)
+%
+% function  [A,rows,cols,entries,rep,field,symm] = mmread(filename)
+%
+%      Reads the contents of the Matrix Market file 'filename'
+%      into the matrix 'A'.  'A' will be either sparse or full,
+%      depending on the Matrix Market format indicated by
+%      'coordinate' (coordinate sparse storage), or
+%      'array' (dense array storage).  The data will be duplicated
+%      as appropriate if symmetry is indicated in the header.
+%
+%      Optionally, size information about the matrix can be 
+%      obtained by using the return values rows, cols, and
+%      entries, where entries is the number of nonzero entries
+%      in the final matrix. Type information can also be retrieved
+%      using the optional return values rep (representation), field,
+%      and symm (symmetry).
+%
+
+mmfile = fopen(filename,'r');
+if ( mmfile == -1 )
+ disp(filename);
+ error('File not found');
+end;
+
+header = fgets(mmfile);
+if (header == -1 )
+  error('Empty file.')
+end
+
+% NOTE: If using a version of Matlab for which strtok is not
+%       defined, substitute 'gettok' for 'strtok' in the 
+%       following lines, and download gettok.m from the
+%       Matrix Market site.    
+[head0,header]   = strtok(header);  % see note above
+[head1,header]   = strtok(header);
+[rep,header]     = strtok(header);
+[field,header]   = strtok(header);
+[symm,header]    = strtok(header);
+head1 = lower(head1);
+rep   = lower(rep);
+field = lower(field);
+symm  = lower(symm);
+if ( length(symm) == 0 )
+   disp(['Not enough words in header line of file ',filename]) 
+   disp('Recognized format: ')
+   disp('%%MatrixMarket matrix representation field symmetry')
+   error('Check header line.')
+end
+if ( ~ strcmp(head0,'%%MatrixMarket') )
+   error('Not a valid MatrixMarket header.')
+end
+if (  ~ strcmp(head1,'matrix') )
+   disp(['This seems to be a MatrixMarket ',head1,' file.']);
+   disp('This function only knows how to read MatrixMarket matrix files.');
+   disp('  ');
+   error('  ');
+end
+
+% Read through comments, ignoring them
+
+commentline = fgets(mmfile);
+while length(commentline) > 0 & commentline(1) == '%',
+  commentline = fgets(mmfile);
+end
+
+% Read size information, then branch according to
+% sparse or dense format
+
+if ( strcmp(rep,'coordinate')) %  read matrix given in sparse 
+                              %  coordinate matrix format
+
+  [sizeinfo,count] = sscanf(commentline,'%d%d%d');
+  while ( count == 0 )
+     commentline =  fgets(mmfile);
+     if (commentline == -1 )
+       error('End-of-file reached before size information was found.')
+     end
+     [sizeinfo,count] = sscanf(commentline,'%d%d%d');
+     if ( count > 0 & count ~= 3 )
+       error('Invalid size specification line.')
+     end
+  end
+  rows = sizeinfo(1);
+  cols = sizeinfo(2);
+  entries = sizeinfo(3);
+  
+  if  ( strcmp(field,'real') )               % real valued entries:
+  
+    [T,count] = fscanf(mmfile,'%f',3);
+    T = [T; fscanf(mmfile,'%f')];
+    if ( size(T) ~= 3*entries )
+       message = ...
+       str2mat('Data file does not contain expected amount of data.',...
+               'Check that number of data lines matches nonzero count.');
+       disp(message);
+       error('Invalid data.');
+    end
+    T = reshape(T,3,entries)';
+    A = sparse(T(:,1), T(:,2), T(:,3), rows , cols);
+  
+  elseif   ( strcmp(field,'complex'))            % complex valued entries:
+  
+    T = fscanf(mmfile,'%f',4);
+    T = [T; fscanf(mmfile,'%f')];
+    if ( size(T) ~= 4*entries )
+       message = ...
+       str2mat('Data file does not contain expected amount of data.',...
+               'Check that number of data lines matches nonzero count.');
+       disp(message);
+       error('Invalid data.');
+    end
+    T = reshape(T,4,entries)';
+    A = sparse(T(:,1), T(:,2), T(:,3) + T(:,4)*sqrt(-1), rows , cols);
+  
+  elseif  ( strcmp(field,'pattern'))    % pattern matrix (no values given):
+  
+    T = fscanf(mmfile,'%f',2);
+    T = [T; fscanf(mmfile,'%f')];
+    if ( size(T) ~= 2*entries )
+       message = ...
+       str2mat('Data file does not contain expected amount of data.',...
+               'Check that number of data lines matches nonzero count.');
+       disp(message);
+       error('Invalid data.');
+    end
+    T = reshape(T,2,entries)';
+    A = sparse(T(:,1), T(:,2), ones(entries,1) , rows , cols);
+
+  end
+
+elseif ( strcmp(rep,'array') ) %  read matrix given in dense 
+                               %  array (column major) format
+
+  [sizeinfo,count] = sscanf(commentline,'%d%d');
+  while ( count == 0 )
+     commentline =  fgets(mmfile);
+     if (commentline == -1 )
+       error('End-of-file reached before size information was found.')
+     end
+     [sizeinfo,count] = sscanf(commentline,'%d%d');
+     if ( count > 0 & count ~= 2 )
+       error('Invalid size specification line.')
+     end
+  end
+  rows = sizeinfo(1);
+  cols = sizeinfo(2);
+  entries = rows*cols;
+  if  ( strcmp(field,'real') )               % real valued entries:
+    A = fscanf(mmfile,'%f',1);
+    A = [A; fscanf(mmfile,'%f')];
+    if ( strcmp(symm,'symmetric') | strcmp(symm,'hermitian') | strcmp(symm,'skew-symmetric') ) 
+      for j=1:cols-1,
+        currenti = j*rows;
+        A = [A(1:currenti); zeros(j,1);A(currenti+1:length(A))];
+      end
+    elseif ( ~ strcmp(symm,'general') )
+      disp('Unrecognized symmetry')
+      disp(symm)
+      disp('Recognized choices:')
+      disp('   symmetric')
+      disp('   hermitian')
+      disp('   skew-symmetric')
+      disp('   general')
+      error('Check symmetry specification in header.');
+    end
+      A = reshape(A,rows,cols);
+  elseif  ( strcmp(field,'complex'))         % complx valued entries:
+    tmpr = fscanf(mmfile,'%f',1);
+    tmpi = fscanf(mmfile,'%f',1);
+    A  = tmpr+tmpi*i;
+    for j=1:entries-1
+      tmpr = fscanf(mmfile,'%f',1);
+      tmpi = fscanf(mmfile,'%f',1);
+      A  = [A; tmpr + tmpi*i];
+    end
+    if ( strcmp(symm,'symmetric') | strcmp(symm,'hermitian') | strcmp(symm,'skew-symmetric') ) 
+      for j=1:cols-1,
+        currenti = j*rows;
+        A = [A(1:currenti); zeros(j,1);A(currenti+1:length(A))];
+      end
+    elseif ( ~ strcmp(symm,'general') )
+      disp('Unrecognized symmetry')
+      disp(symm)
+      disp('Recognized choices:')
+      disp('   symmetric')
+      disp('   hermitian')
+      disp('   skew-symmetric')
+      disp('   general')
+      error('Check symmetry specification in header.');
+    end
+    A = reshape(A,rows,cols);
+  elseif  ( strcmp(field,'pattern'))    % pattern (makes no sense for dense)
+   disp('Matrix type:',field)
+   error('Pattern matrix type invalid for array storage format.');
+  else                                 % Unknown matrix type
+   disp('Matrix type:',field)
+   error('Invalid matrix type specification. Check header against MM documentation.');
+  end
+end
+
+%
+% If symmetric, skew-symmetric or Hermitian, duplicate lower
+% triangular part and modify entries as appropriate:
+%
+
+if ( strcmp(symm,'symmetric') )
+   A = A + A.' - diag(diag(A));
+   entries = nnz(A);
+elseif ( strcmp(symm,'hermitian') )
+   A = A + A' - diag(diag(A));
+   entries = nnz(A);
+elseif ( strcmp(symm,'skew-symmetric') )
+   A = A - A';
+   entries = nnz(A);
+end
+
+fclose(mmfile);
+% Done.

Matlab/pseudores.dat

@@ -0,0 +1,109 @@
+9.646458e+02 
+8.627479e+01 
+1.812252e+01 
+1.716083e+01 
+1.463363e+01 
+4.528704e+00 
+2.557136e+00 
+2.553890e+00 
+2.275236e+00 
+2.218364e+00 
+1.498835e+00 
+5.573258e-01 
+5.459666e-01 
+3.582639e-01 
+2.367821e-01 
+1.475634e-01 
+1.451316e-01 
+1.236858e-01 
+1.005615e-01 
+7.235333e-02 
+7.233631e-02 
+7.221167e-02 
+7.193137e-02 
+7.190673e-02 
+7.167662e-02 
+7.071333e-02 
+6.984837e-02 
+6.915964e-02 
+5.174372e-02 
+5.170931e-02 
+5.074308e-02 
+4.616475e-02 
+3.282722e-02 
+3.171966e-02 
+2.287480e-02 
+2.093466e-02 
+1.168468e-02 
+1.167402e-02 
+1.148494e-02 
+9.699070e-03 
+6.892333e-03 
+5.824841e-03 
+5.329200e-03 
+5.004072e-03 
+4.866720e-03 
+4.769666e-03 
+4.206442e-03 
+4.202169e-03 
+4.200567e-03 
+4.188529e-03 
+4.187999e-03 
+4.187547e-03 
+4.182589e-03 
+4.158763e-03 
+4.087599e-03 
+3.991127e-03 
+3.655560e-03 
+3.410989e-03 
+3.048420e-03 
+3.047078e-03 
+2.956584e-03 
+2.466037e-03 
+2.465243e-03 
+2.459474e-03 
+2.413068e-03 
+2.360192e-03 
+1.785296e-03 
+1.782533e-03 
+1.700346e-03 
+1.522626e-03 
+9.479767e-04 
+9.143709e-04 
+5.507576e-04 
+3.056827e-04 
+3.041580e-04 
+2.992387e-04 
+2.990824e-04 
+2.817873e-04 
+2.587030e-04 
+1.871993e-04 
+1.868418e-04 
+1.851624e-04 
+1.811406e-04 
+1.806046e-04 
+1.778918e-04 
+1.676887e-04 
+1.647738e-04 
+1.646303e-04 
+1.496040e-04 
+1.400422e-04 
+1.297343e-04 
+1.296961e-04 
+1.010297e-04 
+7.573647e-05 
+7.570585e-05 
+7.546339e-05 
+7.489631e-05 
+7.337235e-05 
+6.505731e-05 
+6.288240e-05 
+5.011360e-05 
+4.761092e-05 
+2.725433e-05 
+2.653621e-05 
+1.640640e-05 
+1.640219e-05 
+1.636241e-05 
+1.508365e-05 
+8.775428e-06 

Matlab/solution.m

@@ -0,0 +1,33 @@
+function [b,sol]=solution(A,i_rhs)
+%
+% [b,sol] = solution(A,i_rhs)
+%
+% This function allows to generate a right hand side b for
+% the linear system A*x=b so as that the exact solution is known.
+%
+% Input :   
+%           A       : the coefficient matrix of the linear system A*x=b
+%           i_rhs   : integer (i_rhs = 1 or 2 or 3)
+% if i_rhs = 1 then the entries of the exact solution of the linear system 
+% are generated randomly in [0, 1]
+% if i_rhs = 2 then the entries of the exact solution of the linear 
+% system equals 1.
+% if i_rhs = 3 then the k-th entry of the exact solution of the linear 
+% system equals k.
+%
+% Output :  
+%           b       : the right hand side of the linear system A*x=b
+%           sol     : the exact solution of the linear system A*x=b
+%
+n = length(A)
+    
+if (i_rhs==1)
+   sol=rand(n,1);
+end
+if (i_rhs==2)
+   sol=ones(n,1);
+end
+if (i_rhs==3)
+   sol=[1:n]';
+end
+b=A*sol;

Readme.txt

@@ -0,0 +1,124 @@
+  ***************************************************************************
+  * All the software  contained in this library  is protected by copyright. *
+  * Permission  to use, copy, modify, and  distribute this software for any *
+  * purpose without fee is hereby granted, provided that this entire notice *
+  * is included  in all copies  of any software which is or includes a copy *
+  * or modification  of this software  and in all copies  of the supporting *
+  * documentation for such software.                                        *
+  ***************************************************************************
+  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED *
+  * WARRANTY. IN NO EVENT, NEITHER  THE AUTHORS, NOR THE PUBLISHER, NOR ANY *
+  * MEMBER  OF THE EDITORIAL BOARD OF  THE JOURNAL  "NUMERICAL ALGORITHMS", *
+  * NOR ITS EDITOR-IN-CHIEF, BE  LIABLE FOR ANY ERROR  IN THE SOFTWARE, ANY *
+  * MISUSE  OF IT  OR ANY DAMAGE ARISING OUT OF ITS USE. THE ENTIRE RISK OF *
+  * USING THE SOFTWARE LIES WITH THE PARTY DOING SO.                        *
+  ***************************************************************************
+  * ANY USE  OF THE SOFTWARE  CONSTITUTES  ACCEPTANCE  OF THE TERMS  OF THE *
+  * ABOVE STATEMENT.                                                        *
+  ***************************************************************************
+
+     AUTHORS:
+        S. Duminil, Ph. Marion, H. Sadok
+	 University of Littoral Calais, France
+	 E-mail: duminil@lmpa.univ-littoral.fr, sadok@lmpa.univ-littoral.fr
+		
+	 Mohammed Heyouni
+        University Mohammed Premier Oujda, Maroc
+        E-mail: mohammed.heyouni@gmail.com
+
+      REFERENCE:
+	
+	- CMRH: A new method for solving non symmetric linear systems based on the  	  Hessenberg reduction algorithms, Numer. Algorithms, 20, (1999), pp 303-321
+        - A new implementation of the CMRH method for solving dense linear systems, 
+	  J. Comput. Appl. Math., 213 (2008), pp 387-399
+	- A parallel implementation of the CMRH method for dense linear systems, 
+	  Numer. Algorithms, 63 (2013), 127-142
+
+	SOFWARE REVISION DATE:
+	  v1.0, January 2013
+
+	SOFTWARE LANGUAGE:
+	  Fortran 90 and MATLAB 8.0 (R2012b)
+
+=============================================================================
+SOFTWARE
+=============================================================================
+
+This package provides a set of Fortran and Matlab functions to solve linear systems using CMRH method. 
+
+=============================================================================
+PACKAGE
+=============================================================================
+
+The main directory contains the following files
+
+* README.txt	 :	This file
+
+* Fortran 90 subroutines : 
+
+-makefile	      	: to install and config all fortran programs
+-inputfile.dat      	: In this file, we choose the size of system, the choice of 			  matrix, the choice of solution and the print version
+-initialisation.f90 	: This file creates inputfile.dat 
+-CMRH.f90      	: CMRH program
+-PCMRH.f90		: Parallel CMRH program
+-helsing.f90		: Matrix test stored in matrixfile.dat
+-rbf.f90		: Matrix test stored in matrixfile.dat
+-rand0.f90		: Create a rand vector stored in rand0.dat
+-vecrand.f90		: Create two rand vectors stored in vecrand.dat
+-randvectors.f90  	: Contains some subroutines  
+
+* Matlab programs :
+-main.m  		: main program
+-CMRH.m			: CMRH function
+-matrice.m		: To create some matrix example
+-solution.m		: To create some solution vector example
+-mmread.m		: from matrix market. To store matrix market file in A
+-helsing.m		: Matrix test
+-gemat11.mtx		: A matrix market example
+
+
+=============================================================================
+HOW TO INSTALL
+=============================================================================
+When the archive file containing the package is uncompressed, it creates a directory named CMRHprograms. 
+The directory fortran contains the fortran codes. 
+The directory matlab contains the matlab codes.
+
+* Fortran 90 programs :
+- enter to fortran directory
+- Open makefile
+- Config fortran, mpi, blas and lapack compilers 
+
+* Matlab programs :
+- enter to matlab directory
+- Run main.m
+
+=============================================================================
+HOW TO COMPILE
+=============================================================================
+- make all (to config all programs) or make allseq (to config only sequential programs) 
+- ./initialisation  (to config inputfile.dat)
+- ./CMRH  (to run CMRH program)
+- mprun -n (nb of proc) ./PCMRH (to run PCMRH program)
+
+
+=============================================================================
+NUMERICAL TESTS
+=============================================================================
+FORTRAN PROGRAMS :
+- Run ./initialisation to config the matrix :  
+	-helsing : Matrix arising from Helsing paper (choose 3)
+	-rbf   	 : Matrix arising from radial basis functions (choose 4)
+	-Own matrix : save it in matrixfile.dat  (choose 5)
+- Run ./CMRH or mprun -n (nb of proc) ./PCMRH 
+
+MATLAB PROGRAMS :
+- Run main.m
+- Choose the matrix example.
+- For matrix market test, choose 6 and enter the file name .mtx 
+
+
+=============================================================================
+DELETE FORTRAN FILES
+=============================================================================
+- make clean

git.readme.txt

@@ -0,0 +1,6 @@
+
+
+
+
+1. Récuperez l'archive du serveur
+git clone git@piccolo.univ-littoral.fr:cmrh_NAlgorithm