!
!===============================================================
! Imprimir matriz simetrica en perfil por columnas
!===============================================================
!
!     NOTA: Se omiten todas las tildes y caracteres especiales en los
!           textos para ajustarse al estandar ASCII y favorecer la
!           compatibilidad del archivo con todos los equipos y sistemas.
!
! AUTORES: J. Paris, I. Couceiro, F. Navarrina, I. Colominas
!                                                                 
!        
! FECHA ULTIMA REVISION: 20260201
!
! Copyright: Este programa es propiedad de su autor.
!            No puede utilizarse con fines comerciales sin conocimiento
!            y autorizacion de su autor.
!            Puede utilizarse con fines educativos sin ninguna restriccion.
!
! Definicion de variables de la subrutina:    [I_8 = integer*8, R_8 = real*8]
!                                             [I = input, O = output]
!
!     n  [I_8] [I]   =  dimension de la matriz (numero de puntos de discretizacion)
!     lp [R_8] [I]   =  Vector de dimension n en el que se almacena la posicion en
!                       el vector a de cada una de las diagonales de la matriz
!     a  [R_8] [I]   =  Vector de dimension lp(n) en el que se almacena
!                       la matriz simetrica en perfil columnas
!     b  [R_8] [I]   =  Vector de dimension n en el que se almacena 
!                       inicialmente el termino independiente y luego la solucion
!
      subroutine imprime_matriz_cp(n,a,lp)
      implicit integer*8(i-n)
      implicit real*8(a-h,o-z)
      dimension lp(n),a(lp(n))
      character*10 espacios
      
      espacios='          '
      
      log_mat=10
      
      if (n.gt.50)then
        write(6,*)'El problema es demasiado grande',
     &            ' para escribir su matriz en un archivo'
        write(6,*)'   Para continuar el calculo sin escribir la',
     &            ' matriz pulsa Enter'
        write(6,*)'   Para detener el calculo pulsa CTRL + C'
        read(5,*)
        return
      endif
      
      open(unit=log_mat,file='matriz_cp.txt',status='unknown')
      
  100 format(e10.2,$)
  
      write(log_mat,100)a(1)
      write(log_mat,*)
      write(log_mat,*)
      do i=2,n
        nancho=i-(lp(i)-lp(i-1))
        do j=1,nancho ! Rellenamos con espacios la parte exterior al perfil
          write(log_mat,'(a,$)')espacios
        enddo
        do j=lp(i-1)+1,lp(i)
          write(log_mat,100)a(j)  ! Imprimimos los coeficiente del perfil en su ubicacion real
        enddo
        write(log_mat,*)
        write(log_mat,*)
      enddo
      
      close(log_mat)

      return
      end
      
