c==============================================================================
c=============  LIBRERIA PARA CALCULO MATRICIAL DE ESTRUCTURAS  ===============
c==============================================================================

c______________________________________Libreria Libcmepcge_____________________
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: 20070410
c
c Copyright: Esta libreria es propiedad de su autor.
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===============================================================================
c================================== PRIMER NIVEL ===============================
c===============================================================================

c____________________________________________________________________CAL$Ini_Par
c     Inicializacion de parametros del programa

      Subroutine CAL$Ini_Par(ntype,ndime,npoin,ndofn,nnode,nstre,
     .           nprop,nevab,ntotv)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)
   
c.....Parametros que dependen del tipo de calculo
      if (ntype.eq.2) then        ! ESTRUCTURADAS ARTICULADAS 2D
        ndime=2                   ! 2 Coordenadas por punto        (x,y)
        ndofn=2                   ! 2 Grados de Libertad por Nodo  (u,v)
        nnode=2                   ! 2 Nodos por Elemento
        nstre=1                   ! 1 Esfuerzo                     (N)
                                  ! 1 Deformacion                  (dL)
        nprop=1                   ! 1 Propiedad por Elemento
                                  !      {
                                  !       1-> Rigidez a axil       (EA)
                                  !      }
      else if (ntype.eq.3) then   ! ESTRUCTURADAS ARTICULADAS 3D
        ndime=3                   ! 3 Coordenadas por punto        (x,y,z)
        ndofn=3                   ! 3 Grados de Libertad por Nodo  (u,v,w)
        nnode=2                   ! 2 Nodos por Elemento
        nstre=1                   ! 1 Esfuerzo                     (N)
                                  ! 1 Deformacion                  (dL)
        nprop=1                   ! 1 Propiedad por Elemento
                                  !      {
                                  !       1-> Rigidez a axil       (EA)
                                  !      }
      else if (ntype.eq.4) then   ! ESTRUCTURADAS RETICULADAS 2D
        ndime=2                   ! 2 Coordenadas por punto        (x,y)
        ndofn=3                   ! 1 Grado de Libertad por Nodo   (u,v,gz)
        nnode=2                   ! 2 Nodos por Elemento
        nstre=3                   ! 3 Esfuerzos                    (N,M1,M2)
                                  ! 3 Deformaciones                (dl,g1,g2)
        nprop=2                   ! 2 Propiedades por Elemento
                                  !      {
                                  !       1-> Rigidez a axil       (EA)
                                  !       2-> Rigidez a flexion    (EI)
                                  !      }
      else
        call SRV$Error('CAL$Ini_Par: Opcion no implementada (NTYPE)')
      endif

c.....Parametros de almacenamiento generales
      nevab=nnode*ndofn      ! Grados de libertad de elemento
      ntotv=npoin*ndofn      ! Grados de libertad del sistema

      return
      end

C____________________________________________________________________CAL$Ini_AMC
c     Obtencion de los punteros de los grados de libertad en el almacenamiento
c     por columnas con renumeracion opcional


      Subroutine CAL$Ini_AMC(nelem,nnode,npoin,ndofn,ntotv,lnods,lpamc,
     .           lpntn,nrigi)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension lnods(nnode,nelem),lpamc(ntotv),lpntn(npoin)

      do ipoir=1,npoin           ! Inicializacion
        itot0=(ipoir-1)*ndofn
        lpamc(itot0+1)=0
      enddo

      do ielem=1,nelem           ! Altura de las columnas del perfil
        do inode=1,nnode-1
          itot0=(lpntn(lnods(inode,ielem))-1)*ndofn
          do jnode=inode+1,nnode
            jtot0=(lpntn(lnods(jnode,ielem))-1)*ndofn
            if (itot0.gt.jtot0) then
              kdifr=(itot0-jtot0)
              if (kdifr.gt.lpamc(itot0+1)) lpamc(itot0+1)=kdifr
            else
              kdifr=(jtot0-itot0)
              if (kdifr.gt.lpamc(jtot0+1)) lpamc(jtot0+1)=kdifr
            endif
          enddo
        enddo
      enddo

      ipoir=1                    ! Acumulacion del Perfil
        itot0=(ipoir-1)*ndofn
        kdifr=lpamc(itot0+1)
        idofn=1
          itotv=itot0+idofn
          lpamc(itotv)=               kdifr+idofn
        do idofn=2,ndofn
          itotv=itot0+idofn
          lpamc(itotv)=lpamc(itotv-1)+kdifr+idofn
        enddo
      do ipoir=2,npoin
        itot0=(ipoir-1)*ndofn
        kdifr=lpamc(itot0+1)
        do idofn=1,ndofn
          itotv=itot0+idofn
          lpamc(itotv)=lpamc(itotv-1)+kdifr+idofn
        enddo
      enddo

      nrigi=lpamc(ntotv)         ! Posiciones de memoria que ocupa la matriz
      return
      end

