c==============================================================================
c=====      SUBRUTINAS AUXILIARES DE SISTEMAS DE ECUACIONES LINEALES       ====
c==============================================================================
c
c     Estas subrutinas permiten factorizar una matriz cuadrada simetrica
c     y en banda como producto de 3 (A=L*D*Lt). La matriz del sistema A
c     estara almacenada a partir de una matriz en banda de dimensiones
c     (n,l+1) que a su vez se almacena por filas en un vector v de 
c     dimension n(l+1), siendo n el numero de ecuaciones
c     y l el semiancho de banda.
c     Como resultado, sobre el mismo vector de datos, se almacenan los
c     valores de las matrices L y D adoptando el criterio de que la
c     matriz L tiene coeficientes unitarios en la diagonal.
c     Se incluyen tambien subrutinas para la resolucion de sistemas
c     triangulares superiores, diagonales y triangulares superiores
c     con las matrices L y D resultantes de la factorizacion (que se
c     almacenan sobre la propia matriz inicial).
c
c AUTORES: J. Paris, I. Couceiro, F. Navarrina, I. Colominas
c                                                                 
c        
c FECHA ULTIMA REVISION: 20240319
c
c Copyright: Este programa es propiedad de sus autores.
c            No puede utilizarse con fines comerciales sin conocimiento
c            y autorizacion de su autor.
c            Puede utilizarse con fines educativos sin ninguna restriccion.
c
c===============================================================================
c Definicion de variables de la subrutina:    [I_8 = integer*8, R_8 = real*8]
c                                             [I = input, O = output]
c
c     v [R_8] [I/O] =  Vector de dimensiones n(l+1) 
c                      correspondiente al almacenamiento por filas en
c                      un vector de una matriz almacenada en banda de
c                      dimension (n,l+1) a partir de una matriz
c                      original de dimensiones n x n y semiancho de
c                      banda l (diagonal excluida).
c                      Inicialmente contiene la matriz del sistema. 
c                      Tras las operaciones de factorizacion contiene
c                      la informacion de las matrices L y D resultado
c                      de la factorizacion de Cholesky
c     n [I_8] [I]   =  Dimension de la matriz del problema
c     l [I_8] [I]   =  Semiancho de banda de la matriz (diag. excluida)
c
c     NOTA: Se omiten todas las tildes y caracteres especiales en los
c           textos para ajustarse al estandar ASCII y favorecer la
c           compatibilidad del archivo con todos los equipos y sistemas.
c===============================================================================
c===============================================================================
c     Subrutina para la factorizacion de matriz simetrica almacenada en banda
c===============================================================================
      subroutine fac_simetrica(v,n,l)
      implicit integer*8(i-n)
      implicit real*8(a-h,o-z)
      dimension v(n*(l+1))
      
      do k=1,n-1
        lk1=l*(k+1)
        do i=max(1,k+2-l),k
          il=i*l
          s=0.0d+00
          do j=max(1,k+1-l),i-1
            s = s+v(il+j)*v(lk1+j)
          enddo
          v(lk1+i) = v(lk1+i)-s

        enddo
        do i=max(1,k+1-l),k
          v(lk1+i) = v(lk1+i)/v(i*(l+1))
        enddo
        s=0.0d+00
        do j=max(1,k+1-l),k
          s = s+v(lk1+j)*v(j*(l+1))*v(lk1+j)
        enddo
        v(lk1+(k+1)) = v(lk1+(k+1))-s
      enddo 

      end
c===============================================================================
c     subrutina para la resolucion de sistema L·x=b
c===============================================================================
      subroutine sol_l(v,n,l,x)
      implicit integer*8(i-n)
      implicit real*8(a-h,o-z)
      dimension v(n*(l+1)),x(n)
      
      do i=2,n
        s=0.d+00
        il=i*l
        do j=max(i-l,1),i-1
          s=s+v(il+j)*x(j)
        enddo
        x(i)=x(i)-s
      enddo

      end

c===============================================================================
c     subrutina para la resolucion de sistema D·x=b
c===============================================================================
      subroutine sol_d(v,n,l,x)
      implicit integer*8(i-n)
      implicit real*8(a-h,o-z)
      dimension v(n*(l+1)),x(n)
      
      do i=1,n
         x(i)=x(i)/v(i*(l+1))
      enddo

      end
c===============================================================================
c     subrutina para la resolucion de sistema L^T·x=b
c===============================================================================
      subroutine sol_lt(v,n,l,x)
      implicit integer*8(i-n)
      implicit real*8(a-h,o-z)
      dimension v(n*(l+1)),x(n)
      
      do i=n,2,-1
        il=i*l
         do j=max(i-l,1),i-1
           x(j)=x(j)-v(il+j)*x(i)
         enddo
      enddo
 
      end
c===============================================================================
c===============================================================================
