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

c______________________________________SUBRUTINA Solver_PCG_ER_________________
c
c AUTORES: F. Navarrina, I. Colominas, H. Gomez, M. Casteleiro
c                                                                 _
c        Dep. de Metodos Matematicos y de Rep., ETSICCP de La Coruna, UDC.
c
c FECHA ULTIMA REVISION: 20070320
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
c------------------------------------------------------------------------------
c
c OBJETIVO:
c
c    Solucion general del sistema de ecuaciones :
c
c                                         -   -   -
c                                     A . x = b + r
c                                     ~            v
c    con las condiciones adicionales
c                                     _     _
c                                     x   = p       
c                                      v     v
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           r = reacciones en gdl coaccionados                      [incognita]
c            v
c           p = valores prescritos de los gdl coaccionados          [dato]
c            v
c------------------------------------------------------------------------------
c 
c DESCRIPCION GENERAL:
c
c       La matriz A es simetrica y puede ser singular.
c
c       La diagonal principal de la matriz A se encuentra ensamblada, y se
c    almacena en los primeros elementos del vector rgmat de forma que
c
c          a            <--->  rgmat(lpii),   lpii=lpamc(itotv)=itotv
c           itotv,itotv                       itotv=1,...,ntotv
c       
c       El residuo se recalcula a nivel elemental cuando es necesario,
c    de forma que el componente del residuo del elemento ielem
c
c           ielem
c          r           ,                      ievab=1,...,nevab
c           ievab                             nevab=ndofn*nnode
c
c    corresponde al grado de libertad idofn del nodo inode, siendo
c
c           ievab=(inode-1)*ndofn+idofn,
c
c    y se ensambla en el coeficiente r      del residuo global,
c                                     itotv
c    donde
c
c           itotv=(ipoin-1)*ndofn+idofn, ipoin=lnods(inode,ielem),
c
c       Se asume que al imponerle al sistema las condiciones de vinculacion
c    (gdl coaccionados) se obtiene una matriz definida positiva (no singular).
c       
c       El sistema resultante A'. x = b' se resuelve mediante el metodo 
c    de Gradientes Conjugados con Precondicionamiento Diagonal. 
c
c       Simultaneamente se calculan las reacciones en los gdl coaccionados.
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     ntype = [Entrada]: codigo de tipo de calculo
c     nnode = [Entrada]: numero de nodos por elemento
c     nelem = [Entrada]: numero de elementos
c     nprop = [Entrada]: numero de propiedades
c     nmats = [Entrada]: numero de materiales
c     nstre = [Entrada]: numero de deformaciones/tensiones por punto
c     ndime = [Entrada]: numero de dimensiones del espacio material
c     ndofn = [Entrada]: numero de gdl por nodo
c     npoin = [Entrada]: numero de nodos de la malla
c     nevab = [Entrada]: numero de gdl por elemento = ndofn*nnode
c     ntotv = [Entrada]: numero de ecuaciones =ndofn*npoin 
c     lnods = [Entrada]: lnods(in,ie)->num. global del nodo in del elemento ie
c     coord = [Entrada]: coord(id,ip)->coord. id del nodo global ip
c     elcod = [Trabajo]: elcod(id,in)->coord. id del nodo local in
c     matno = [Entrada]: matno(ie)->numero de material del elemento i
c     props = [Entrada]: props(ip,im)->propiedad ip del material im
c     nvfix = [Entrada]: numero de gdl coaccionados.
c     ngfix = [Entrada]: lista de gdl coaccionados.
c     presc = [Entrada]: lista de valores prescritos en gdl coaccionados.
c     lpamc = [Entrada]: vector de punteros de los elementos de la diagonal
c                        principal de la matriz A en el vector rgmat
c                        (normalmente lpamc(itotv)=itotv).
c     rgmat = [Entrada]: vector de dimension ntotv en el que se encuentra
c                        almacenada la diagonal de la matriz A.
c     gload = [Entrada]: vector de dimension ntotv en el que se encuentran
c                        almacenados los terminos independientes.
c     eldis = [Trabajo]: vector de trabajo de nevab componentes (desplaz.)
c     eldef = [Trabajo]: vector de trabajo de nstre componentes (deformaciones)
c     elstr = [Trabajo]: vector de trabajo de nstre componentes (tensiones)
c     eload = [Trabajo]: vector de trabajo de nevab componentes (fuerzas)
c     elres = [Trabajo]: vector de trabajo de nevab componentes (residuos)
c     dmatx = [Trabajo]: matriz dmatx(nstre,nstre)
c     bmatx = [Trabajo]: matriz bmatx(nstre,nevab)
c     niter = [Entrada]: numero maximo de iteraciones
c     nactu = [Entrada]: cada cuantas iteraciones se recalcula el residuo
c                        (0=>nunca)
c     epslx = [Entrada]: error absoluto admisible en x  (norma del maximo)
c     rltvx = [Entrada]: error relativo admisible en x  (norma del maximo)
c     epslr = [Entrada]: error absoluto admisible en r' (norma del maximo)
c                        siendo r' el residuo           (r' = b' - A' x)
c     rltvr = [Entrada]: error relativo admisible en r' (norma del maximo)
c                        en relacion con b' 
c     wvect = [Trabajo]: vector de trabajo de ntotv componentes
c     svect = [Trabajo]: vector de trabajo de ntotv componentes
c     zvect = [Trabajo]: vector de trabajo de ntotv componentes
c     xdisp = [Entrada]: aproximacion inicial del vector solucion
c     xdisp = [Salida] : vector solucion
c     resdl = [Salida] : residuo r' = b' - A' x
c     treac = [Salida] : valores de las reacciones.
c     rnorm = [Salida] : norma del maximo del residuo r'
c     kiter = [Salida] : numero de iteraciones realizadas
c     cnvgc = [Salida] : ha convergido el algoritmo (.TRUE.) o no (.FALSE)
c
c------------------------------------------------------------------------------

      subroutine SLE$Solver_PCG_ER(ntype,nnode,nelem,nprop,nmats,
     .                             nstre,ndime,ndofn,npoin,
     .                             nevab,ntotv,
     .                             lnods,coord,elcod,matno,props,
     .                             nvfix,ngfix,presc,
     .                             lpamc,rgmat,gload,
     .                             eldis,eldef,elstr,eload,elres,
     .                             dmatx,bmatx,
     .                             niter,nactu,
     .                             epslx,rltvx,epslr,rltvr,
     .                             wvect,svect,zvect,
     .                             xdisp,resdl,treac,
     .                             rnorm,kiter,cnvgc)

