! commentaire
program Cyclo
	! Déclaration de variables
	implicit none

	real(kind=8), dimension(:),allocatable::b,u
	real(8)::alpha,beta
	integer, parameter::l=13,n=2**(l+1)-1,m=2000
	!integer, parameter::l=3,n=2**(l+1)-1,m=15
	integer::i,j
	real(8) :: hx, hy, hyCarre
	
	! Maillage
	real(kind=8), dimension(0:m+1) :: xg
	real(kind=8), dimension(0:n+1) :: yg

	! Cpu time
	real(kind=8) :: debut, fin, dummy

	! Initialisation du maillage et calcul du pas
	xg=0.0d0
	yg=0.0d0
	hx=1/dble(m+1)
	hy=1/dble(n+1)

	call cpu_time(dummy)
	
	! création du maillage
	do i=0,m+1
		xg(i)=dble(i)*hx
	enddo
	do i=0,n+1
		yg(i)=dble(i)*hy
	enddo

	! Allocation des vecteurs
	allocate(b(n*m))
	allocate(u(n*m))

	print *,'calcul de b'
	
	! Calcul de b
	b=0
	hyCarre=hy*hy
	do i=1,m
		do j=1,n
			b(i+m*(j-1))=(hyCarre)*get_rhs(xg(i),yg(j))
		end do
	end do
	
	!calcul de beta et alpha
	alpha=hyCarre/(hx*hx)
	beta=2+2*alpha

	print *,'Reduction'
	!appel de la reduction cyclique
	call cpu_time(debut)
	call reduction(b,n*m,m,-alpha,beta)
	call cpu_time(fin)
	
	call solution(xg,yg,u,m,n)
	print *, sqrt(dot_product(b-u,b-u))
	print *, (fin - debut), ' s'
	
	deallocate(b)
	deallocate(u)

contains
function get_rhs(x,y)
	implicit none
	real(8) :: x,y,get_rhs
	get_rhs=-2*x*(x-1)-2*y*(y-1)
end function

end program Cyclo

subroutine solution(xg,yg,u,m,n)
	implicit none
	real(8), dimension(0:m+1) :: xg
	real(8), dimension(0:n+1) :: yg
	real(8), dimension(m*n) :: u
	
	integer :: i,j,m,n
	
	do i=1,m
	do j=1,n
		u(i+m*(j-1))=xg(i)*yg(j)*(xg(i)-1)*(yg(j)-1)
	enddo
	enddo

end subroutine


subroutine factoAk(k,b,m,moinsAlpha,beta)
	implicit none

	!argument de  la fonction
	real(8),intent(in)::moinsAlpha,beta
	integer,intent(in)::k,m
	real(8),dimension(m),intent(inout)::b

	!variables locales
	integer::j,jmax
	real(8)::pidiv,lambda
	real(8),dimension(m)::z

	!calcul d'une constante d'optimisation de temps de calcul
	pidiv=2*atan(x=1.0d0)/2**k

	do j=1,2**k
		lambda = 2*cos((2*j-1)*pidiv)
		call tridiag(moinsAlpha,beta - lambda,m,b(1:m))
		
	end do


end subroutine