C____________________________________________________________________CAL$Ini_AEM
c     Inicializacion del almacenamiento de la diagonal principal (ensamblada)
c     y de las partes triangulares inferiores de las matrices elementales
c     (sin ensamblar)

      Subroutine CAL$Ini_AEM(nelem,ntotv,nevab,lpamc,nrigi)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension lpamc(ntotv)

      do itotv=1,ntotv           ! Diagonal principal ensamblada
        lpamc(itotv)=itotv
      enddo

      nrigi=ntotv+nelem*((nevab*(nevab-1))/2) ! Posiciones de memoria en total

      return
      end

c___________________________________________________________________CAL$Proc_Rig
c     Procesamiento de la matriz de rigidez del elemento IELEM
c     K_NIV = 0 -> Calculo de matrices D, B, S=D.B
c           = 1 -> Calculo de matrices D, B, S=D.B, K=Bt.S
c           = 2 -> Calculo de matrices D, B, S=D.B, K=Bt.S y Ensamblaje CP
c           = 3 -> Calculo de matrices D, B, S=D.B, K=Bt.S y Almacenamiento EM

      Subroutine CAL$Proc_Rig(k_niv,ntype,ielem,nelem,nnode,
     .                        npoin,ndime,ndofn,nstre,nmats,
     .                        nprop,nevab,ntotv,nrigi,lnods,
     .                        coord,elcod,matno,props,dmatx,
     .                        bmatx,smatx,estif,lpamc,lpntn,
     .                        rgmat)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension lnods(nnode,nelem),matno(nelem),props(nprop,nmats)

c.....Seleccion de coordenadas nodales del elemento
      call SRV$Sel_Subvector(nnode,lnods(1,ielem),ndime,npoin,coord,
     .     elcod)

c.....Selecion del tipo de material
      imats=matno(ielem)

c.....Generacion de matriz de rigidez del elemento
      call CAL$Gen_K(ntype,nnode,ndime,ndofn,nstre,
     .               nprop,nevab,elcod,props(1,imats),
     .               dmatx,bmatx,smatx)

      if (k_niv.eq.0) return

c.....Calculo de matriz de rigidez del elemento
      dvolu=1.d+00
      call CAL$IN_K(nstre,nevab,1,dvolu,bmatx,smatx,estif)

      if (k_niv.eq.1) return

c.....Ensamblaje de la matriz de rigidez de elemento en la general
      if (k_niv.eq.2) then
        call CAL$Ensamb_K(nnode,ndofn,npoin,ntotv,nrigi,lnods(1,ielem),
     .       estif,lpamc,lpntn,rgmat)
      else if (k_niv.eq.3) then
        call CAL$Ensamb_Kd(nnode,ndofn,npoin,ntotv,nrigi,lnods(1,ielem),
     .       estif,lpamc,rgmat)
        call CAL$Almacenar_EM(nevab,ntotv,nrigi,ielem,estif,rgmat)
      endif

      return
      end

c__________________________________________________________________CAL$Proc_FPun
c     Procesamiento de la fuente puntual ILOAD de valor PFORC(ndofn)
c     aplicada en el nodo IPOIN, teniendo en cuenta los coeficientes de
c     seguridad FLOAD(ncase,iload)

      Subroutine CAL$Proc_FPun(ipoin,iload,ndofn,npoin,ntotv,
     .           nload,ncase,pforc,fload,lpntn,gload)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension pforc(ndofn),fload(ncase,nload),gload(ntotv,ncase)

      cero=0.d+00