d     implicit real*8 (a-h,o-z)

      logical cnvgc,recal

      cero=0.d+00
      uno =1.d+00

c.....Cambio de signo a punteros de gdl restringidos
      if (nvfix.gt.0) call SLE$MarcarCoacciones(ntotv,lpamc,nvfix,ngfix)

c.....Comprobacion
      call SLE$VerificarPivots(ntotv,lpamc,rgmat)

c.....Inicializacion
      call SLE$NormaMaximo(ntotv,lpamc,gload,gnorm)
      rnmax=gnorm*rltvr
      if (epslr.gt.rnmax) rnmax=epslr

      if (nvfix.gt.0)
     .  call SLE$ValoresPrescritos(nvfix,ngfix,presc,ntotv,xdisp)

      call SLE$CopiarVector(ntotv,gload,resdl)
      do ielem=1,nelem
        call CAL$Proc_Resdl(2    ,ntype,ielem,nelem,nnode,
     .                      npoin,ndime,ndofn,nstre,nmats,
     .                      nprop,nevab,ntotv,lnods,
     .                      coord,elcod,matno,props,
     .                      eldis,eldef,elstr,eload,elres,
     .                      dmatx,bmatx,xdisp,resdl)
      enddo
      call SLE$InicializarVector(ntotv,zvect)
      call SLE$InicializarVector(ntotv,svect)
      sas=uno
      cnvgc=.FALSE.
      kiter=0

c.....ITERACIONES
      do while ((kiter.lt.niter).and.(.not.cnvgc))

c.......Direccion conjugada (s)
        call SLE$PreconDiagonal(ntotv,lpamc,rgmat,resdl,wvect)
        call SLE$ProductoEscalar(ntotv,zvect,wvect,zw)
        beta=zw/sas
        call SLE$ActualizarVector(ntotv,wvect,beta,svect,svect)

c.......Factor de avance (alfa)
        call SLE$InicializarVector(ntotv,zvect)
        do ielem=1,nelem
          call CAL$Proc_Resdl(2    ,ntype,ielem,nelem,nnode,
     .                        npoin,ndime,ndofn,nstre,nmats,
     .                        nprop,nevab,ntotv,lnods,
     .                        coord,elcod,matno,props,
     .                        eldis,eldef,elstr,eload,elres,
     .                        dmatx,bmatx,svect,zvect)
        enddo
        call SLE$ProductoEscalar(ntotv,zvect,svect,sas)
        if (sas.eq.cero) sas=uno
        sas=-sas
        call SLE$ProductoEscalar(ntotv,resdl,svect,rs)
        alfa=rs/sas

