program CMRH
!***************************************************************
! CMRH program
!
! Created by H. Sadok and M. Heyouni -1999
! Updated by S. Duminil dec. 2012
!***************************************************************

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