c.....ensamblaje
      do icase=1,ncase
        if (fload(icase,iload).ne.cero) then
          call CAL$Ensamb_F(    1,ndofn,npoin,ipoin,pforc,
     .                      fload(icase,iload),lpntn,gload(1,icase))
        endif
      enddo

      return
      end

c__________________________________________________________________CAL$Proc_Tens
c     Procesamiento de intensidades en el elemento IELEM

      Subroutine CAL$Proc_Tens(ielem,nelem,ncase,nnode,ndofn,npoin,
     .           nstre,nevab,lnods,lpntn,xdisp,eldis,smatx,strsg)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension lnods(nnode,nelem),xdisp(ndofn,npoin,ncase),
     .          strsg(nstre,ncase)

      do icase=1,ncase

c.......selecciona desplazamientos del elemento
        call SRV$Sel_Subvector_2(nnode,lnods(1,ielem),lpntn,ndofn,npoin,
     .                           xdisp(1,1,icase),eldis)

c.......calculo de tensiones
        call CAL$Cal_Tens(nstre,nevab,
     .                    smatx,eldis,strsg(1,icase))

      enddo

      return
      end


c===============================================================================
c=============================== SEGUNDO NIVEL =================================
c===============================================================================

c______________________________________________________________________CAL$Gen_K
c     Generacion de la matriz de rigidez de elemento
   
      Subroutine CAL$Gen_K(ntype,nnode,ndime,ndofn,nstre,
     .                     nprop,nevab,elcod,props,
     .                     dmatx,bmatx,smatx)


      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension props(nprop)

c.....calculo de la matriz D
      call CAL$Cal_D(ntype,
     .               nprop,props,nnode,ndime,elcod,
     .               nstre,dmatx)

c.....generacion de la matriz B
      call CAL$Cal_B(ntype,
     .               nnode,ndime,elcod,
     .               nevab,nstre,bmatx)

c.....calculo de la matriz S=D.B
      call CAL$Cal_S(nstre,nevab,dmatx,bmatx,smatx)

      return
      end

c===============================================================================
c=========================== INTEGRACION NUMERICA ==============================
c===============================================================================

c_______________________________________________________________________CAL$IN_K
c                                                                      T
c     Integracion numerica de la matriz de rigidez del elemento: K  = B  . S
c                                                                ~e   ~    ~
c     (se calcula solo la triangular inferior, por simetria)

      Subroutine CAL$IN_K(nstre,nevab,ngaue,dvolu,bmatx,smatx,estif)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension dvolu(ngaue),bmatx(nstre,nevab,ngaue),
     .          smatx(nstre,nevab,ngaue),estif(nevab,nevab)

      cero=0.d+00

      do ievab=1,nevab
        do jevab=ievab,nevab
          estif(jevab,ievab)=cero
          do igaue=1,ngaue
            p=cero
            do istre=1,nstre
              p=p+bmatx(istre,ievab,igaue)*smatx(istre,jevab,igaue)
            enddo
            estif(jevab,ievab)=estif(jevab,ievab)+dvolu(igaue)*p
          enddo
        enddo
      enddo

      return
      end

c===============================================================================
c=========================== FORMULACION MATRICIAL =============================
c===============================================================================

c______________________________________________________________________CAL$Cal_D
c     Calculo de la matriz D
c                          ~

      Subroutine CAL$Cal_D(ntype,
     .                     nprop,props,nnode,ndime,elcod,
     .                     nstre,dmatx)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension props(nprop),elcod(ndime,nnode),dmatx(nstre,nstre)

      cero=0.d+00

      if (ntype.eq.2) then
        ea=props(1)
        dx=elcod(1,2)-elcod(1,1)
        dy=elcod(2,2)-elcod(2,1)
        rl=dsqrt(dx*dx+dy*dy)
        dmatx(1,1)=ea/rl
      else if (ntype.eq.3) then
        ea=props(1)
        dx=elcod(1,2)-elcod(1,1)
        dy=elcod(2,2)-elcod(2,1)
        dz=elcod(3,2)-elcod(3,1)
        rl=dsqrt(dx*dx+dy*dy+dz*dz)
        dmatx(1,1)=ea/rl
      else if (ntype.eq.4) then
        ea=props(1)
        ei=props(2)
        dx=elcod(1,2)-elcod(1,1)
        dy=elcod(2,2)-elcod(2,1)
        rl=dsqrt(dx*dx+dy*dy)
        dmatx(1,1)=ea/rl
        dmatx(2,1)=cero
        dmatx(3,1)=cero
        dmatx(1,2)=cero
        dmatx(2,2)=4.d+00*ei/rl
        dmatx(3,2)=2.d+00*ei/rl
        dmatx(1,3)=cero
        dmatx(2,3)=dmatx(3,2)
        dmatx(3,3)=dmatx(2,2)
      else
        call SRV$Error('CAL$CAL$Cal_D: Opcion no implementada (NTYPE)')
      endif

      return
      end