!
!===============================================================
! Factorizacion LDL^T con matriz almacenada en perfil por columnas
!===============================================================
!
!     NOTA: Se omiten todas las tildes y caracteres especiales en los
!           textos para ajustarse al estandar ASCII y favorecer la
!           compatibilidad del archivo con todos los equipos y sistemas.
!
! AUTORES: J. Paris, I. Couceiro, F. Navarrina, I. Colominas
!                                                                 
!        
! FECHA ULTIMA REVISION: 20260205
!
! Copyright: Este programa es propiedad de su autor.
!            No puede utilizarse con fines comerciales sin conocimiento
!            y autorizacion de su autor.
!            Puede utilizarse con fines educativos sin ninguna restriccion.
!
! Definicion de variables de la subrutina:    [I_8 = integer*8, R_8 = real*8]
!                                             [I = input, O = output]
!
!     n  [I_8] [I]   =  dimension de la matriz (numero de puntos de discretizacion)
!     lp [R_8] [I]   =  Vector de dimension n en el que se almacena la posicion en
!                       el vector a de cada una de las diagonales de la matriz
!     a  [R_8] [I/O] =  Vector de dimension lp(n) en el que se almacena
!                       la matriz simetrica en perfil columnas
!
      subroutine ldlt_cp(n,a,lp)
      implicit integer*8(i-n)
      implicit real*8(a-h,o-z)
      dimension lp(n),a(lp(n))
      
      do k=1,n-1
        lbk=lp(k+1)-lp(k)-1  ! Ancho del perfil de la fila k+1 exluyendo la diagonal
        lpk0=lp(k+1)-(k+1)
        kaux=k+1-lbk
        do i=kaux+1,k
          lbi= lp(i)-lp(i-1)-1 ! Ancho del perfil de la fila i exluyendo la diagonal,
                               ! ya que si el vector lp esta bien definido el bucle
                               ! empieza como minimo en i=2. Si i=1 no sirve esta ecuacion   
          lpij0=lp(i)-i
          s=0.d+00
          do j=max(i-lbi,kaux),i-1
            lpij=lpij0+j  ! Ubicacion de a(i,j)
            lpkj=lpk0+j  ! Ubicacion de a(k+1,j)
            s=s+a(lpij)*a(lpkj)
          enddo
          lpki=lpk0+i  ! Ubicacion de a(k+1,i)
          a(lpki)=a(lpki)-s
        enddo
        
        do i=k+1-lbk,k
          lpki=lpk0+i  ! Ubicacion de a(k+1,i)
          lpii=lp(i)  ! Ubicacion de a(i,i)
          a(lpki)=a(lpki)/a(lpii)
        enddo
        
        s=0.d+00
        do j=k+1-lbk,k
          lpjj=lp(j)  ! Ubicacion de a(j,j)
          lpkj=lpk0+j  ! Ubicacion de a(k+1,j)
          s=s+a(lpjj)*a(lpkj)**2
        enddo
        lpkk=lp(k+1)  ! Ubicacion de a(k+1,k+1)
        a(lpkk)=a(lpkk)-s
      enddo
      
      return
      end
!
!===============================================================
! Sistema triangular inferior con matriz en perfil por columnas
!===============================================================
!
!     NOTA: Se omiten todas las tildes y caracteres especiales en los
!           textos para ajustarse al estandar ASCII y favorecer la
!           compatibilidad del archivo con todos los equipos y sistemas.
!
! AUTORES: J. Paris, I. Couceiro, F. Navarrina, I. Colominas
!                                                                 
!        
! FECHA ULTIMA REVISION: 20260205
!
! Copyright: Este programa es propiedad de su autor.
!            No puede utilizarse con fines comerciales sin conocimiento
!            y autorizacion de su autor.
!            Puede utilizarse con fines educativos sin ninguna restriccion.
!
! Definicion de variables de la subrutina:    [I_8 = integer*8, R_8 = real*8]
!                                             [I = input, O = output]
!
!     n  [I_8] [I]   =  dimension de la matriz (numero de puntos de discretizacion)
!     lp [R_8] [I]   =  Vector de dimension n en el que se almacena la posicion en
!                       el vector a de cada una de las diagonales de la matriz
!     a  [R_8] [I/O] =  Vector de dimension lp(n) en el que se almacena
!                       la matriz simetrica en perfil columnas ya factorizada.
!     b  [R_8] [I/O] =  Vector de dimension n que almacena inicialmente
!                       el vector de terminos independientes y devuelve
!                       el vector solucion del sistema.
!
      
      subroutine lx_b_cp(n,a,lp,b)
      implicit integer*8(i-n)
      implicit real*8(a-h,o-z)
      dimension lp(n),a(lp(n)),b(n)


      do i=2,n
        s=0.d+00
        lbi= lp(i)-lp(i-1)-1 ! Ancho del perfil de la fila i exluyendo la diagonal
        lpi0=lp(i)-i
        do j=i-lbi,i-1
          lpij=lpi0+j  ! Ubicacion de a(i,j)
          s=s+a(lpij)*b(j)
        enddo
        b(i)=b(i)-s
      enddo
      
      return
      end

