C ---------------------------------------------------------
      subroutine creartindpL(b,tbl)
c
c Proposito: Crear termino independiente para el problema de 
c            obtencion del campo electrico en una region cuadrada
c            con un potencial inicial:
c              100V
c          --------
c          |      | 100V
c    100V  |      | 
c          --------
c            0V
c
c input: tbl: Orden de la malla
c output: b: termino independiente
c
      implicit none
      integer tbl,i
      real*8 b(tbl*tbl)
      do i=1,tbl*tbl
         b(i)=0.d0
      enddo
      do i=tbl,tbl*tbl,tbl
         b(i)=100.d0
      enddo
      do i=1,tbl*tbl,tbl
         b(i)=b(i)+100.d0
      enddo
      do i=tbl*tbl-tbl+1,tbl*tbl
         b(i)=b(i)+100.d0
      enddo
      return
      end
C ---------------------------------------------------------
      subroutine crearmatrizL(A,m,lda,tbl)
c
c proposito: crear la matriz de Laplace
c
c input: m: numero de filas y columnas de A
c        lda: dimension principal del array A
c        tbl: Orden de la malla
c output: A: matriz de Laplace
c
      implicit none
      integer m,lda,tbl,i,j,k
      real*8 A(lda,m)

C Bloque 1
      A(1,1)=4
      A(2,1)=-1
      A(1+tbl,1)=-1
      do i=2,tbl-1
         A(i-1,i)=-1
         A(i,i)=4
         A(i+1,i)=-1
         A(i+tbl,i)=-1
      enddo
      A(tbl-1,tbl)=-1
      A(tbl,tbl)=4
      A(2*tbl,tbl)=-1
C Bloque intermedios   
      do k=2,tbl-1
         A((k-1)*tbl+1-tbl,(k-1)*tbl+1)=-1
         A((k-1)*tbl+1,(k-1)*tbl+1)=4
         A((k-1)*tbl+2,(k-1)*tbl+1)=-1
         A((k-1)*tbl+1+tbl,(k-1)*tbl+1)=-1
         do i=2,tbl-1
            A((k-1)*tbl+i-tbl,(k-1)*tbl+i)=-1
            A((k-1)*tbl+i-1,(k-1)*tbl+i)=-1
            A((k-1)*tbl+i,(k-1)*tbl+i)=4
            A((k-1)*tbl+i+1,(k-1)*tbl+i)=-1
            A((k-1)*tbl+i+tbl,(k-1)*tbl+i)=-1
         enddo
         A((k-1)*tbl,(k-1)*tbl+tbl)=-1
         A((k-1)*tbl+tbl-1,(k-1)*tbl+tbl)=-1
         A((k-1)*tbl+tbl,(k-1)*tbl+tbl)=4
         A((k-1)*tbl+2*tbl,(k-1)*tbl+tbl)=-1
      enddo
C Bloque final 
      A(tbl*(tbl-1)+1-tbl,tbl*(tbl-1)+1)=-1
      A(tbl*(tbl-1)+1,tbl*(tbl-1)+1)=4
      A(tbl*(tbl-1)+2,tbl*(tbl-1)+1)=-1
      do i=2,tbl-1
         A(tbl*(tbl-1)+i-tbl,tbl*(tbl-1)+i)=-1
         A(tbl*(tbl-1)+i-1,tbl*(tbl-1)+i)=-1
         A(tbl*(tbl-1)+i,tbl*(tbl-1)+i)=4
         A(tbl*(tbl-1)+i+1,tbl*(tbl-1)+i)=-1
      enddo
      A(tbl*(tbl-1),tbl*(tbl-1)+tbl)=-1
      A(tbl*(tbl-1)+tbl-1,tbl*(tbl-1)+tbl)=-1
      A(tbl*(tbl-1)+tbl,tbl*(tbl-1)+tbl)=4
      return
      end
C ---------------------------------------------------------
      function norma(x,y,n)
c 
c proposito: Calcular la norma_1 de la diferencia de dos vectores
c
c input: x: vector 1
c        y: vector 2
c        n: tamanio de los vectores
c output: norma: norma_1(x-y)
c
      implicit none
      integer n,i
      real*8 x(n),y(n),norma
      norma=0.d0
      do i=1,n
         norma=norma+abs(x(i)-y(i))
      enddo
      return 
      end
C ---------------------------------------------------------
      subroutine solveLU(A,n,b,x,lda)
c 
c proposito: Calcular la solucion del sistema Ax=b, donde
c            se conoce la descomposicion LU de A y esta 
c            almacenada en la propia A.
c
c input: A: array que contiene la descomposicion LU
c        n: numero de filas y columna del sistema a resolver
c        b: termino independiente
c        lda: dimension principal del array A
c output x: solucion del sistema
c
      implicit none
      integer n,maxm,i,j,lda
      parameter (maxm=1001)
      real*8 A(lda,n),x(n),b(n),y(maxm)
     
C     solucion de Ly=b (barrido de columnas)
      do j=1,n
         y(j)=b(j)
         do i=j+1,n
            b(i)=b(i)-y(j)*A(i,j)      
         enddo
      enddo
C     solucion de Ux=y (barrido de columnas)
      do j=n,1,-1
         x(j)=y(j)/A(j,j)
         do i=1,j-1
            y(i)=y(i)-x(j)*A(i,j)
         enddo
      enddo
      return
      end
C ---------------------------------------------------------
      subroutine LU(A,n,lda)
c 
c proposito: Calcular la descomposicion LU de una matriz A
c
c input: A: matriz a descomponer
c        n: numero de filas y columnas de A
c        lda: dimension principal del array A
c output: A: la descomposicion LU se almacena en la propia A
c
      implicit none
      integer n,lda,k,s,i,j
      real*8 A(lda,n)

      do k=1,n-1
         do s=k+1,n
            A(s,k)=A(s,k)/A(k,k)
         enddo
         do j=k+1,n
            do i=k+1,n
               A(i,j)=A(i,j)-A(i,k)*A(k,j)
            enddo
         enddo
      enddo
      return
      end
C ---------------------------------------------------------
      subroutine newtindp(A,b,x,m,n,s,lda,bnew)
c
c proposito: calcular el nuevo termino independiente en cada iteracion
c            de un Jacobi por bloques
c
c input: A: array correspondiente a m filas y n columnas consecutivas
c        n: numero de filas
c        m: numero de columnas
c        s: entero utilizado para conocer el bloque diagonal
c           el que va de la columna s-n+1 a la s
c        lda: dimension principal del array A
c        x: vector iterado actual de dimension m
c        b: termino independiente del sistema
c output: bnew: nuevo termino independiente de diemension n
c
      implicit none
      integer m,n,s,lda,i,j
      real*8 A(lda,m),b(n),x(m),bnew(n)
      
      do i=1,n
         bnew(i)=0.d0
      enddo
      do j=1,s-n
         do i=1,n
            bnew(i)=bnew(i)+A(i,j)*x(j)
         enddo
      enddo
      do j=s+1,m
         do i=1,n
            bnew(i)=bnew(i)+A(i,j)*x(j)
         enddo
      enddo
      do i=1,n
         bnew(i)=b(i)-bnew(i)
      enddo
      return
      end