c______________________________________________________________________CAL$Cal_B
c     Calculo de la matriz B
c                          ~

      Subroutine CAL$Cal_B(ntype,
     .                     nnode,ndime,elcod,
     .                     nevab,nstre,bmatx)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension elcod(ndime,nnode),bmatx(nstre,nevab)

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

      if (ntype.eq.2) then
        dx=elcod(1,2)-elcod(1,1)
        dy=elcod(2,2)-elcod(2,1)
        rl=dsqrt(dx*dx+dy*dy)
        bmatx(1,1)=-dx/rl
        bmatx(1,2)=-dy/rl
        bmatx(1,3)=-bmatx(1,1)
        bmatx(1,4)=-bmatx(1,2)
      else if (ntype.eq.3) then
        dx=elcod(1,2)-elcod(1,1)
        dy=elcod(2,2)-elcod(2,1)
        dz=elcod(3,2)-elcod(3,1)
        rl=dsqrt(dx*dx+dy*dy+dz*dz)
        bmatx(1,1)=-dx/rl
        bmatx(1,2)=-dy/rl
        bmatx(1,3)=-dz/rl
        bmatx(1,4)=-bmatx(1,1)
        bmatx(1,5)=-bmatx(1,2)
        bmatx(1,6)=-bmatx(1,3)
      else if (ntype.eq.4) then
        dx=elcod(1,2)-elcod(1,1)
        dy=elcod(2,2)-elcod(2,1)
        rl=dsqrt(dx*dx+dy*dy)
        bmatx(1,1)=-dx/rl
        bmatx(2,1)=-dy/rl/rl
        bmatx(3,1)= bmatx(2,1)
        bmatx(1,2)=-dy/rl
        bmatx(2,2)= dx/rl/rl
        bmatx(3,2)= bmatx(2,2)
        bmatx(1,3)= cero
        bmatx(2,3)= uno
        bmatx(3,3)= cero
        bmatx(1,4)=-bmatx(1,1)
        bmatx(2,4)=-bmatx(2,1)
        bmatx(3,4)=-bmatx(3,1)
        bmatx(1,5)=-bmatx(1,2)
        bmatx(2,5)=-bmatx(2,2)
        bmatx(3,5)=-bmatx(3,2)
        bmatx(1,6)= cero
        bmatx(2,6)= cero
        bmatx(3,6)= uno
      else
        call SRV$Error('CAL$CAL$Cal_B: Opcion no implementada (NTYPE)')
      endif

      return
      end

c______________________________________________________________________CAL$Cal_S
c     Calculo de la matriz S = D . B
c                          ~   ~   ~

      Subroutine CAL$Cal_S(nstre,nevab,dmatx,bmatx,smatx)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

c.....Producto: SMATX(nstre,nevab)=DMATX(nstre,nstre)*BMATX(nstre,nevab)
      call SRV$Mult_M_M(nstre,nstre,nevab,dmatx,bmatx,smatx)

      return
      end

c___________________________________________________________________CAL$Cal_Tens
c                                 _     _
c     Calculo de las intensidades s = S a
c                                     ~

      Subroutine CAL$Cal_Tens(nstre,nevab,smatx,eldis,strsg)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension smatx(nstre,nevab),strsg(nstre)

c.....Producto: STRSG(nstre)=SMATX(nstre,nevab)*ELDIS(nevab)
      call SRV$Mult_M_M(nstre,nevab,1,smatx,eldis,strsg)

      return
      end


c===============================================================================
c============================= GESTION DE MEMORIA ==============================
c===============================================================================

