c==============================================================================
c============    SOLUCION DE SISTEMAS DE ECUACIONES LINEALES    ===============
c==============================================================================

c______________________________________SUBRUTINA Solver_GMRES_________________
c
c AUTORES: J. Paris, F. Navarrina, H. Gomez, X. Nogueira, I. Colominas, M. Casteleiro
c                                                                _
c        Dep. de Metodos Matematicos y de Rep., ETSICCP de A Coruna, UDC.
c
c FECHA ULTIMA REVISION: 20110510
c
c Copyright: Esta subrutina es propiedad de sus autores.
c            No puede utilizarse con fines comerciales sin conocimiento
c            y autorizacion de sus autores.
c            Puede utilizarse con fines educativos sin ninguna restriccion.
c            Los comentarios que se anaden en este programa utilizan
c            los caracteres del codigo ASCII estandar. Por lo tanto, no
c            se utilizan ni la letra \~n ni las tildes.
c
c------------------------------------------------------------------------------
c
c OBJETIVO:
c
c    Solucion general del sistema de ecuaciones :
c
c                                         -   -
c                                     A . x = b
c                                     ~            
c    siendo:
c           A = matriz de coeficientes  (matriz de rigidez)         [dato]
c           x = solucion                (vector de desplazamientos) [incognita]
c           b = terminos independientes (vector de fuerzas nodales) [dato]
c------------------------------------------------------------------------------
c 
c DESCRIPCION GENERAL:
c
c       La matriz A es una matriz cualquiera no singular.
c
c       La matriz se almacena como matriz en perfil por columnas.
c
c       El sistema resultante A'. x = b' se resuelve mediante el metodo 
c       GMRES con Precondicionamiento Diagonal. 
c
c Nota: En este tipo de metodos, no es posible aprovechar los calculos
c       para resolver varios sistemas de ecuaciones simultaneamente.
c
c------------------------------------------------------------------------------
c
c VARIABLES:
c
c     ntotv = [Entrada] (INT*8) : Numero de ecuaciones del sistema
c     ntotr1= [Entrada] (INT*8) : Numero de terminos no nulos de la matriz del sistema (triangular inferior).
c     ntotr2= [Entrada] (INT*8) : Numero de terminos no nulos de la matriz del sistema (triangular superior).
c     nvfix = [Entrada] (INT*8) : Numero de incognitas con valores prescritos.
c     niter = [Entrada] (INT*8) : Numero maximo de iteraciones. Dimension del subespacio
c                                 de Krylov utilizado para el calculo. (niter <= ntotv)
c     nrestart = [Ent.] (INT*8) : Numero maximo de veces que se reinicia el metodo GMRES.
c     nactu = [Entrada] (INT*8) : Separacion entre iteraciones en las que se imprime
c                                 por pantalla el residuo. (nactu>0)
c     epslr = [Entrada] (REAL*8): Error absoluto admisible en r' (norma del maximo)
c                                 siendo r' el residuo           (r' = b' - A' x)
c     rltvr = [Entrada] (REAL*8): Error relativo admisible en r' (norma del maximo)
c                                 en relacion con b'
c     lp1   = [Entrada] (INT*8) : Vector de direccionamiento de los coeficientes de la parte
c                                 inferior de la matriz almacenada en perfil por columnas.
c     lp2   = [Entrada] (INT*8) : Vector de direccionamiento de los coeficientes de la parte
c                                 superior de la matriz almacenada en perfil por columnas.
c     rgmat_l=[Entrada] (REAL*8): Matriz de coeficientes del sistema de ecuaciones
c                                 almacenada en perfil por columnas (parte inferior)
c     rgmat_u=[Entrada] (REAL*8): Matriz de coeficientes del sistema de ecuaciones
c                                 almacenada en perfil por columnas (parte superior)
c     gload = [Entrada] (REAL*8): Vector de terminos independientes.
c     xdisp = [Ent/Sal] (REAL*8): Aproximacion inicial del vector solucion. Vector solucion
c     resdl = [Salida]  (REAL*8): Residuo r' = D^(-1) (b - A x)
c     ngfix = [Entrada] (INT*8) : lista de gdl coaccionados.
c     presc = [Entrada] (REAL*8): lista de valores prescritos en gdl coaccionados.
c     treac = [Salida ] (REAL*8): valores de las reacciones.
c     wrksp = [Ent/Sal] (REAL*8): Vector de trabajo del GMRES de dimension ((ntotv+3)*(niter+3))
c     residuo=[Salida]  (REAL*8): Norma 2 del vector residuo resdl.
c     kiter = [Salida]  (INT*8) : Numero de iteraciones realizadas en cada restart.
c     irestart = [Sal]  (INT*8) : Numero de restarts completados del metodo GMRES cuando
c                                 alcanzamos convergencia.
c     ierr  = [Salida]  (INT*8) : Variable que indica si el algoritmo ha funcionado o no
c                                 ierr=0   ! Todo ha ido bien
c                                 ierr=1   ! Numero maximo de iteraciones y restarts alcanzado
c                                 ierr=2   ! Algoritmo estancado
c     cputime=[Salida]  (REAL*8): Variable que indica el tiempo consumido en el calculo
c                                 de la solucion del S.E.L.
c
c------------------------------------------------------------------------------

      subroutine SLE_Solver_GMRES_CP(ntotv,ntotr1,ntotr2,nvfix,
     &                          niter,nrestart,nactu,epslr,rltvr,
     &                          lp1,lp2,rgmat_l,rgmat_u,gload,xdisp,
     &                          resdl,
     &                          ngfix,presc,treac,
     &                          wrksp,
     &                          residuo,kiter,irestart,ierr,cputime)
      implicit integer*8 (i-n)
      implicit real*8 (a-h,o-z)
      dimension wrksp((ntotv+4)*(niter+3))