c.......Actualizacion de la solucion (x)
        call SLE$NormaMaximo(ntotv,lpamc,svect,snorm)
        dx=abs(alfa)*snorm
        call SLE$NormaMaximo(ntotv,lpamc,xdisp,xnorm)
        dxmax=xnorm*rltvx
        if (epslx.gt.dxmax) dxmax=epslx
        cnvgc=(dx.le.dxmax)
        call SLE$ActualizarVector(ntotv,xdisp,alfa,svect,xdisp)

c.......Actualizacion del residuo (r)
        if (nactu.eq.0) then
          recal=.FALSE.
        else
          recal=(mod(kiter+1,nactu).eq.0)
        endif
        if (.not.recal) then
          call SLE$ActualizarVector(ntotv,resdl,alfa,zvect,resdl)
          call SLE$NormaMaximo(ntotv,lpamc,resdl,rnorm)
          cnvgc=(cnvgc.and.(rnorm.le.rnmax))
        endif
        if (recal.or.cnvgc) then
          call SLE$CopiarVector(ntotv,gload,resdl)
          do ielem=1,nelem
            call CAL$Proc_Resdl(2    ,ntype,ielem,nelem,nnode,
     .                          npoin,ndime,ndofn,nstre,nmats,
     .                          nprop,nevab,ntotv,lnods,
     .                          coord,elcod,matno,props,
     .                          eldis,eldef,elstr,eload,elres,
     .                          dmatx,bmatx,xdisp,resdl)
          enddo
          call SLE$NormaMaximo(ntotv,lpamc,resdl,rnorm)
          cnvgc=(cnvgc.and.(rnorm.le.rnmax))
        endif

        kiter=kiter+1
      enddo

c.....Reacciones
      if (nvfix.gt.0)
     .  call SLE$ObtenerReacciones(ntotv,resdl,nvfix,ngfix,treac)

c.....Restauracion del signo a punteros de gdl restringidos
      if (nvfix.gt.0) call SLE$MarcarCoacciones(ntotv,lpamc,nvfix,ngfix)

      return
      end


c______________________________________SUBRUTINA MarcarCoacciones______________
c
c     Cambio de signo de punteros de gdl restringidos
c
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension de la matriz A
c     lp  = [Entrada]: vector de punteros de los elementos de la diagonal.
c     nv  = [Entrada]: numero de gdl coaccionados.
c     ng  = [Entrada]: lista de gdl coaccionados.
c     -------------------------------------------------------------------------

      subroutine SLE$MarcarCoacciones(n,lp,nv,ng)
      dimension  lp(n),ng(nv)

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

      return
      end


c______________________________________SUBRUTINA VerificarPivots_______________
c
c     Verifica que los gdl no coaccionados tengan un pivot no nulo
c
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension de la matriz A
c     lp  = [Entrada]: vector de punteros de los elementos de la diagonal.
c                      (puntero negativo => el gdl esta coaccionado).
c     a   = [Entrada]: matriz de coeficientes (simetrica, semidefinida +)
c                      almacenada por filas en Column Profile.
c     -------------------------------------------------------------------------

      subroutine SLE$VerificarPivots(n,lp,a)

d     implicit real*8 (a-h,o-z)
      dimension lp(n),a(1)

      cero=0.d+00

      do i=1,n
        if (lp(i).gt.0) then
          lpii=iabs(lp(i))
          if (a(lpii).eq.cero) call SLE$FError(' pivot nulo en fila:',i)
        endif
      enddo

      return
      end


c______________________________________SUBRUTINA ValoresPrescritos_____________
c
c     Inicializa los valores prescritos
c
c     -------------------------------------------------------------------------
c     nv  = [Entrada]: numero de gdl coaccionados.
c     ng  = [Entrada]: lista de gdl coaccionados.
c     pr  = [Entrada]: lista de valores prescritos en los gdl coaccionados.
c     n   = [Entrada]: dimension del problema
c     x   = [Entrada]: aproximacion inicial al vector solucion
c           [Salida] : idem, inicializada con los valores prescritos
c     -------------------------------------------------------------------------

      subroutine SLE$ValoresPrescritos(nv,ng,pr,n,x)

d     implicit real*8 (a-h,o-z)
      dimension ng(nv),pr(nv),x(n)

      do iv=1,nv
        i=ng(iv)
        x(i)=pr(iv)
      enddo

      return
      end


