! Ejemplo del Metodo de Diferencias Finitas para resolver la ecuacion diferencial parcial:
!    Uxx=-Pi*Pi*cos(Pix)
!        xi<=U<=xf
!        U(xi)=vi y U(xf)=vf
! con solucion analitica: 
!     cos(pix)
!
! Para conocer mas del metodo ver:
!    http://132.248.182.159/acl/MDF
!
! Autor: Antonio Carrillo Ledesma
!    http://mmc.geofisica.unam.mx/acl


	program mdf1D
		implicit none
		integer i, N, nn 
		real*16, allocatable :: A(:,:), b(:), x(:)
		real*16 xi, xf, vi, vf, h, R, P, Q, y

		xi = 0.0                    ! Inicio del diminio
		xf = 1.0                    ! Fin del dominio
		vi = 1.0                    ! Valor en la frontera xi
		vf = -1.0                   ! Valor en la frontera xf
		nn = 11                     ! Particion
		N = nn - 2                  ! Nodos interiores
		h = (xf - xi) / (nn-1)   ! Incremento en la malla

    		allocate (A(N,N), b(N), x(N))
                ! Stencil
		R = 1. / (h * h)
		P = -2. / (h * h)
		Q = 1. / (h * h)

		! Primer renglon de la matriz A y vector b
		A(1,1)=P
		A(2,1)=Q
		call ladoDerecho(xi,y)
		b(1)=y-vi*R
		! Renglones intermedios de la matriz A y vector b
		do i=2,N-1 
			A(i-1,i)=R
			A(i,i)=P
			A(i+1,i)=Q
			call ladoDerecho(xi+h*(i),y)
			b(i)= y
		end do
		! Renglon final de la matriz A y vector b
		A(N-1,N)=R
		A(N,N)=P
		call ladoDerecho(xi+h*N,y)
		b(N)= y -vf*Q

		call gaussSeidel(A, x, b, N, 1000)
                print*, "A: ", A
                print*, "b: ", b
		print*, "x: ", x

                call jacobi(A, x, b, N, 1000)
                print*, "A: ", A
                print*, "b: ", b
                print*, "x: ", x

	end program


	subroutine ladoDerecho(x,y)	
		real*16, intent(in) :: x
		real*16, intent(inout) :: y
		real*16 pi

		pi = 3.1415926535897932384626433832;
		y = -pi * pi * cos(pi * x);
	end subroutine


	subroutine gaussSeidel(a, x, b, nx, iter)
		implicit none
		integer, intent(in) :: nx, iter
		real*16, intent(in) :: a(nx,nx), b(nx)
		real*16, intent(inout) :: x(nx)
		integer i, j, m
		real*16 sum
		
		do m = 1, iter
			do i = 1, nx
				sum = 0.0
				do j = 1, nx
					if (i .NE. j) then
						sum = sum + a(i,j) * x(j)
					end if
				end do
				x(i) = (1.0 / a(i,i)) * (b(i) - sum)
			end do
		end do
	end subroutine


        subroutine jacobi(a, x, b, nx, iter)
                implicit none
                integer, intent(in) :: nx, iter
                real*16, intent(in) :: a(nx,nx), b(nx)
                real*16, intent(inout) :: x(nx)
                real*16, allocatable :: xt(:)
                integer i, j, m
                real*16 sum

                allocate (xt(nx))
                do m = 1, iter
                        do i = 1, nx
                                sum = 0.0
                                do j = 1, nx
                                        if (i .NE. j) then
                                                sum = sum + a(i,j)*x(j)
                                        end if
                                end do
                                xt(i) = (1.0 / a(i,i)) * (b(i) - sum)
                        end do
                        do i = 1, nx
                                x(i) = xt(i)
                        end do
                end do
        end subroutine