PROGRAM BUCKLE
C======================================================================
C AN FEM SOLVER FOR BUCKLING PROBLEM WITH MOMENT AT BOTH ENDS OF BEAM
C        EQUATION: D2U/DXDX + ALPHA*U=0; ALPHA = P/(EI)
C P=APPLIED FORCE AT ENDS OF BEAM, E=YOUNG MODULUS, I=2ND MOMENT OF
C INERTIA. RL=LENGTH OF ELEMENT, IBTYPE(1)=BOUNDARY CONDITION AT THE
C LEFT END OF BEAM, IBTYPE(2)=BC AT RIGHT END. IBTYPE(I)=1 FOR FIXED 
C DISPLACEMENT, IBTYPE(I)=2 FOR PRESCRIBED SLOPE AT THE END.
C MAY 16, 1994   EIJI FUKUMORI
C======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=2, MXE=10, MXN=MXE+1 )
      DIMENSION NODEX(MXE,ND),EI(MXE),X(MXN),A(MXN,MXN),RHS(MXN),
     * IBTYPE(2), BV(2), STIFF(ND,ND)
C======================================================================
C (1) READING OF DATA
      CALL INPUT ( MXE,MXN,ND,P,NE,NNODE,NODEX,EI,X,IBTYPE,BV )
C======================================================================
C (2) CONSTRUCTION OF FEM-MATRIX EQUATION
      CALL MATRIX ( MXE,MXN,ND,P,NE,NNODE,STIFF,NODEX,EI,X,A,RHS )
C======================================================================
C (3) IMPLEMENTATION OF BOUNDARY CONDITION
      CALL FORM ( MXN, NNODE, A, RHS, IBTYPE, BV )
C======================================================================
C (4) SOLVE FOR UNKNOWN VARIABLES
      CALL SYSTEM ( MXN, NNODE, A, RHS )
C======================================================================
C (5) PRINTING RESULTS
      DO I = 1 , NNODE
      WRITE(*,*)' NODAL # =',I, '     DISPLACEMENT =', RHS(I)
      END DO
      STOP' NORMAL TERMINATION'
      END
C
C
      SUBROUTINE INPUT ( MXE,MXN,ND,P,NE,NNODE,NODEX,EI,X,IBTYPE,BV )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),EI(MXE),X(MXN),IBTYPE(2),BV(2)
      OPEN ( 1,FILE='BUCKLE.DAT', STATUS='OLD')
      READ(1,*) P
      READ(1,*) NE
      DO I = 1 , NE
      READ(1,*) IEL, (NODEX(IEL,J),J=1,ND), EI(IEL)
      END DO
      NNODE = NE + 1
      DO I = 1 , NNODE
      READ(1,*) NODE, X(NODE)
      END DO
      READ(1,*) IBTYPE(1), BV(1)
      READ(1,*) IBTYPE(2), BV(2)
      CLOSE (1)
      RETURN
      END
C
C
      SUBROUTINE MATRIX (MXE,MXN,ND,P,NE,NNODE,STIFF,NODEX,EI,X,A,RHS)
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),EI(MXE),X(MXN),A(MXN,MXN),RHS(MXN),
     * STIFF(ND,ND)
      DO I = 1 , NNODE
      DO J = 1 , NNODE
      A(I,J) = 0.
      END DO
      RHS(I) = 0.
      END DO
      DO IEL = 1 , NE
      I = NODEX(IEL,1)
      J = NODEX(IEL,2)
      RL = X(J) - X(I)
      ALPHA = P / EI(I)
      STIFF(1,1) = -1./RL + ALPHA*RL/3.
      STIFF(1,2) =  1./RL + ALPHA*RL/6.
      STIFF(2,1) =  1./RL + ALPHA*RL/6.
      STIFF(2,2) = -1./RL + ALPHA*RL/3.
      A(I,I) = A(I,I) + STIFF(1,1)
      A(I,J) = A(I,J) + STIFF(1,2)
      A(J,I) = A(J,I) + STIFF(2,1)
      A(J,J) = A(J,J) + STIFF(2,2)
      END DO
      RETURN
      END
C
C
      SUBROUTINE FORM ( MXN, NNODE, A, RHS, IBTYPE, BV )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION IBTYPE(2),BV(2),A(MXN,MXN), RHS(MXN)
      IF ( IBTYPE(1) .EQ. 1 ) THEN
      A(1,1) = 1.
      A(1,2) = 0.
      RHS(1) = BV(1)
      RHS(2) = RHS(2) - BV(1)*A(2,1)
      A(2,1) = 0.
      ELSE
      RHS(1) = RHS(1) - BV(1)
      END IF
      IF ( IBTYPE(2) .EQ. 1 ) THEN
      A(NNODE,NNODE) = 1.
      A(NNODE,NNODE-1) = 0.
      RHS(NNODE) = BV(2)
      RHS(NNODE-1) = RHS(NNODE-1) - BV(2)*A(NNODE-1,NNODE)
      A(NNODE-1,NNODE) = 0.
      ELSE
      RHS(NNODE) = RHS(NNODE) - BV(2)
      END IF
      RETURN
      END
C
C
      SUBROUTINE SYSTEM ( MXN, N , A , C )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION A (MXN,MXN) , C (MXN)
      N1 = N - 1
      DO K = 1, N1
      L = K + 1
       DO I = L , N
        P = A (I,K) / A (K,K)
        IF ( P .NE. 0. ) THEN
         DO J = L , N
          A (I,J) = A (I,J) - P * A ( K , J )
         END DO
         C ( I ) = C ( I) - P * C ( K )
        END IF
       END DO
      END DO
C------- BACK SUBSTITUTION ---
      C (N) = C (N) / A (N,N)
      DO K = 1 , N1
       I = N - K
       L = I + 1
       P = C ( I )
       DO J = L , N
        P = P - C (J) * A (I,J)
       END DO
       C ( I ) = P / A (I,I)
      END DO
      RETURN
      END