c___________________________________________________________________CAL$Ensamb_K
c     Ensamblaje de la matriz de rigidez de elemento (ESTIF)
c     en la matriz de rigidez global (RGMAT) en perfil (LPAMC)
c     considerando la renumeracion nodal (LPNTN)

      Subroutine CAL$Ensamb_K(nnode,ndofn,npoin,ntotv,nrigi,
     .                        lnods,estif,lpamc,lpntn,rgmat)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)
   
      dimension lnods(nnode),estif(ndofn,nnode,ndofn,nnode),
     .          lpamc(ntotv),lpntn(npoin),rgmat(nrigi)

      do inode=1,nnode
        ipoir=lpntn(lnods(inode))
        do idofn=1,ndofn
          itotv=(ipoir-1)*ndofn+idofn

          jnode=inode
            jpoir=ipoir
            do jdofn=idofn,ndofn
              jtotv=(jpoir-1)*ndofn+jdofn
                lpjj=iabs(lpamc(jtotv))
                lpj0=lpjj-jtotv
                lp=lpj0+itotv
              rgmat(lp)=rgmat(lp)+estif(jdofn,jnode,idofn,inode)
            enddo

          do jnode=inode+1,nnode
            jpoir=lpntn(lnods(jnode))
            do jdofn=1,ndofn
              jtotv=(jpoir-1)*ndofn+jdofn
              if (jtotv.le.itotv) then
                lpii=iabs(lpamc(itotv))
                lpi0=lpii-itotv
                lp=lpi0+jtotv
              else
                lpjj=iabs(lpamc(jtotv))
                lpj0=lpjj-jtotv
                lp=lpj0+itotv
              endif
              rgmat(lp)=rgmat(lp)+estif(jdofn,jnode,idofn,inode)
            enddo
          enddo

        enddo
      enddo

      return
      end

c__________________________________________________________________CAL$Ensamb_Kd
c     Ensamblaje de diagonal de la matriz de rigidez de elemento (ESTIF)
c     en la diagonal de la matriz de rigidez global, que se almacena en
c     los primeros componentes del vector RGMAT

      Subroutine CAL$Ensamb_Kd(nnode,ndofn,npoin,ntotv,nrigi,
     .                        lnods,estif,lpamc,rgmat)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)
   
      dimension lnods(nnode),estif(ndofn,nnode,ndofn,nnode),
     .          lpamc(ntotv),rgmat(nrigi)

      do inode=1,nnode
        ipoin=lnods(inode)
        do idofn=1,ndofn
          itotv=(ipoin-1)*ndofn+idofn
          lpii=iabs(lpamc(itotv))
          rgmat(lpii)=rgmat(lpii)+estif(idofn,inode,idofn,inode)
        enddo
      enddo

      return
      end

c_______________________________________________________________CAL$Almacenar_EM
c     Ensamblaje de la diagonal de la matriz de rigidez de elemento (ESTIF)
c     en la matriz de rigidez global (RGMAT), y almacenamiento secuencial
c     de la parte restante (triangular inferior) de la matriz de rigidez
c     de elemento

      Subroutine CAL$Almacenar_EM(nevab,ntotv,nrigi,ielem,estif,rgmat)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)
   
      dimension estif(nevab,nevab),rgmat(nrigi)

      lp0=ntotv+(ielem-1)*((nevab*(nevab-1))/2)
      do jevab=2,nevab
        do ievab=1,jevab-1
          lpij=lp0+((jevab-1)*(jevab-2))/2+ievab
          rgmat(lpij)=estif(jevab,ievab)
        enddo
      enddo

      return
      end


c___________________________________________________________________CAL$Ensamb_F
c     Ensamblaje del vector de fuerzas de elemento (ELOAD)
c     en el vector de fuerzas global (GLOAD)
c     considerando la renumeracion nodal (LPNTN)
c     y teniendo en cuenta el coeficiente de seguridad (FLOAD)

      Subroutine CAL$Ensamb_F(nnode,ndofn,npoin,
     .                        lnods,eload,fload,lpntn,gload)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)
   
      dimension lnods(nnode),eload(ndofn,nnode),lpntn(npoin),
     .          gload(ndofn,npoin)

      do inode=1,nnode
        ipoir=lpntn(lnods(inode))
        do idofn=1,ndofn
          gload(idofn,ipoir)=gload(idofn,ipoir)+
     .                       eload(idofn,inode)*fload
        enddo
      enddo

      return
      end

