c_______________________________________________________________________________
c
c     Calculo de la carga critica de pandeo de una viga biapoyada
c
c     Sea una viga biapoyada de longitud L y rigidez a flexión EI.
c     La viga está sometida a una carga axial de valor P.
c     Sea v(chi) la deflexion en cada punto x (entre 0 y L), donde chi = x / L
c     es la variable espacial adimensionalizada.
c
c     El pandeo se produce cuando existen configuraciones de equilibrio en las
c     que las deflexiones v no son identicamente nulas, aunque no se
c     hayan aplicado cargas perpendiculares al eje de la viga.
c
c     Se llaman cargas críticas a los valores de la carga P que producen pandeo.
c
c     Para determinar las posibles cargas críticas se realiza un analisis de
c     segundo orden, esto es: se supone que los desplazamientos y sus gradientes
c     son pequeños, pero se plantea el equilibrio sobre la configuracion
c     deformada.
c
c     La ecuacion diferencial que describe el fenomeno es
c 
c                           2                                2     2
c          v"(chi)= - lambda  v(chi),   0 < chi < 1,   lambda  = PL /EI > 0
c
c          v(0) = 0, v(1) = 0
c
c     Esta ecuación solo admite solucion para determinados valores de lambda.
c
c     Para cada uno de ellos, el correspondiente valor de P es una posible carga
c     critica y la solucion asociada v(chi) describe la correspondiente forma de
c     pandeo (esto es: la forma que adquiere el eje de la viga al deformarse). 
c
c     Para resolver el problema anterior se discretiza el intervalo [0,1] en
c     n subintervalos iguales
c
c          [0, 1] = U [chi   , chi ],     chi = x / L = i / n,     i=1,___,n,
c                         i-1     i          i   i
c
c     lo que permite sustituir la derivada segunda en los nodos interiores
c     (i=1,___,n-1) por las aproximaciones en diferencias de segundo orden
c
c                                                              2      2
c          v"(chi ) = (- v(chi   ) + 2 v(chi ) - v(chi   )) / h  + O(h ),   
c                i            i-1           i         i+1
c
c     donde
c
c          h = chi  - chi    = 1 / n,     i = 1,___, n.
c                 i      i-1 
c
c     Despreciando los errores locales de truncamiento (de segundo orden),
c     teniendo en cuenta las condiciones de contorno e introduciendo las
c     aproximaciones
c
c          ^
c          v  ~ v(chi ),     i = 0,___,n,
c           i        i
c
c     se obtiene el problema de valores y vectores propios
c
c
c                                  |  2 -1          |
c                                  | -1  2 -1       |
c            _      _              |   .  .  .      |              2    2
c          A u = mu u,   con   A = |     .  .  .    |,    mu = (1/n ) PL /EI > 0
c          ~                   ~   |       .  .  .  |
c                                  |       -1  2 -1 |
c                                  |          -1  2 |
c
c                                    _
c     cuya solucion es el autovector u y su autovalor asociado mu, donde
c
c              | u    |
c              |  1   |
c          _   | .    |         ^           2       2         2  2
c          u = | .    |,   u  = v ,   mu = h  lambda  = (1/n ) PL /EI > 0.
c              | .    |     i    i
c              | u    |
c              |  n-1 |
c
c     El problema tiene (n-1) soluciones. Los autovectores u son los modos de
c     pandeo (las componentes son los desplazamientos de los nodos). Los
c     Los correspondientes autovalores mu determinan las cargas criticas P
c     asociadas a cada uno de los modos.
c
c     Nos interesa la solucion correspondiente a la carga critica mas baja.
c
c     Aplicando el método de MISES a la matriz inversa de A obtendremos el
c     valor de (1/mu), que se corresponde con el menor valor de mu.
c
c     El algoritmo se aplica de la siguiente forma:
c               _      _
c          Dado v ,   |v | = 1  
c                0      0
c                      _     -1 _           _ t
c          Calculamos: w  = A   v ,   rho = v   w ,
c                       0        0       0   0   0
c
c          e iteramos hasta convergencia con
c                      _     _     _      _       -1 _                _ t  _
c                      v   = w  / |w |,   w    = A   v   ,   rho    = v    w   .
c                       k+1   k     k      k+1   ~    k+1       k+1    k+1  k+1
c
c     de forma que     
c          _          _                                     2     2
c          u =   lim  v ,   (1/mu) =   lim  rho ,   P = mu n  EI/L .
c               k->oo  k              k->oo    k
c
c                    _    -1 _
c     Los productos  w = A   v  se calculan resoviendo el sistema de ecuaciones
c                        ~
c            _   _
c          A w = v.
c          ~
c
c_______________________________________________________________________________

      implicit real*8 (a-h,o-z)
      parameter(DCERO=0.d+00,
     .          PI=3.1415926535897932d+00)

      logical cnvgc
      allocatable v(:),w(:),v_l(:),a(:)