!
!  Asignamos posiciones de memoria para cada variable auxiliar del GMRES
!
      i_v=1
      i_z=i_v+(ntotv+1)
      i_y=i_z+(ntotv+1)
      i_xh=i_y+(niter+1)
      i_ang=i_xh+(ntotv+1)*(niter+1)

      call SLE_Solver_GMRES_PD_CP(ntotv,ntotr1,ntotr2,nvfix,
     &                         niter,nrestart,nactu,epslr,rltvr,
     &                         lp1,lp2,rgmat_l,rgmat_u,gload,xdisp,
     &                         resdl,
     &                         ngfix,presc,treac,
     &                         wrksp(i_v),wrksp(i_z),
     &                         wrksp(i_y),wrksp(i_xh),wrksp(i_ang),
     &                         residuo,kiter,irestart,ierr,cputime)
     
      return
      end

c
c------------------------------------------------------------------------------
c
c   VARIABLES AUXILIARES:
c
c     vvect = [ Aux.  ] (REAL*8): Vector que almacena los distintos vectores v
c                                 que se van calculando en el Metodo GMRES
c     zvect = [ Aux.  ] (REAL*8): Vector auxiliar que se utiliza para almacenar
c                                 resultados intermedios a lo largo del proceso
c     y     = [ Aux.  ] (REAL*8): Vector de terminos independientes del sistema
c                                 triangular resultante de aplicar triangularizacion
c                                 a la matriz de Hessemberg del problema.
c                                 Como resultado se obtiene la solucion del sistema
c                                 triangular, necesaria para calcular la solucion
c                                 del sistema de ecuaciones original.
c     xh    = [ Aux.  ] (REAL*8): Matriz de dimension ((ntotv+1),(niter+1)) en la
c                                 que se almacena en la parte triangular superior los
c                                 coeficientes de la matriz de Hessemberg original y
c                                 posteriormente los coeficientes obtenidos tras aplicar
c                                 triangularizacion.
c                                 La parte inferior de la misma almacena los vectores
c                                 w necesarios para aplicar los proyectores de Householder.
c                                 En el coeficiente de la diagonal se almacena inicialmente
c                                 el ultimo termino de cada vector h que se anula tras
c                                 cada triangularizacion. Asi, al final el termino de la
c                                 diagonal puede almacenar el primer coeficiente del
c                                 vector w.
c     ang_tri=[ Aux.  ] (REAL*8): Vector que almacena el angulo que es necesario
c                                 aplicar para triangular la matriz de Hessemberg.
c                                 Se almacenan todos desde el principio porque a medida
c                                 que se aumenta la dimension del subespacio es preciso
c                                 aplicar todas las triangularizaciones anteriores al nuevo
c                                 vector h que se obtenga (la nueva columna de la matriz de
c                                 Hessemberg).

      subroutine SLE_Solver_GMRES_PD_CP(ntotv,ntotr1,ntotr2,nvfix,
     &                            niter,nrestart,nactu,epslr,rltvr,
     &                            lp1,lp2,rgmat_l,rgmat_u,gload,xdisp,
     &                            resdl,
     &                            ngfix,presc,treac,
     &                            vvect,zvect,
     &                            y,xh,ang_tri,
     &                            residuo,kiter,irestart,ierr,cputime)
      implicit integer*8 (i-n)
      implicit real*8 (a-h,o-z)
      dimension lp1(ntotv),lp2(ntotv),rgmat_l(ntotr1),rgmat_u(ntotr2),
     &          gload(ntotv),
     &          vvect(ntotv+1),zvect(ntotv+1),
     &          y(niter+1),xh(ntotv+1,niter+1),
     &          ang_tri(2*(niter+1)),xdisp(ntotv),resdl(ntotv),
     &          ngfix(nvfix),presc(nvfix),treac(nvfix)

      logical cnvgc
!
! Inicializamos el contador de tiempo
!      
      call cpu_time(time)
      
      nuno=1 ! Lo utilizamos para imponer que el residuo es nulo en las variables prescritas
      ncero=0 ! Lo utilizamos para indicar que no se tienen en cuenta las vbles. prescritas
!
!  Comprobamos que no hay diagonales nulas para utilizar precondicionamiento
!
      call SLE_CompruebaPivots(ntotv,ntotr1,ntotr2,rgmat_l,rgmat_u,
     &                         lp1,lp2)
!
!   Desmarcamos las variables prescritas
!
      call SLE_MarcarCoacciones(ntotv,nvfix,lp1,lp2,ngfix)
!
!   Imponemos los valores de las variables prescritas
!
      call SLE_ValoresPrescritos(ntotv,nvfix,ngfix,xdisp,presc)
!
!   Calculamos el vector residuo.
!
      call SLE_CalcularResiduo(nuno,ntotv,ntotr1,ntotr2,nvfix,lp1,lp2,
     &               rgmat_l,rgmat_u,xdisp,resdl,ngfix,presc,gload)

      call SLE_ResiduosPrescritos(ntotv,nvfix,ngfix,resdl)
    
! Calculamos la norma del residuo en el primer paso

      call SLE_Norma2(ntotv,resdl,residuoini)
      
      write(6,*)' Residuo inicial',(residuoini)

      call SLE_ActVector(nuno,residuoini,residuo) ! Lo definimos para irlo actualizando

! Calculamos la norma del vector de terminos independientes

      call SLE_Norma2(ntotv,gload,bnorma)
!
!  Arrancamos el bucle en restarts.
!
      do irestart=0,nrestart
      
        if (irestart.ge.1)then
          write(6,'(/a,i3/)')'Calculo reiniciado. Restart no.:',irestart
        endif
!
!  Inicializacion de las variables auxiliares
!
      call SLE_IniVector(ntotv+1,zvect)
      call SLE_IniVector(ntotv+1,vvect)
      call SLE_IniVector(2*(niter+1),ang_tri)
      call SLE_IniVector((ntotv+1)*(niter+1),xh)
      call SLE_IniVector(niter+1,y)
!
!  Copiamos en la primera columna de la matriz de Hessember el vector residuo
!
      call SLE_ActVector(ntotv,resdl,xh(1,1))
!
!  Inicializamos variables para cada restart
!
      ierr=-1
      cnvgc=.FALSE.
      kiter=0