C_________________________________________________________________CAL$Imp_Vincul
c     Formacion del vector de grados de libertad prescritos (NGFIX)
c     a partir de las condiciones de vinculacion (NOFIX)
c     considerando la renumeracion nodal (LPNTN)

      Subroutine CAL$Imp_Vincul(nvfix,ndofn,npoin,nofix,ngfix,lpntn)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension nofix(2,nvfix),ngfix(nvfix),lpntn(npoin)

      do ivfix=1,nvfix
        ipoin=nofix(1,ivfix)
        idofn=nofix(2,ivfix)
        ipoir=lpntn(ipoin)
        itotv=(ipoir-1)*ndofn+idofn
        ngfix(ivfix)=itotv
      enddo

      return
      end


c===============================================================================
c======================== Subrutinas Generales de Servicio =====================
c===============================================================================

c_________________________________________________________________SRV$Ini_Vector
c     Inicializacion de un vector a cero

      Subroutine SRV$Ini_Vector(n,v)

      implicit integer*4 (i-n)
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_________________________________________________________________SRV$Cop_Vector
c     Copia un vector en otro

      Subroutine SRV$Cop_Vector(n,v,w)

      implicit integer*4 (i-n)
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______________________________________________________________SRV$Sel_Subvector
c     Seleccion de un subvector mediante un redireccionamiento

      Subroutine SRV$Sel_Subvector(nn,ln,nd,np,co,el)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension ln(nn),co(nd,np),el(nd,nn)

      do in=1,nn
        ip=ln(in)
        do id=1,nd
          el(id,in)=co(id,ip)
        enddo
      enddo

      return
      end

c____________________________________________________________SRV$Sel_Subvector_2
c     Seleccion de un subvector mediante dos redireccionamientos consecutivos

      Subroutine SRV$Sel_Subvector_2(nn,ln,ln2,nd,np,co,el)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension ln(nn),ln2(np),co(nd,np),el(nd,nn)

      do in=1,nn
        ip=ln2(ln(in))
        do id=1,nd
          el(id,in)=co(id,ip)
        enddo
      enddo

      return
      end


c___________________________________________________________________SRV$Mult_M_M
c     Producto Matricial C         = A         . B
c                        ~ (n x m)   ~ (n x l)   ~ (l x m)

      Subroutine SRV$Mult_M_M(n,l,m,a,b,c)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension a(n,l),b(l,m),c(n,m)

      cero=0.d+00

      do j=1,m
        do i=1,n
          c(i,j)=cero
          do k=1,l
            c(i,j)=c(i,j)+a(i,k)*b(k,j)
          enddo  
        enddo
      enddo

      return
      end

c__________________________________________________________________SRV$Mult_M_Mt
c                                                  T
c     Producto Matricial C         = A         . B
c                        ~ (n x m)   ~ (n x l)   ~ (m x l)

      Subroutine SRV$Mult_M_Mt(n,l,m,a,b,c)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension a(n,l),b(m,l),c(n,m)

      cero=0.d+00

      do j=1,m
        do i=1,n
          c(i,j)=cero
          do k=1,l
            c(i,j)=c(i,j)+a(i,k)*b(j,k)
          enddo  
        enddo
      enddo

      return
      end

c__________________________________________________________________SRV$Mult_Mt_M
c                                      T
c     Producto Matricial C         = A         . B
c                        ~ (n x m)   ~ (l x n)   ~ (l x m)

      Subroutine SRV$Mult_Mt_M(n,l,m,a,b,c)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      dimension a(l,n),b(l,m),c(n,m)

      cero=0.d+00

      do j=1,m
        do i=1,n
          c(i,j)=cero
          do k=1,l
            c(i,j)=c(i,j)+a(k,i)*b(k,j)
          enddo  
        enddo
      enddo

      return
      end

c______________________________________________________________________SRV$Error
c     Impresion de mensaje de error y detencion del programa

      Subroutine SRV$Error(mensaje)

      implicit integer*4 (i-n)
d     implicit real*8    (a-h,o-z)

      character*(*) mensaje

      write(6,100) mensaje
 100  format(/' ERROR > ',a/)

      stop
      end