subroutine reduction(b,taille,m,moinsAlpha,beta)
	implicit none
	
	!parametres
	real(8),intent(in)::moinsAlpha,beta
	integer,intent(in)::taille,m
	real(8),dimension(taille),intent(inout)::b
	
	!variables locales
	integer::n,l,k,i,j,inc,oldinc,imax,im,r,i1,i2,ie, ind1,ind2
	real(8),dimension(:),allocatable::p,q
	real(8),dimension(:),allocatable::ptemp
	real(8)::unsura
	!variables locales
	integer::jmax
	real(8)::pidiv
	      
	print*,'initialisation de la réduction'
	allocate(p(taille))
	allocate(q(taille))
	allocate(ptemp(m))
	
	j = 1
	jmax = 2**k
	pidiv=2*atan(x=1.0d0)/jmax	
	
	!initialisation des variables locales
	n=taille/m
	l=log(x=n+1.0)/log(x=2.0)-1.0
	imax=2**(l+1)
	q=0
	p=b
	ptemp = 0	
	
	print*,'initialisation de p et q---------'
	!Premiere Boucle : initialisation de p et q
	do i=1,n
		call tridiag(moinsAlpha,beta,m,p((i-1)*m+1:i*m))
	end do
	
	print*,'p apres initialisation---------'
	do i=1,n
		!print *,p((i-1)*m+1:i*m)
	enddo	
	
	do i=1,2**l-1
		i1=(2*i-1)*m+1
		i2=2*i*m
		q(i1:i2)= b(i1-m:i2-m) + b(i1+m:i2+m) + 2*p(i1:i2)
	end do

	print*,'q apres initialisation ---------'
	do i=1,n
		!print *,q((i-1)*m+1:i*m)
	enddo	

	!--- Deuxieme boucle : calcul de p et de q

	do k=2,l
		r=2**k
		
		do i=1,2**(l-k+1)-1
			i1=(r*i-1)*m+1
			i2=r*i*m
			ie=m*r/2
			
			!print *,i1,',',i2,',',ie
			
			ptemp(1:m) = q(i1:i2)+p(i1-ie:i2-ie)+p(i1+ie:i2+ie)

			call factoAk(k-1,ptemp,m,moinsAlpha,beta)	
				
			p(i1:i2) = p(i1:i2) + ptemp
			q(i1:i2) = 2*p(i1:i2)+q(i1-ie:i2-ie)+q(i1+ie:i2+ie)	
			
		end do
	end do

	print*,'---------p apres reduction'
	do i=1,n
		!print *,p((i-1)*m+1:i*m)
	enddo	

	print*,'---------q apres reduction'
	do i=1,n
		!print *,q((i-1)*m+1:i*m)
	enddo	
	!on a plus besoin de b a partir d'ici

	b=0
	ptemp = q((2**l-1)*m+1:(2**l)*m)
	call factoAk(l,ptemp,m,moinsAlpha,beta)
	b((2**l-1)*m+1:(2**l)*m)=p((2**l-1)*m+1:(2**l)*m)+ptemp
	
	print*,'---------y'
!	print *, "y=", (2**l-1)*m+1, (2**l)*m

	print*,'---------x'
!	print *, "x=", b((2**l-1)*m+1:2**l*m)
	
	do k=l-1,0,-1
		r=2**k
		im=2**(l-k+1)-1
		do i=1,im,2
			i1=(r*i-1)*m+1
			i2=r*i*m
			ie=m*r
			if(i == 1) then
				ptemp = b(i1+ie:i2+ie)+q(i1:i2)
			else
				if (i==im) then
					ptemp = b(i1-ie:i2-ie)+q(i1:i2)
				else
					ptemp = b(i1+ie:i2+ie)+b(i1-ie:i2-ie)+q(i1:i2)
				endif
			endif
		
			call factoAk(k,ptemp,m,moinsAlpha,beta)
			b(i1:i2)=p(i1:i2)+ptemp
		end do
	end do

	!print *, "x=", b((2**l-1)*m+1:2**l*m)

end subroutine


subroutine tridiag(diagValue,tridiagValue,m,b)
	implicit none
	real(8),intent(in)::diagValue, tridiagValue
	integer,intent(in)::m
	real(8),dimension(m),intent(inout)::b
	! d, diag
	! e, sous/sur-diag
	real(8)::r
	integer::i
	real(8),dimension(m)::e, d

	! factorisation
	d = tridiagValue
	e = diagValue
	do i=2,m
		r = e(i-1)
		e(i-1) = dble(r / d(i-1))
		d(i) = d(i) - r*e(i-1)
		!print *,'\t1\t{',d
	end do

	! phase de substitution directe
	do i=2,m
		b(i) = b(i) - e(i-1)*b(i-1)
		!print *,'\t2\t{',b
	end do
	b(m) = b(m) / d(m)
	
	! substitution inverse
	do i=m-1,1,-1
		b(i) = b(i)/d(i)-e(i)*b(i+1)
		!print *,'\t3\t{',b
	end do
	
end subroutine