c.....Lectura de datos
      write(6,'()')

    1 write(6,'(a,$)') 'n = '
      read(5,*) n
      if (n.lt.2) goto 1

c.....Solicitud de memoria
      allocate(v(n-1),STAT=ist)
      if (ist.ne.0) call ErrorFatal('No hay memoria para vector v',1)
      allocate(w(n-1),STAT=ist)
      if (ist.ne.0) call ErrorFatal('No hay memoria para vector w',2)
      allocate(v_l(n-1),STAT=ist)
      if (ist.ne.0) call ErrorFatal('No hay memoria para vector v_l',3)
      allocate(a(2*(n-1)),STAT=ist)
      if (ist.ne.0) call ErrorFatal('No hay memoria para matriz a',4)

c.....Vector inicial (aleatorio)
      call Generar(n-1,v)       ! generacion
      call Normalizar(n-1,v,v)  ! normalizacion

c.....Ensamblaje de la matriz A
      call EM(n-1,a)

c.....Factorizacion de Cholesky
      call FCh(n-1,a)           ! A = L Lt

c.....Iteracion 0
      call SSEL(n-1,a,v,w)      ! solucion de A w = v mediante sustituciones
      call PE(n-1,v,w,rho)      ! producto escalar rho = v · w
      k=0
      cnvgc=.FALSE.

c.....Siguientes iteraciones
      do while (.not.cnvgc)
        rho_l=rho                 ! rho_l = valor anterior de rho
        call CopiaV(n-1,v,v_l)    ! v_l   = valor anterior del vector v
        call Normalizar(n-1,w,v)  ! normalizacion
        call SSEL(n-1,a,v,w)      ! solucion de A w = v mediante sustituciones
        call PE(n-1,v,w,rho)      ! producto escalar rho = v · w
        k=k+1
        call Convergencia(n-1,v_l,v,rho_l,rho,cnvgc)
      enddo

c.....Impresion del resultado
      fmu=1.d+00/rho
      p=fmu*dble(n)**2           ! P = mu n^2 (EI/L^2)
      err=1.d+00-p/PI**2         ! solucion exacta: P = pi^2 (EI/L^2)
      write(6,100) k,p,err*100.d+00
  100 format(/,'Numero de iteraciones   = ',i9,
     .       /,'Carga critica de pandeo = ',f20.12,' EI/L^2',
     .       /,'Error relativo          = ',f20.12,' %'/)

      read(5,'()')
      call exit(0)

      end


c________________________________________________________________________Generar
c
c     Generacion de un vector de numeros aleatorios
c_______________________________________________________________________________

      Subroutine Generar(n,v)
      implicit real*8 (a-h,o-z)

      dimension v(n)

c.....Obtencion de la semilla
    1 write(6,'(a,$)') 'Semilla (No. entero > 0) ........... iseed = '
      read(5,*) iseed
      if (iseed.le.0) goto 1

c.....Generacion de vector aleatorio
      call srand(iseed)
      do i=1,n
        v(i)=dble(rand())
      enddo

      return
      end


c_____________________________________________________________________Normalizar
c
c     Normalizacion de un vector     ( v  =  w / |w| )
c_______________________________________________________________________________

      Subroutine Normalizar(n,w,v)
      implicit real*8 (a-h,o-z)

      dimension w(n),v(n)

c.....Calculo del modulo
      call PE(n,w,w,wm)
      wm=sqrt(wm)

c.....Normalizacion
      do i=1,n
        v(i)=w(i)/wm
      enddo

      return
      end


