program Sedimentation
        !---------------------------
        ! Déclaration de variables
        !---------------------------
        !Pas de variables implictes
        implicit none

        !Deux seuls vecteurs necessaires à passer en paramètres
        !de la fonction reduction (la solution est mise dans b)
        real(kind=8), dimension(:),allocatable::b,u,pprec
        !paramètres des blocs
        real(8)::alpha,beta,rho,condition
        !paramètres du domaines
        integer, parameter::l=8,n=2**(l+1)-1,m=500
        !indices de boucle
        integer::i,j,iteration,maxIteration
        !pas du domaine et de temps
        real(8) :: hx, hy, hyCarre, dt

        ! 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=10.0d0/dble(m+1)
        hy=10.0d0/dble(n+1)

        !on initialise le temps cpu
        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

        !------------
        !       Réduction
        !------------

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

        print *,'m= ',m,'n= ',n

        ! Calcul de b
        b=0
        dt=1000.0d0
        rho=1.0d0
        hyCarre=hy*hy
        pprec=0
        iteration=0
        maxIteration=100
        condition=1e-3

!calcul de beta et alpha
        alpha = hyCarre/(hx*hx)
        beta  = 2+2*alpha+rho*hyCarre/dt

        !appel de la reduction cyclique
        call cpu_time(debut)

        do while(condition .ge. 1e-5 .and. iteration < maxIteration)
                do i=1,m
                        do j=1,n
                                b(i+m*(j-1))=hyCarre*(second_membre(xg(i),yg(j))+pprec(i+m*(j-1))*rho/dt)
                        end do
                end do

                !print *,sqrt(dot_product(b,b))
                call reduction(b,n*m,m,-alpha,beta)

                condition = sqrt(dot_product(b-pprec,b-pprec)/dot_product(b,b))
                pprec = b
                iteration = iteration + 1
        end do

        print *,iteration

        call cpu_time(fin)
        !print *, 'Fin du calcul'

        !calcul de la solution de l'erreur et de la durée de calcul
        call solution(xg,yg,u,m,n)

        print *, "Norme de la solution b= ", sqrt(dot_product(b,b))
        print *, "Norme de la solution u= ", sqrt(dot_product(u,u))
		 print *, "Norme euclidienne= ", sqrt(dot_product(b-u,b-u))
        print *, "Norme H1", sqrt( dot_product(b,b)+dot_product(b-pprec,b-pprec)/(dt*dt))
        !print *, (fin - debut), ' s'

        !on libère la mémoire
        deallocate(b)
        deallocate(u)

contains
!------------------------------------------------------------
! second_membre : calcul la valeur du terme de droite
!                 pour l'EDP en 2D
! Paramètres :
!               x,y les coordonnées du point où calculer cette fonction
! Retour:
!               la valeur de cette fonction
!-------------------------------------------------------------
function second_membre(x,y)
        implicit none
        real(8) :: x,y,second_membre
        second_membre=-2*x*(x-10)-2*y*(y-10)
end function

end program Sedimentation

!------------------------------------------------------------
! solution : calcul la valeur de la solution au prolème
!            d'EDP en 2D
! Paramètres :
!   xg,yg le  maillage du domaine
!   u     le  vecteur  de retour
!   m,n   les dimensions du domaines
!-------------------------------------------------------------
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), intent(out) :: u

        integer :: i,j,m,n

        !on calcule u en chaque noeud du maillage
 u=0
        do i=1,m
                do j=1,n
                        u(i+m*(j-1)) = xg(i)*yg(j)*(xg(i)-10)*(yg(j)-10)
                enddo
        enddo

end subroutine

!-------------------------------------------------------------
! tridiag :  résolution d'un système
! Paramètres :
!   diagValue,tridiagValue les paramètres du bloc
!   b le second membre et le vecteur de retour
!               m la dimension de b
!-------------------------------------------------------------
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)
        end do

        ! phase de substitution directe
        do i=2,m
                b(i) = b(i) - e(i-1)*b(i-1)
        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)
        end do

end subroutine

!-------------------------------------------------------------
! factoAk : factorisation du calcul de la puissance kieme de a
!               dans la resolution des sous systèmes (stabilisation)
! Paramètres :
!               k l'indice de la puissance de a
!   b le second membre et le vecteur de retour
!               m la dimension de b
!   moinsAlpha,beta les paramètres du bloc
!-------------------------------------------------------------
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

        !-------------------------------------------
        ! Résolution des sous systèmes tridiagonaux
        !-------------------------------------------
        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


!-------------------------------------------------------------
! reduction : procédure principale de la réduction cyclique
! Paramètres :
!   b le second membre et le vecteur de retour
!               taille,m les dimensions de b sachant que taille=m*n
!   moinsAlpha,beta les paramètres des matrices blocs
!-------------------------------------------------------------
subroutine reduction(b,taille,m,moinsAlpha,beta)
        !--------------------------
        ! Déclaration de variables
        !--------------------------
        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, lPow
        !variables locales
        integer::jmax
        real(8)::pidiv

        !--------------------------------
        ! Initialisation de la réduction
        !--------------------------------

        !allocation des vecteurs temporaires
        allocate(p(taille))
        allocate(q(taille))
        allocate(ptemp(m))

        !initialisation des variables locales
        j = 1
        !jmax = 2**k
        !pidiv=2*atan(x=1.0d0)/jmax
        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
!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

        do i=1,2**l-1
                i1=(i+i-1)*m+1
                i2 = i1 + m - 1
                q(i1:i2)= b(i1-m:i2-m) + b(i1+m:i2+m) + 2*p(i1:i2)
        end do

        !-----------------------------------------
        ! 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
                        i2 = i1 + m - 1
                        ie=m*r/2

                        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
!--- on a plus besoin de b a partir d'ici

        !------------------------------------
        ! Calcul de la solution (mise dans b)
        !------------------------------------

        b=0
        lPow = 2**l
        ptemp = q((lPow-1)*m+1:(lPow)*m)
        call factoAk(l,ptemp,m,moinsAlpha,beta)
        b((lPow-1)*m+1:(lPow)*m)=p((lPow-1)*m+1:(lPow)*m)+ptemp


        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

        p=0
        q=0
        ptemp=0
        deallocate(p)
        deallocate(q)
        deallocate(ptemp)

end subroutine