!
!===============================================================
! Sistema diagonal con matriz en perfil por columnas
!===============================================================
!
!     NOTA: Se omiten todas las tildes y caracteres especiales en los
!           textos para ajustarse al estandar ASCII y favorecer la
!           compatibilidad del archivo con todos los equipos y sistemas.
!
! AUTORES: J. Paris, I. Couceiro, F. Navarrina, I. Colominas
!                                                                 
!        
! FECHA ULTIMA REVISION: 20260205
!
! Copyright: Este programa es propiedad de su autor.
!            No puede utilizarse con fines comerciales sin conocimiento
!            y autorizacion de su autor.
!            Puede utilizarse con fines educativos sin ninguna restriccion.
!
! Definicion de variables de la subrutina:    [I_8 = integer*8, R_8 = real*8]
!                                             [I = input, O = output]
!
!     n  [I_8] [I]   =  dimension de la matriz (numero de puntos de discretizacion)
!     lp [R_8] [I]   =  Vector de dimension n en el que se almacena la posicion en
!                       el vector a de cada una de las diagonales de la matriz
!     a  [R_8] [I/O] =  Vector de dimension lp(n) en el que se almacena
!                       la matriz simetrica en perfil columnas ya factorizada.
!     b  [R_8] [I/O] =  Vector de dimension n que almacena inicialmente
!                       el vector de terminos independientes y devuelve
!                       el vector solucion del sistema.
!
      subroutine dx_b_cp(n,a,lp,b)
      implicit integer*8(i-n)
      implicit real*8(a-h,o-z)
      dimension lp(n),a(lp(n)),b(n)

      do i=1,n
        lpii=lp(i)  ! Ubicacion de a(i,i)
        b(i)=b(i)/a(lpii)
      enddo
      
      return
      end

!
!===============================================================
! Sistema triangular superior con matriz en perfil por columnas
!===============================================================
!
!     NOTA: Se omiten todas las tildes y caracteres especiales en los
!           textos para ajustarse al estandar ASCII y favorecer la
!           compatibilidad del archivo con todos los equipos y sistemas.
!
! AUTORES: J. Paris, I. Couceiro, F. Navarrina, I. Colominas
!                                                                 
!        
! FECHA ULTIMA REVISION: 20260205
!
! Copyright: Este programa es propiedad de su autor.
!            No puede utilizarse con fines comerciales sin conocimiento
!            y autorizacion de su autor.
!            Puede utilizarse con fines educativos sin ninguna restriccion.
!
! Definicion de variables de la subrutina:    [I_8 = integer*8, R_8 = real*8]
!                                             [I = input, O = output]
!
!     n  [I_8] [I]   =  dimension de la matriz (numero de puntos de discretizacion)
!     lp [R_8] [I]   =  Vector de dimension n en el que se almacena la posicion en
!                       el vector a de cada una de las diagonales de la matriz
!     a  [R_8] [I/O] =  Vector de dimension lp(n) en el que se almacena
!                       la matriz simetrica en perfil columnas ya factorizada.
!     b  [R_8] [I/O] =  Vector de dimension n que almacena inicialmente
!                       el vector de terminos independientes y devuelve
!                       el vector solucion del sistema.
!
      subroutine ltx_b_cp(n,a,lp,b)
      implicit integer*8(i-n)
      implicit real*8(a-h,o-z)
      dimension lp(n),a(lp(n)),b(n)

      do i=n,2,-1
        lbi= lp(i)-lp(i-1)-1 ! Ancho del perfil de la fila i exluyendo la diagonal
        lpi0=lp(i)-i
        do j=i-lbi,i-1
          lpij=lpi0+j  ! Ubicacion de a(i,j)
          b(j)=b(j)-a(lpij)*b(i) 
        enddo
      enddo

      return
      end 