!
!
!   Comienza el bucle do while en iteraciones (ampliaciones del subespacio de Krylov
!
!
      do while ((kiter.lt.niter+1).and.(.not.cnvgc))

        kiter=kiter+1
!
!  Guardamos el valor de la columna kiter de la matriz de Hessemberg en el vector z (Idem. Saad 2000)
!
        call SLE_ActVector(ntotv,xh(1,kiter),zvect)
!
! Almacenamos en beta el coeficiente de la diagonal de w porque de momento ese coeficiente
! de la matriz de Hessemberg aun no es nulo porque no hemos aplicado triangularizacion
!
        call SLE_Calcula_beta(ntotv,niter,kiter,xh,beta)
!
!  Normalizamos el vector w obtenido (que se almacena en la parte inferior de xh y en la
!  variable beta temporalmente
!
        call SLE_Normaliza_w(ntotv,niter,kiter,beta,xh)
!
!  Una vez calculado el vector del Proyector de Householder se aplica la transformacion al vector z
!  y se obtiene una nueva columna de la Matriz de Hessemberg. En la practica los coeficientes de la
!  la matriz no se modifican salvo el termino de la diagonal que es el unico que hay que calcular.
!  La matriz de Hessemberg se almacena en la parte superior de la matriz xh.
!
        call SLE_Calcula_h(ntotv,niter,kiter,beta,xh,zvect)  ! La parte inferior de xh y beta --> w
                                                             ! La parte superior de xh --> h
!
!  Si estamos en la primera iteracion almacenamos el coeficiente xh(1,1) en el vector de terminos
!  independientes y (que luego modificaremos de acuerdo con las triangularizaciones de la matriz
!  de Hessemberg.
!
        if (kiter.eq.1) y(1)=xh(1,1)   ! Guardamos el valor del residuo obtenido
!
!  En la primera iteracion no se triangulariza nada porque la matriz de Hessemberg comienza en la
!  columna 2 de xh. A partir de la primera se aplican las triangularizaciones anteriores y la ultima
!  (que se define como aquella que anula el coeficiente de la diagonal (Yousef Saad 2000)
!  La aplicacion de las operaciones de triangularizacion sobre el vector de terminos independientes
!  nos arroja como resultado y(kiter), cuyo valor absoluto proporciona el residuo en esa iteracion.
!
        if (kiter.gt.1)then
          call SLE_triangulariza(ntotv,niter,kiter,xh,ang_tri,
     &                         y,residuo)
        endif
!
!  Una vez triangularizada la columna de la matriz de Hessemberg nos queda libre el termino de la
!  diagonal de xh (que hemos anulado con la triangularizacion) y en esa posicion almacenamos la
!  primera componente del vector w. Las siguientes componentes no nulas estaran a continuacion de esta.
!
        call SLE_almacena_beta(ntotv,niter,kiter,xh,beta)
!
!  Si el numero de iteraciones es el adecuado imprimimos por pantalla el residuo actual.
!
        if (mod(dfloat(kiter-1),dfloat(nactu)).eq.0)then
          write(6,'(a,i6,a,e15.6)')'Iteracion: ',kiter-1,
     &                         ' ---> Residuo: ',dabs(residuo)
        endif
!
!  Si el residuo cumple la condicion de convergencia finalizamos el método GMRES
!
        if ((dabs(residuo).lt.dabs(epslr)).and.
     &      (dabs(residuo).lt.dabs(rltvr*residuoini)).and.
!     &      (dabs(residuo).lt.dabs(rltvr*bnorma)).and.
     &      (residuo.ne.0.d+00)    ! Lo hacemos para reiniciar cuando iter=ntotv
     &       )then
          cnvgc=.true.
          ierr=0
          if (mod(dfloat(kiter-1),dfloat(nactu)).ne.0)then ! Para evitar escribirlo 2 veces
            write(6,'(a,i6,a,e15.6)')'Iteracion: ',kiter-1,
     &                           ' ---> Residuo: ',dabs(residuo)
          endif
          write(6,'(/a/)')' Se ha alcanzado convergencia'
          exit
        elseif(kiter.eq.ntotv+1)then
          cnvgc=.true.    ! Lo establecemos para que se detenga el algoritmo
          ierr=3
        endif
!
!  Si el angulo de la triangularizacion es nulo quiere decir que en esa iteracion no hemos
!  conseguido reducir el residuo --> El metodo se estanca (stagnation). Esto puede ocurrir
!  si se utiliza GMRES con restart
!
        if (kiter.gt.1)then
          if (ang_tri(2*kiter-3).eq.0.d+00.and.ierr.ne.0)then
            ierr=2  ! El metodo se estanca
            write(6,*)' El metodo se estanca'
            exit
          endif
        endif
!
!  A continuacion comprobamos si se alcanza el numero maximo de iteraciones.
!
        if (kiter.eq.(niter+1).and.ierr.ne.0)then
          ierr=3     ! Reiniciamos
          write(6,'(a/)')' Numero de iteraciones alcanzado'
          if (irestart.eq.nrestart)ierr=1
          exit
        endif
!
!  Si no hemos llegado a la ultima iteracion, calculamos el vector ortornormal del subespacio de Krylov.
!
        if (kiter.lt.niter+1)then

          call SLE_calcula_v(ntotv,niter,kiter,xh,vvect)
!
!  Y con este vector ortonormal obtenemos el nuevo vector h que aumenta en una dimension
!  la matriz de Hessemberg antes de aplicar el proyector de Householder
!
          call SLE_calcula_z(ntotv,ntotr1,ntotr2,nvfix,niter,kiter,xh,
     &               vvect,zvect,
     &               rgmat_l,rgmat_u,lp1,lp2,ngfix,presc)
        endif

      enddo
!
!  Restamos en una unidad el indice de iteracion porque el metodo anterior se aplica hasta una dimension
!  menor que el numero de iteracion utilizada.
!
      kiter=kiter-1
!
!  Resolvemos el sistema triangular que minimiza el residuo. El sistema estara definido por la
!  matriz de Hessemberg ya triangularizada y por el vector y tambien modificado por las operaciones
!  de triangularizacion
!
      if (ierr.ne.-1)then   ! Lo imponemos por seguridad aunque no seria necesario comprobarlo
        call SLE_Uxb(ntotv,niter,kiter,xh,y)
        
        write(6,'(a,i8)')' kiter:',kiter
!
!  Una vez resuelto el sistema triangular se propone la nueva solucion. Si ierr=0 sera la solucion
!  del sistema de ecuaciones original. Si ierr=1 obtendremos una nueva solucion para reiniciar el
!  metodo GMRES.
!
        call SLE_calcula_solucion(ntotv,niter,nvfix,kiter,xh,y,
     &                          zvect,xdisp,ngfix)
      endif
!
!  Con la nueva solucion (si es que la hay) se calcula el nuevo residuo con el que reiniciar
!  el metodo si es necesario
!
      call SLE_CalcularResiduo(nuno,ntotv,ntotr1,ntotr2,nvfix,lp1,lp2,
     &                   rgmat_l,rgmat_u,xdisp,resdl,ngfix,presc,gload)
     
      call SLE_ResiduosPrescritos(ntotv,nvfix,ngfix,resdl)

      call SLE_Norma2(ntotv,resdl,residuo)
      
      write(6,'(a,e15.6)')' Residuo real ',residuo
      
      if (irestart.eq.nrestart.and.ierr.eq.1)then
        ierr=1
        write(6,'(/2a/)')'No. de iteraciones y de ',
     &                'restarts maximo alcanzado'
      endif
!
!  Si el metodo se estanca o si hemos alcanzado convergencia finalizamos el calculo
!
      if(ierr.eq.2.or.ierr.eq.0)exit
!
!   Finalizamos el bucle en restarts
!
      enddo    ! Fin del bucle en restarts
!
!  Corregimos en numero de irestart por si hemos alcanzado el maximo en el bucle en restarts.
!
      if (irestart.gt.nrestart)irestart=nrestart
      
!  Calculamos el residuo real con todas las ecuaciones. Asi en el vector resdl
!  obtendremos los valores de las reacciones que almacenamos en treac.
!
      call SLE_CalcularResiduo(ncero,ntotv,ntotr1,ntotr2,nvfix,lp1,lp2,
     &                  rgmat_l,rgmat_u,xdisp,resdl,ngfix,presc,gload)
!
!  Calculamos las reacciones
!
      call SLE_CalculaReacciones(ntotv,nvfix,ngfix,resdl,treac,gload)
!
!   Desmarcamos las variables prescritas
!
      call SLE_MarcarCoacciones(ntotv,nvfix,lp1,lp2,ngfix)
!
!  Obtenemos de nuevo informacion sobre el tiempo de calculo
!
      call cpu_time(cputime)

      cputime = cputime - time
!
!  Finalizamos el Metodo.
!
      return
      end

c_________________________________SUBRUTINA CompruebaPivots___________________
c
c     Comprueba que no hay ceros en la diagonal para aplicar precondicionamiento
c
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension del problema
c     nr1 = [Entrada]: numero de elementos no nulos de la matriz (inferior)
c     nr2 = [Entrada]: numero de elementos no nulos de la matriz (superior)
c     r   = [Entrada]: matriz del sistema de ecuaciones
c     lp1 = [Ent/Sal]: vector de posicion de fila de la matriz de rigidez
c     lp2 = [Ent/Sal]: vector de posicion de columna de la matriz de rigidez
c     
c     -------------------------------------------------------------------------

      subroutine SLE_CompruebaPivots(n,nr1,nr2,r_l,r_u,lp1,lp2)
      implicit integer*8(i-n)
      implicit real*8 (a-h,o-z)
      dimension  lp1(n),lp2(n),r_l(nr1),r_u(nr2)
      
      do i=1,n
        if (r_l(lp1(i)).eq.0.d+00.or.r_u(lp2(i)).eq.0.d+00)then
          write(6,'(a,i9)')' Pivot nulo en ecuacion:',i
          stop
        endif
      enddo
      
      return
      end
      
c_________________________________SUBRUTINA MarcarCoacciones___________________
c
c     Cambia de signo el puntero de las componentes de la matriz de rigidez afectadas por
c     variables prescritas
c
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension del problema
c     nr  = [Entrada]: numero de elementos no nulos de la matriz
c     nv  = [Entrada]: numero de variables prescritas del problema
c     lp1 = [Ent/Sal]: vector de posicion de fila de la matriz de rigidez
c     lp2 = [Ent/Sal]: vector de posicion de columna de la matriz de rigidez
c     ng  = [Entrada]: vector que indica el numero de variable prescrita
c     
c     --------------------------------------------------------------------

      subroutine SLE_MarcarCoacciones(n,nv,lp1,lp2,ng)
      implicit integer*8(i-n)
      implicit real*8 (a-h,o-z)
      dimension  lp1(n),lp2(n),ng(nv)

c.....bucle sobre coacciones
      do iv=1,nv
        lp1(ng(iv))=-lp1(ng(iv))
        lp2(ng(iv))=-lp2(ng(iv))
      enddo

      return
      end
      
c_________________________________SUBRUTINA ValoresPrescritos___________________
c
c     Impone de forma directa los valores prescritos de las variables
c
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension del problema
c     nv  = [Entrada]: numero de variables prescritas del problema
c     ng  = [Entrada]: vector que indica el numero de variable prescrita
c     x   = [Ent/Sal]: vector de variables del problema
c     pr  = [Ent/Sal]: vector que contiene los valores de las variables prescritas
c     
c     -----------------------------------------------------------------------

      subroutine SLE_ValoresPrescritos(n,nv,ng,x,pr)
      implicit integer*8(i-n)
      implicit real*8 (a-h,o-z)
      dimension ng(nv),x(n),pr(nv)
      
      do i=1,nv
        x(ng(i))=pr(i)
      enddo
      
      return
      end
      
c______________________________________SUBRUTINA CalcularResiduo_______________
c                           _   _     _
c     Calculo del residuo   r = b - A x
c                                   ~
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension de la matriz A
c     nr1 = [Entrada]: numero de terminos no nulos de la matriz del sistema (parte inferior)
c     nr2 = [Entrada]: numero de terminos no nulos de la matriz del sistema (parte superior)
c     nv  = [Entrada]: numero de variables prescritas del problema
c     lp1 = [Entrada]: vector de indice de fila de la matriz del sistema
c     lp2 = [Entrada]: vector de indice de columna de la matriz del sistema
c     r_l = [Entrada]: matriz del sistema en perfil por columnas (parte inferior)
c     r_u = [Entrada]: matriz del sistema en perfil por columnas (parte superior)
c     v   = [Entrada]: vector por el que se multiplica la matriz del sistema
c     res = [Entrada]: vector residuo
c     ng  = [Entrada]: vector que indica el numero de variable prescrita
c     pr  = [Ent/Sal]: vector que contiene los valores de las variables prescritas
c     b   = [Entrada]: vector de terminos independientes
c     -------------------------------------------------------------------------

      subroutine SLE_CalcularResiduo(kniv,n,nr1,nr2,nv,lp1,lp2,r_l,r_u,
     &                               v,res,ng,pr,b)
      implicit integer*8 (i-n)
      implicit real*8 (a-h,o-z)
      dimension res(n),b(n),ng(nv),pr(nv),v(n)
!
!  Realizamos el producto matriz vector
!
      call SLE_MatrizVector(kniv,n,nr1,nr2,nv,lp1,lp2,r_l,r_u,
     &                      v,ng,pr,res)
!
!  Calculamos el vector residuo
!
      do i=1,n
        res(i) = b(i) - res(i)
      enddo
!
!  Aplicamos el precondicionamiento diagonal
!
      if (kniv.eq.1)then
        call SLE_Precondicionamiento(n,nr1,nr2,lp1,lp2,r_l,r_u,res)
      endif

      return
      end

c______________________________________SUBRUTINA MatrizVector__________________
c
c     Realiza el producto de la matriz del sistema por un vector.
c     En este caso la matriz esta almacenada en perfil por columnas
c
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension del problema
c     nr1 = [entrada]: numero de elementos no nulos de la matriz (parte inferior)
c     nr2 = [entrada]: numero de elementos no nulos de la matriz (parte superior)
c     nv  = [Entrada]: numero de variables prescritas del problema
c     lp1 = [entrada]: vector de posicion de fila de la matriz de rigidez
c     lp2 = [entrada]: vector de posicion de columna de la matriz de rigidez
c     r_l = [entrada]: matriz del sistema en perfil por columnas (parte inferior)
c     r_u = [entrada]: matriz del sistema en perfil por columnas (parte superior)
c     ng  = [Entrada]: vector que indica el numero de variable prescrita
c     pr  = [Ent/Sal]: vector que contiene los valores de las variables prescritas
c     v   = [Entrada]: vector dato
c     rxv = [Salida]:  vector resultado del producto
c     
c     -------------------------------------------------------------------------

      subroutine SLE_MatrizVector(kniv,n,nr1,nr2,nv,lp1,lp2,r_l,r_u,
     &                      v,ng,pr,rxv)
      implicit integer*8(i-n)
      implicit real*8 (a-h,o-z)
      dimension lp1(n),lp2(n),r_l(nr1),r_u(nr2),v(n),rxv(n),
     &          ng(nv),pr(nv)

      do i=1,n
        rxv(i)=0.d+00
      enddo
! Parte inferior
      rxv(1)=r_l(1)*v(1)
      do k=2,n
        do j=abs(lp1(k-1))+1-abs(lp1(k)),0
          rxv(k)=rxv(k)+r_l(abs(lp1(k))+j)*v(k+j)
        enddo
      enddo
! Parte superior
      do k=2,n
        do j=abs(lp2(k-1))+1-abs(lp2(k)),-1   ! En este caso no operamos con la diagonal
          rxv(k+j)=rxv(k+j)+r_u(abs(lp2(k))+j)*v(k)
        enddo
      enddo
      
      if(kniv.eq.1)then
        do i=1,nv
          rxv(ng(i))=v(ng(i))
        enddo
      endif

      return
      end

c______________________________________SUBRUTINA Precondicionamiento_______________
c                           _          _
c     Calculo del residuo   r = P^{-1} r
c                               ~
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension de la matriz A
c     nr1 = [Entrada]: numero de terminos no nulos de la parte inferior de la matriz
c     nr2 = [Entrada]: numero de terminos no nulos de la parte superior de la matriz
c     lp1 = [Entrada]: vector de indice de fila de la matriz del sistema
c     lp2 = [Entrada]: vector de indice de columna de la matriz del sistema
c     r_l = [Entrada]: matriz de rigidez en perfil por columnas (parte inferior)
c     r_u = [Entrada]: matriz de rigidez en perfil por columnas (parte superior)
c     v   = [Entrada]: vector por el que se multiplica la matriz del sistema
c     -------------------------------------------------------------------------

      subroutine SLE_Precondicionamiento(n,nr1,nr2,lp1,lp2,r_l,r_u,v)
      implicit integer*8 (i-n)
      implicit real*8 (a-h,o-z)
      dimension lp1(n),lp2(n),r_l(nr1),r_u(nr2),v(n)

      do k=1,n
        if(lp1(k).gt.0)then
          v(k)=v(k)/r_l(lp1(k))
        endif
      enddo

      return
      end

c_________________________________SUBRUTINA ResiduosPrescritos___________________
c
c     Impone de forma directa los valores prescritos del residuo
c
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension del problema
c     nv  = [Entrada]: numero de variables prescritas del problema
c     ng  = [Entrada]: vector que indica el numero de variable prescrita
c     res = [Ent/Sal]: vector residuo del problema
c     
c     -------------------------------------------------------------------------

      subroutine SLE_ResiduosPrescritos(n,nv,ng,res)
      implicit integer*8(i-n)
      implicit real*8 (a-h,o-z)
      dimension ng(nv),res(n)
      
      do i=1,nv
        res(ng(i))=0.d+00
      enddo
      
      return
      end   


c______________________________________SUBRUTINA IniVector_____________
c                          _   _
c     Inicializa un vector v = 0
c
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension del problema
c     v   = [Salida] : vector inicializado a cero

c     -------------------------------------------------------------------------

      subroutine SLE_IniVector(n,v)
      implicit integer*8 (i-n)
      implicit real*8 (a-h,o-z)
      dimension v(n)

      cero=0.d+00

      do i=1,n
        v(i)=cero
      enddo

      return
      end
      
c______________________________________SUBRUTINA IniVector_____________
c                          
c     Copia un vector en otro
c
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension del problema
c     v   = [Entrada]: vector dato
c     c   = [Salida] : vector copia 

c     -------------------------------------------------------------------------

      subroutine SLE_ActVector(n,v,c)
      implicit integer*8 (i-n)
      implicit real*8 (a-h,o-z)
      dimension v(n), c(n)

      do i=1,n
        c(i)=v(i)
      enddo

      return
      end
      
c______________________________________SUBRUTINA Norma2_______________________
c
c     Calcula la norma 2 de un vector
c
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension del vector
c     v   = [Entrada]: vector v del que se calcula la norma
c     vn  = [Salida] : norma dos del vector v

c     -------------------------------------------------------------------
            
      subroutine SLE_Norma2(n,v,vmod)
      implicit integer*8(i-n)
      implicit real*8 (a-h,o-z)
      dimension v(n)
      
      vmod=0.d+00
      
      do i=1,n
        vmod=vmod+v(i)*v(i)
      enddo
      
      if (vmod.gt.0.d+00)then
        vmod=dsqrt(vmod)
      elseif(vmod.eq.0.d+00)then
        vmod=0.d+00 ! Mantenemos esta opcion porque cuando esto ocurra
                    ! el vector sera nulo y no se operara en las restantes
                    ! partes del algoritmo
      else
        write(6,*)' Vector nulo'
        stop
      endif
      
      return
      end

c___________________________SUBRUTINA Calcula_beta_______________________
c
c     Calcula la norma 2 de un vector
c
c     -------------------------------------------------------------------
c     n     = [Entrada]: Dimension del vector
c     nit   = [Entrada]: Dimension maxima del subespacio de Krylov utilizado
c     kiter = [Entrada]: Dimension actual del subespacio de Krylov
c     xh    = [Ent/Sal]: Matriz de Hessemberg y vectores w de proyeccion de Householder
c     b     = [Salida] : Variable en la que se almacena la w(kiter)

c     -------------------------------------------------------------------

      subroutine SLE_Calcula_beta(n,nit,kiter,xh,b)
      implicit integer*8(i-n)
      implicit real*8 (a-h,o-z)
      dimension xh(n+1,nit+1)

      vn=0.d+00

      do i=kiter,n
        vn = vn + xh(i,kiter) * xh(i,kiter)
      enddo

      if (vn.gt.0.d+00)then
        b=dsqrt(vn)
      elseif(vn.eq.0.d+00)then
        b=0.d+00
      else
        write(6,*)' ERROR: Vector negativo 2'
        stop
      endif

      return
      end

c______________________________________SUBRUTINA Normaliza_w_______________________
c
c     Normaliza el vector w del algoritmo de Householder que se encuentra almacenado
c     en la parte inferior de "xh" y con la diagonal en "b"
c
c     -------------------------------------------------------------------------
c     n     = [Entrada]: dimension del vector
c     nit   = [Entrada]: Dimension maxima del subespacio de Krylov utilizado
c     kiter = [Entrada]: Dimension actual del subespacio de Krylov
c     b     = [Salida] : Variable en la que se almacena la w(kiter)
c     xh    = [Ent/Sal]: Matriz de Hessemberg y vectores w de proyeccion de Householder

c     -------------------------------------------------------------------

      subroutine SLE_Normaliza_w(n,nit,kiter,b,xh)
      implicit integer*8 (i-n)
      implicit real*8(a-h,o-z)
      dimension xh(n+1,nit+1)

      xmod=2.d+00*b*b+2.d+00*b*abs(xh(kiter,kiter))  ! Es equivalente a la parte comentada mas abajo

      if(xh(kiter,kiter).gt.0.d+00)then
        b=+b+xh(kiter,kiter)
      else
        b=-b+xh(kiter,kiter)
      endif

!      xmod=0.d+00
!      do i=kiter+1,n      ! El termino de la diagonal esta almacenado en b
!        xmod=xmod+xh(i,kiter)*xh(i,kiter)
!      enddo
!      
!      xmod=xmod + b*b
      
      if (xmod.gt.0.d+00)then
        xmod=dsqrt(xmod)
      elseif(xmod.eq.0.d+00)then
        continue
      else
        write(6,*)' ERROR: Vector nulo o negativo'
        stop
      endif
      
      if (xmod.gt.0.d+00)then
        b=b/xmod
        do i=kiter+1,n
          xh(i,kiter)=xh(i,kiter)/xmod
        enddo
      endif
      
      return
      end

c______________________________________SUBRUTINA Calcula_h_______________________
c
c     Calculamos el valor de h a partir del valor de z y del valor de w almacenado 
c     en la parte inferior de la matriz xh y en beta. El resultado se almacena en la
c     parte superior de xh
c
c     -------------------------------------------------------------------------
c     n     = [Entrada]: dimension del sistema de ecuaciones
c     nit   = [entrada]: dimension del subespacio de krylov maximo
c     kiter = [entrada]: Indice de iteracion
c     xh    = [Ent/Sal]: matriz que almacena h en la parte superior y w en la inferior
c     b     = [Entrada]: vector que almacena las diagonales de w
c     z     = [entrada]: vector que almacena el vector z antes de aplicar Householder.
c                        Una vez aplicado se almacena en la parte superior de xh.

c     -------------------------------------------------------------------
      subroutine SLE_Calcula_h(n,nit,kiter,b,xh,z)
      implicit integer*8 (i-n)
      implicit real*8(a-h,o-z)
      dimension xh(n+1,nit+1),z(n+1)
!
!  Calculamos el valor de sigma (Yousef Saad 2000)
!
      sigma=b*z(kiter)
      do i=kiter+1,n
        sigma=sigma+xh(i,kiter)*z(i)
      enddo
!
!  Actualizamos la componente de la diagonal, que es la unica que cambia
!
      xh(kiter,kiter)=z(kiter)-2.d+00*sigma*b 
      
      return
      end      
      
c______________________________________SUBRUTINA Calcula_v_______________________
c
c     Calculamos el valor de v a partir de todos los w almacenados en la parte
c     inferior de la matriz xh y en beta
c
c     -------------------------------------------------------------------------
c     n     = [Entrada]: dimension del sistema de ecuaciones
c     nit   = [entrada]: dimension del subespacio de krylov maximo
c     kiter = [entrada]: Indice de iteracion
c     xh    = [Ent/Sal]: matriz que almacena h en la parte superior y w en la inferior
c     v     = [Ent/Sal]: vector que almacena el resultado de v

c     -------------------------------------------------------------------
      subroutine SLE_calcula_v(n,nit,kiter,xh,v)
      implicit integer*8 (i-n)
      implicit real*8(a-h,o-z)
      dimension xh(n+1,nit+1),v(n)
      
!  Inicializamos el vector v

      do i=1,n
        v(i)=0.d+00
      enddo
      
      v(kiter)=1.d+00
      
!  Aplicamos P_1 P_2...P_j e_j de forma recurrente

      do i=kiter,1,-1
        
! Aqui calculariamos ahora P_j·v de forma recurrente

        sigma=0.d+00
        do ik=i,n
          sigma=sigma+v(ik)*xh(ik,i)
        enddo
        
        do ik=i,n
          v(ik)=v(ik)-2.d+00*sigma*xh(ik,i)
        enddo

      enddo

      return
      end

c______________________________________SUBRUTINA Calcula_z_______________________
c
c     Calculamos el valor de z que se almacenara como una nueva columna de la matriz
c     xh. Se partira de "v" de rgmat a partir de todos los w almacenados en la parte
c     inferior de la matriz xh y en beta
c
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension del sistema de ecuaciones
c     nr1 = [entrada]: Numero de terminos no nulos de la matriz r (parte inferior)
c     nr2 = [entrada]: Numero de terminos no nulos de la matriz r (parte superior)
c     nv  = [Entrada]: numero de variables prescritas del problema
c     nit = [Entrada]: Dimension del subespacio de krylov utilizado
c     kiter = [Entrada]: Dimension actual del subespacio de Krylov
c     xh  = [Ent/Sal]: Matriz de Hessemberg en la parte superior y vector w en la inferior
c     v   = [Ent/Sal]: vector que almacena el resultado de v
c     w   = [Auxilia]: vector auxiliar que almacena w a partir de xh y beta
c     r_l = [Entrada]: vector que almacena una matriz en perfil por columnas (parte inferior)
c     r_u = [Entrada]: vector que almacena una matriz en perfil por columnas (parte superior)
c     lp1 = [Entrada]: Vector que indica la fila de los elementos almacenados en r
c     lp2 = [Entrada]: Vector que indica la columna de los elementos almacenados en r
c     ng  = [Entrada]: vector que indica el numero de variable prescrita
c     pr  = [Ent/Sal]: vector que contiene los valores de las variables prescritas
c     -------------------------------------------------------------------

      subroutine SLE_calcula_z(n,nr1,nr2,nv,nit,kiter,xh,v,w,
     &                     r_l,r_u,lp1,lp2,ng,pr)     
      implicit integer*8 (i-n)
      implicit real*8(a-h,o-z)
      dimension xh(n+1,nit+1),v(n),w(n),ng(nv),pr(nv)
      
      nuno=1
!
! Hacemos el producto de la matriz del sistema por el vector v
! calculado anteriormente
!
      call SLE_MatrizVector(nuno,n,nr1,nr2,nv,lp1,lp2,r_l,r_u,
     &                      v,ng,pr,w)
!
! Aplicamos precondicionamiento diagonal
!
      call SLE_Precondicionamiento(n,nr1,nr2,lp1,lp2,r_l,r_u,w)
!
! Actualizamos el vector w en el que hacemos el producto anterior
!
      call SLE_ActVector(n,w,xh(1,kiter+1)) 
      
      do i=1,kiter
!
! Calculamos el valor de sigma para aplicar el proyector de Householder
!
        sigma=0.d+00
      
        do ik=i,n
          sigma=sigma+xh(ik,kiter+1)*xh(ik,i)
        enddo
!
! Una vez calculado sigma aplicamos la transformacion de Householder
! para calcular z que se almacena provisionalmente sobre la siguiente
! columna de la matriz xh
!
        do ik=i,n
          xh(ik,kiter+1)=xh(ik,kiter+1)-2.d+00*sigma*xh(ik,i)
        enddo
        
      enddo

      return
      end

c______________________________________SUBRUTINA triangulariza________________
c
c    
c
c     -------------------------------------------------------------------------
c     n     = [Entrada]: dimension del sistema de ecuaciones
c     nit   = [Entrada]: dimension del subespacio
c     Kiter = [Entrada]: Numero de iteracion (tamano del subespacio)
c     xh    = [Ent/Sal]: Matriz de Hessemberg en la parte superior
c                        que vamos a triangularizar
c     ang_tri=[Ent/Sal]: vector que almacena el angulo con el que
c                        triangularizamos en cada iteración
c     y     = [Entrada]: vector que almacena el termino independiente
c                        del sistema que estamos triangularizando.
c     residuo= [Salida]: variable que almacena la norma del residuo
c                        obtenida como resultado de aplicar la triangularizacion
c                        de la matriz de Hessemberg

c     ------------------------------------------------------------------------

      subroutine SLE_triangulariza(n,nit,kiter,xh,ang_tri,
     &           y,residuo)
      implicit integer*8 (i-n)
      implicit real*8(a-h,o-z)
      parameter (pi=3.141592653589793238d+00)
      dimension xh(n+1,nit+1),ang_tri(2,nit+1),y(nit+1) ! La ultima componente la
                                                      ! utilizamos para el residuo
!
!  Triangularizamos la nueva columna hasta la ultima triangularizacion
!
      do i=1,kiter-2
        aux1= ang_tri(1,i)*xh(i  ,kiter)+
     &        ang_tri(2,i)*xh(i+1,kiter)
        aux2=-ang_tri(2,i)*xh(i  ,kiter)+
     &        ang_tri(1,i)*xh(i+1,kiter)
        xh(i  ,kiter)=aux1
        xh(i+1,kiter)=aux2
      enddo
!
!  Calculamos el nuevo angulo para aplicar la triangularizacion
!
      ang_tri(1,kiter-1)=xh(kiter-1,kiter)/      ! esto es el coseno
     &      dsqrt(xh(kiter-1,kiter)**2+xh(kiter,kiter)**2)
      ang_tri(2,kiter-1)=xh(kiter,kiter)/    ! esto es el seno
     &      dsqrt(xh(kiter-1,kiter)**2+xh(kiter,kiter)**2)
!
!  Triangularizamos la ultima posicion con el nuevo angulo calculado
!
      i=kiter-1
      aux1= ang_tri(1,i)*xh(i  ,kiter)+
     &      ang_tri(2,i)*xh(i+1,kiter)
      aux2=-ang_tri(2,i)*xh(i  ,kiter)+
     &      ang_tri(1,i)*xh(i+1,kiter)
      xh(i  ,kiter)=aux1
      xh(i+1,kiter)=aux2
!
!  Triangularizamos el vector de terminos independientes
!
      aux1= ang_tri(1,kiter-1)*y(kiter-1)+
     &      ang_tri(2,kiter-1)*y(kiter)
      aux2=-ang_tri(2,kiter-1)*y(kiter-1)+
     &      ang_tri(1,kiter-1)*y(kiter)
      y(kiter-1) = aux1
      y(kiter)   = aux2
      
      residuo=aux2   !   Actualizamos el residuo con el valor obtenido

      return
      end
c_____________________________SUBRUTINA SLE_almacena_beta______________________
c
c
c     -------------------------------------------------------------------------
c     n    = [Entrada]: dimension del sistema de ecuaciones
c     nit  = [Entrada]: dimension del subespacio
c     Kiter= [Entrada]: Numero de iteracion (tamano del subespacio)
c     xh   = [Ent/Sal]: Matriz de Hessemberg en la parte superior.
c                       La parte inferior incluida la diagonal almacenara
c                       los vectores w
c     b    = [Entrada]: valor que toma el termino kiter en la iteracion kiter
c                       del vector w

c     -------------------------------------------------------------------------

      subroutine SLE_almacena_beta(n,nit,kiter,xh,b)
      implicit integer*8 (i-n)
      implicit real*8(a-h,o-z)
      dimension xh(n+1,nit+1)
!
!   Guardamos el valor de b en la matriz xh una vez que hemos triangularizado
!
      xh(kiter,kiter)=b

      return
      end

c______________________________________SUBRUTINA Uxb___________________________
c
c    
c
c     -------------------------------------------------------------------------
c     n    = [Entrada]: dimension del sistema de ecuaciones
c     nit  = [Entrada]: dimension del subespacio
c     Kiter= [Entrada]: Numero de iteracion (tamano del subespacio)
c     xh   = [Ent/Sal]: Matriz de Hessemberg en la parte superior
c                       que vamos a triangularizar
c     y    = [Entrada]: vector que almacena el termino independiente
c                       del sistema que estamos triangularizando.

c     -------------------------------------------------------------------------

      subroutine SLE_Uxb(n,nit,kiter,xh,y)
      implicit integer*8 (i-n)
      implicit real*8(a-h,o-z)
      dimension xh(n+1,nit+1),y(nit+1)
!
!  Resolvemos un sistema con matriz triangular superior llena
!
      do i=kiter,1,-1
        y(i)=y(i)/xh(i,i+1)
        do j=i-1,1,-1
          y(j)=y(j)-xh(j,i+1)*y(i)
        enddo
      enddo
      
      return
      end
c____________________________SUBRUTINA Calcula_solucion_____________________
c
c
c
c     ----------------------------------------------------------------------
c     n    = [Entrada]: dimension del sistema de ecuaciones
c     nit  = [Entrada]: dimension del subespacio
c     nv  = [Entrada]: numero de variables prescritas del problema
c     Kiter= [Entrada]: Numero de iteracion (tamano del subespacio)
c     xh   = [Ent/Sal]: Matriz de Hessemberg en la parte superior
c                       que vamos a triangularizar
c     y    = [Entrada]: vector que almacena el termino independiente
c                       del sistema que estamos triangularizando.
c     z    = [ Aux.  ]: vector auxiliar que utilizamos para calcular
c                       la nueva solucion
c     x    = [ Aux.  ]: vector solucion que obtenemos tras aplicar
c                       el metodo
c     ng  = [Entrada]: vector que indica el numero de variable prescrita
c---------------------------------------------------------------------------

      subroutine SLE_calcula_solucion(n,nit,nv,kiter,xh,y,z,x,ng)
      implicit integer*8 (i-n)
      implicit real*8(a-h,o-z)
      dimension xh(n+1,nit+1),y(nit+1),z(n),x(n),ng(nv)

      call SLE_IniVector(n,z)
!
!  Aplicamos el metodo propuesto por Yousef Saad 2000 para calcular la nueva solucion
!
      do j=kiter,1,-1   ! Los w empiezan una columna antes que los h
      
        z(j)=z(j)+y(j)
                
        sigma=0.d+00
        
        do i=j,n
          sigma=sigma+xh(i,j)*z(i)
        enddo
        
        do i=j,n
          z(i)=z(i)-2.d+00*sigma*xh(i,j)
        enddo
        
      enddo
      
      do i=1,nv
        z(ng(i))=0.d+00
      enddo
      
      do i=1,n
        x(i)=x(i)+z(i)
      enddo
        
      return
      end
      
c____________________________SUBRUTINA CalculaReacciones_____________________
c
c
c
c     ----------------------------------------------------------------------
c     n    = [Entrada]: dimension del sistema de ecuaciones
c     nv  = [Entrada]: numero de variables prescritas del problema
c     ng  = [Entrada]: vector que indica el numero de variable prescrita
c     res = [Entrada]: vector residuo
c     treac=[Ent/Sal]: vector que contiene los valores de las reacciones

c---------------------------------------------------------------------------

      subroutine SLE_CalculaReacciones(n,nv,ng,res,tr,b)
      implicit integer*8 (i-n)
      implicit real*8(a-h,o-z)
      dimension ng(nv),res(n),tr(nv),b(n)

      do i=1,nv
        tr(i)=-res(ng(i))+b(ng(i))
      enddo
      
      return
      end
      