c_____________________________________________________________________________EM
c
c     A es tridiagonal simétrica. Se almacena su parte inferior en banda.
c_______________________________________________________________________________

      Subroutine EM(n,a)
      implicit real*8 (a-h,o-z)
      dimension a(0:1,n)

        a(1,1)= 2.d+00
      do i=2,n
        a(0,i)=-1.d+00
        a(1,i)= 2.d+00
      enddo

      return
      end


c____________________________________________________________________________FCh
c
c     Factorizacion de Cholesky T = L Lt
c
c     T es tridiagonal simétrica. Se almacena su parte inferior en banda.
c     L se almacena sobre T.
c_______________________________________________________________________________

      Subroutine FCh(n,t)
      implicit real*8 (a-h,o-z)
      dimension t(0:1,n)

      t(1,1)=sqrt(t(1,1))
      do i=2,n
        t(0,i)=t(0,i)/t(1,i-1)
        t(1,i)=sqrt(t(1,i)-t(0,i)*t(0,i))
      enddo

      return
      end
  

c___________________________________________________________________________SSEL
c
c     Solucion de T w = v, con T = L Lt
c
c     T es tridiagonal simétrica. Se almacena su parte inferior en banda.
c     La factorizacion fue resuelta previamente. L ocupa el lugar de T.
c_______________________________________________________________________________

      Subroutine SSEL(n,t,v,w)
      implicit real*8 (a-h,o-z)
      dimension t(0:1,n),v(n),w(n)

c.....Sustitucion hacia delante
      w(1)=v(1)/t(1,1)
      do i=2,n
        w(i)=(v(i)-t(0,i)*w(i-1))/t(1,i)
      enddo

c.....Sustitucion hacia atras
      w(n)=w(n)/t(1,n)
      do i=n-1,1,-1
        w(i)=(w(i)-t(0,i+1)*w(i+1))/t(1,i)
      enddo

      return
      end


c_____________________________________________________________________________PE
c
c     Calculo del producto escalar p = v · w
c_______________________________________________________________________________

      Subroutine PE(n,v,w,p)
      implicit real*8 (a-h,o-z)
      parameter(DCERO=0.d+00)

      dimension v(n),w(n)

      p=DCERO
      do i=1,n
        p=p+v(i)*w(i)
      enddo
 
      return
      end


c_________________________________________________________________________CopiaV
c
c     Copia de un vector v2 = v1
c_______________________________________________________________________________

      Subroutine CopiaV(n,v1,v2)
      implicit real*8 (a-h,o-z)

      dimension v1(n),v2(n)

      do i=1,n
        v2(i)=v1(i)
      enddo

      return
      end


c___________________________________________________________________Convergencia
c
c     Chequeo de convergencia en el autovector y en el autovalor
c_______________________________________________________________________________

      Subroutine Convergencia(n,v_l,v,rho_l,rho,cnvgc)
      implicit real*8 (a-h,o-z)
      parameter(DCERO=0.d+00)

      logical cnvgc
      dimension v_l(n),v(n)

      epsl=1.d-20             !error absoluto maximo admisible
      rltv=1.d-12             !error relativo maximo admisible

c.....Convergencia en el autovalor
      dr  =abs(rho-rho_l)
      drmx=abs(rho)*rltv
      if (epsl.gt.drmx) drmx=epsl
      cnvgc=(dr.le.drmx)
      if (.not.cnvgc) return

c.....Convergencia en el autovector (norma del maximo)
      dnorm=DCERO
      vnorm=DCERO
      do i=1,n
        di=abs(v(i)-v_l(i))
        if (di.gt.dnorm) dnorm=di   ! norma del maximo de (v - v_l)
        vi=abs(v(i))
        if (vi.gt.vnorm) vnorm=vi   ! norma del maximo de v
      enddo
      dvmx=vnorm*rltv
      if (epsl.gt.dvmx) dvmx=epsl
      cnvgc=(dnorm.le.dvmx)

      return
      end


c_____________________________________________________________________ErrorFatal
c
c     Impresion de mensaje de error y detencion del programa
c_______________________________________________________________________________


      Subroutine ErrorFatal(mensaje,ierr)
      character*(*) mensaje

      write(6,100) mensaje
  100 format(' Error: ',a)

      call exit(ierr)
      end
 