c______________________________________SUBRUTINA ObtenerReacciones_____________
c
c     Obtiene las reacciones a partir del residuo
c
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension del problema
c     r   = [Entrada]: residuo
c     nv  = [Entrada]: numero de gdl coaccionados.
c     ng  = [Entrada]: lista de gdl coaccionados.
c     tr  = [Salida] : reacciones en los gdl coaccionados.
c     -------------------------------------------------------------------------

      subroutine SLE$ObtenerReacciones(n,r,nv,ng,tr)

d     implicit real*8 (a-h,o-z)
      dimension r(n),ng(nv),tr(nv)

      cero=0.d+00

      do iv=1,nv
        i=ng(iv)
        tr(iv)=-r(i)
        r(i)=cero
      enddo

      return
      end


c______________________________________SUBRUTINA InicializarVector_____________
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$InicializarVector(n,v)

d     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 CopiarVector__________________
c                     _   _
c     Copia un vector w = v
c
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension del problema
c     v   = [Entrada]: vector dato
c     w   = [Salida] : copia del vector anterior

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

      subroutine SLE$CopiarVector(n,v,w)

d     implicit real*8 (a-h,o-z)
      dimension v(n),w(n)


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

      return
      end


c______________________________________SUBRUTINA PreconDiagonal________________
c                                            _     -1 _
c     Precondiciona diagonalmente el residuo w' = D   r'
c                                                 ~A
c     ignorando los gdl coaccionados
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension de la matriz A
c     lp  = [Entrada]: vector de punteros de los elementos de la diagonal
c                      (puntero negativo => el gdl esta coaccionado).
c     a   = [Entrada]: matriz de coeficientes (simetrica, semidefinida +)
c                      almacenada por filas en Column Profile.
c     r   = [Entrada]: vector residuo 
c     w   = [Salida] : vector residuo precondicionado

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

      subroutine SLE$PreconDiagonal(n,lp,a,r,w)

d     implicit real*8 (a-h,o-z)
      dimension lp(n),a(1),r(n),w(n)

      cero=0.d+00

      do i=1,n
        if (lp(i).gt.0) then
          lpii=iabs(lp(i))
          w(i)=r(i)/a(lpii)
        else
          w(i)=cero
        endif
      enddo

      return
      end


c______________________________________SUBRUTINA ActualizarVector______________
c                                                _   _       _
c     Realiza actualizaciones vectoriales tipo   y = x + t * u
c
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension de los vectores
c     x   = [Entrada]: vector inicial x
c     t   = [Entrada]: factor de avance
c     u   = [Entrada]: vector de avance
c     y   = [Salida] : vector actualizado

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

      subroutine SLE$ActualizarVector(n,x,t,u,y)

d     implicit real*8 (a-h,o-z)
      dimension x(n),u(n),y(n)

      do i=1,n
        y(i)=x(i)+t*u(i)
      enddo

      return
      end


c______________________________________SUBRUTINA ProductoEscalar_______________
c                                                  _    _
c     Calcula el producto escalar de dos vectores, v  y v
c                                                   1    2
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension de los vectores
c     v1  = [Entrada]: vector v1
c     v2  = [Entrada]: vector v2
c     pe  = [Salida] : producto escalar

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

      subroutine SLE$ProductoEscalar(n,v1,v2,pe)

d     implicit real*8 (a-h,o-z)
      dimension v1(n),v2(n)

      cero=0.d+00

      pe=cero

      do i=1,n
        pe=pe+v1(i)*v2(i)
      enddo

      return
      end

c______________________________________SUBRUTINA NormaMaximo___________________
c
c     Calcula la norma del maximo de un vector, ignorando los gdl coaccionados
c
c     -------------------------------------------------------------------------
c     n   = [Entrada]: dimension del vector
c     lp  = [Entrada]: vector de punteros
c                      (puntero negativo => el gdl esta coaccionado).
c     v   = [Entrada]: vector v
c     vn  = [Salida] : norma del maximo del vector v

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

      subroutine SLE$NormaMaximo(n,lp,v,vn)

d     implicit real*8 (a-h,o-z)
      dimension lp(n),v(n)

      cero=0.d+00

      vn=cero

      do i=1,n
        if (lp(i).gt.0) then
          if (abs(v(i)).gt.vn) vn=abs(v(i))
        endif
      enddo

      return
      end


c______________________________________SUBRUTINA FError________________________

      subroutine SLE$FError(alfa,i)
      character alfa*(*)

      write(6,'(2a,i5)') ' SLE> ',alfa,i

      stop
      end
