program f90_EGV_BFRM
    !*************************************************************
    !2D Frame Buckling Analysis
    !*************************************************************
    implicit none
    integer::i,j,k,l,m,kk,ne,n,ia,ja
    character::strcom*50                !Comment
    integer::NODT                       !Number of nodes
    integer::NELT                       !Number of elements
    integer::MATEL                      !Number of material sets
    integer::KOX                        !Number of fixed nodes in x-direction
    integer::KOY                        !Number of fixed nodes in y-direction
    integer::KOZ                        !Number of fixed nodes in rotation
    integer::NF                         !Number of loaded nodes
    integer,allocatable::kakom(:,:)     !Element connectivity
    integer,allocatable::matno(:)       !Material set number
    real(8),allocatable::wm(:,:)        !work for material properties
    real(8),allocatable::Em(:)          !Elastic modulus of element
    real(8),allocatable::AA(:)          !Section area of element
    real(8),allocatable::AI(:)          !Moment of inertia of element
    real(8),allocatable::x(:)           !x-coordinate
    real(8),allocatable::y(:)           !y-coordinate
    integer,allocatable::nokx(:)        !Fixed node number in x-direction
    integer,allocatable::noky(:)        !Fixed node number in y-direction
    integer,allocatable::nokz(:)        !Fixed node number in rotation
    real(8),allocatable::f(:)           !External force vector
    integer,allocatable::mindex(:)      !Index for treatment of fixed nodes
    real(8),allocatable::ftvec(:)       !Total load vector
    real(8),allocatable::tsm(:,:)       !Total stiffness matrix
    real(8),allocatable::dis(:)         !Displacement
    real(8),allocatable::reac(:)        !work array
    real(8),allocatable::sresul(:,:)    !Stress resultants
    real(8),allocatable::sm(:,:)        !Element stiffnes matrix
    real(8),allocatable::hen(:,:)       !Transform matrix
    real(8),allocatable::fevec(:)       !Element load vector
    integer::nt                         !Number of degrees of freedom of all nodes (NODT x nhen)
    integer::nod=2                      !Number of nodes per element
    integer::nhen=3                     !Number of degrees of freedom per node
    integer::mm                         !Number of degree of freedom of reduced FE equation
    integer::ib                         !Band width of reduced FE equation
    real(8)::fact0=1.0D-30              !0 judgement
    character::fnameR*50,fnameW*50      !Input file name, output file name

    real(8),allocatable::smg(:,:)       !Element geometric stiffness matrix
    real(8),allocatable::tsg(:,:)       !Total geometric stiffness matrix
    real(8),allocatable::ev(:)          !Eigenvalue
    real(8),allocatable::vec(:,:)       !Eigenvector
    real(8),allocatable::work(:)        !work for eigenvector
    integer::negv                       !Number of outputted eigenvalues
    integer::kst                        !Index for positive eigenvalue

    character :: linebuf*1000,dummy*100
    character :: fmt1*200,fmt2*200,fmt3*200,fmt4*200
    real(4) :: t1,t2

    fmt1="(i5,2(',',e15.7),3(',',e15.7),3(',',i2))"                  !Node
    fmt2="(i5,2(',',i5),3(',',e15.7),',',i5)"                        !Element
    fmt3="(i5,2(',',e15.7),3(',',e15.7),3(',',e15.7),3(',',e15.7))"  !Displacement & reaction
    fmt4="(i5,3(',',e15.7),3(',',e15.7))"                            !Stress resultant

    call getarg(1,dummy)  !Number of outputted eigenvalues
    call getarg(2,fnameR) !Input file name
    call getarg(3,fnameW) !Output file name
    read(dummy,*) negv    !To convert character 'dummy' to integer 'negv'

    !To start the measurement of calculation time
    t1=DIFT()
    !***************************
    !data input
    !***************************
    open(11,file=fnameR,status='old')
    !/*--------------------------------------------------------------------------------------*/
    !/*Comment*/
    !/*--------------------------------------------------------------------------------------*/
    read(11,'(a)') strcom
    !------------------------------------------------------------------------------------------
    !Basic values for analysis
    !------------------------------------------------------------------------------------------
    read(11,*) NODT,NELT,MATEL,KOX,KOY,KOZ,NF
    !------------------------------------------------------------------------------------------
    !To declare array size
    !------------------------------------------------------------------------------------------
    nt=NODT*nhen
    allocate(kakom(1:NELT,1:nod))
    allocate(matno(1:NELT))
    allocate(wm(1:NELT,1:3))
    allocate(Em(1:NELT))
    allocate(AA(1:NELT))
    allocate(AI(1:NELT))
    allocate(x(1:NODT))
    allocate(y(1:NODT))
    if(1<=KOX)allocate(nokx(1:KOX))
    if(1<=KOY)allocate(noky(1:KOY))
    if(1<=KOZ)allocate(nokz(1:KOZ))
    allocate(f(1:nt))
    allocate(mindex(1:nt))
    allocate(ftvec(1:nt))
    allocate(tsm(1:nt,1:nt))
    allocate(dis(1:nt))
    allocate(reac(1:nt))
    allocate(sresul(1:NELT,1:nod*nhen))
    allocate(sm(1:nod*nhen,1:nod*nhen))
    allocate(hen(1:nod*nhen,1:nod*nhen))
    allocate(fevec(1:nod*nhen))

    allocate(smg(1:nod*nhen,1:nod*nhen))
    allocate(tsg(1:nt,1:nt))
    allocate(ev(1:nt))
    allocate(vec(1:nt,1:nt))
    allocate(work(1:nt))

    !------------------------------------------------------------------------------------------
    !Material properties (Em,AA,AI)
    !------------------------------------------------------------------------------------------
    do i=1,MATEL
        read(11,*) wm(i,1),wm(i,2),wm(i,3)
    end do
    !------------------------------------------------------------------------------------------
    !Element-nodes relationship, material set number (node-1,node-2,matno)
    !------------------------------------------------------------------------------------------
    do ne=1,NELT
        read(11,*) kakom(ne,1),kakom(ne,2),matno(ne)
        do i=1,MATEL
            if(i==matno(ne))then
                Em(ne)    =wm(i,1)
                AA(ne)   =wm(i,2)
                AI(ne)   =wm(i,3)
            end if
        end do
    end do
    deallocate(wm)
    !------------------------------------------------------------------------------------------
    !(x,y) coordinates of node
    !------------------------------------------------------------------------------------------
    do i=1,NODT
        read(11,*) x(i),y(i)
    end do
    !------------------------------------------------------------------------------------------
    !Fixed node number
    !------------------------------------------------------------------------------------------
    if(0<KOX)then
        do i=1,KOX
            read(11,*) nokx(i)
        end do
    end if
    if(0<KOY)then
        do i=1,KOY
            read(11,*) noky(i)
        end do
    end if
    if(0<KOZ)then
        do i=1,KOZ
            read(11,*) nokz(i)
        end do
    end if
    !------------------------------------------------------------------------------------------
    !Loaded node number, load in x-direction, load in y-direction, moment
    !------------------------------------------------------------------------------------------
    do i=1,nt
        f(i)=0.0D0
    end do
    if(0<NF)then
        do i=1,NF
            read(11,*) n,f(3*n-2),f(3*n-1),f(3*n)
        end do
    end if
    close(11)

    !*****************************************************
    !To print input data
    !*****************************************************
    open(12, file=fnameW, status='replace')
    write(12,'(a)') strcom
    write(12,'(a)') 'NODT,NELT,MATEL,KOX,KOY,KOZ,NF'
        write(linebuf,*) NODT,',',NELT,',',MATEL,',',KOX,',',KOY,',',KOZ,',',NF
        call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    write(12,'(a)') '*node characteristics'
    write(12,'(a)') 'node,x,y,fx,fy,fz,fix-x,fix-y,fix-z'
    do i=1,NODT
        k=0
        l=0
        m=0
        if(0<KOX)then
            do j=1,KOX
                if(i==nokx(j))k=1
            end do
        end if
        if(0<KOY)then
            do j=1,KOY
                if(i==noky(j))l=1
            end do
        end if
        if(0<KOZ)then
            do j=1,KOZ
                if(i==nokz(j))m=1
            end do
        end if
            write(linebuf,FMT=fmt1) i,x(i),y(i),f(3*i-2),f(3*i-1),f(3*i),k,l,m
            call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    write(12,'(a)') '*element characteristics'
    write(12,'(a)') 'element,node-1,node-2,E,A,I,matno'
    do ne=1,NELT
            write(linebuf,FMT=fmt2) ne,kakom(ne,1),kakom(ne,2),Em(ne),AA(ne),AI(ne),matno(ne)
            call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    close(12)
    !*****************************************************
    !Main calculation
    !*****************************************************
    !-----------------------------------------------------
    !Initialization of matrices and vectors
    !-----------------------------------------------------
    do i=1,nt
        ftvec(i)=0.0D0
        reac(i)=0.0D0
        do j=1,nt
            tsm(i,j)=0.0D0
        end do
    end do
    !-----------------------------------------------------
    !To assemble matrix
    !-----------------------------------------------------
    do ne=1,NELT
        !Stiffness matrix
        call SM_FRM(ne,kakom,x,y,Em,AA,AI,sm,hen,NELT,NODT)
        call ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT)
    end do
    !-----------------------------------------------------
    !To set load vector
    !-----------------------------------------------------
    do i=1,nt
        ftvec(i)=f(i)
    end do
    !-----------------------------------------------------
    !Treatment of fixed nodes
    !-----------------------------------------------------
    do i=1,nt
        mindex(i)=i
        reac(i)=ftvec(i)
    end do
    if(0<KOX)then
        do i=1,KOX
            n=3*nokx(i)-2
            mindex(n)=0
        end do
    end if
    if(0<KOY)then
        do i=1,KOY
            n=3*noky(i)-1
            mindex(n)=0
        end do
    end if
    if(0<KOZ)then
        do i=1,KOZ
            n=3*nokz(i)
            mindex(n)=0
        end do
    end if
    mm=0
    do i=1,nt
        if(mindex(i)>0)then
            mm=mm+1
            mindex(mm)=i
        end if
    end do
    do i=1,mm
        ia=mindex(i)
        reac(i)=reac(ia)
        do j=1,mm
            ja=mindex(j)
            tsm(i,j)=tsm(ia,ja)
        end do
    end do
    !-----------------------------------------------------
    !To memory band matrix
    !-----------------------------------------------------
    ib=IBAND(mm,tsm,fact0,nt)
    call BAND(mm,ib,tsm,nt)
    !-----------------------------------------------------
    !To verify diagonal element & to stop if necessary
    !-----------------------------------------------------
    do i=1,mm
        if(abs(tsm(i,1))<fact0)stop
    end do
    !-----------------------------------------------------
    !Cholesky method
    !-----------------------------------------------------
    call CHBAND10(mm,ib,tsm,nt)!Triangular matrix
    call CHBAND11(mm,ib,tsm,reac,nt)!Forward elimination & backward substitution

    !*****************************************************
    !Culculation of stress resultant
    !*****************************************************
    do i=1,nt
        dis(i)=0.0D0
    end do
    do i=1,mm
        ia=mindex(i)
        dis(ia)=reac(i)
    end do
    do ne=1,NELT
        call SM_FRM(ne,kakom,x,y,Em,AA,AI,sm,hen,NELT,NODT)
        call CALST_FRM(ne,kakom,dis,sm,hen,reac,sresul,NELT,NODT,nt)
    end do
    !*****************************************************
    !To print the results of small load model
    !*****************************************************
    open(12, file=fnameW, access='append')
    write(12,'(a)') '*displacements and forces'
    write(12,'(a)') 'node,coord-x,coord-y,dis-x,dis-y,dis-z,reac-x,reac-y,reac-z,ftvec-x,ftvec-y,ftvec-z'
    do i=1,NODT
        write(linebuf,FMT=fmt3) i,x(i),y(i),dis(3*i-2),dis(3*i-1),dis(3*i),&
            reac(3*i-2),reac(3*i-1),reac(3*i),ftvec(3*i-2),ftvec(3*i-1),ftvec(3*i)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    write(12,'(a)') '*stress resultants'
    write(12,'(a)') 'element,Ni,Si,Mi,Nj,Sj,Mj'
    do ne=1,NELT
        write(linebuf,FMT=fmt4) ne,sresul(ne,1),sresul(ne,2),sresul(ne,3),sresul(ne,4),sresul(ne,5),sresul(ne,6)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    close(12)
    !*****************************************************
    !Eigenvalue analysis
    !*****************************************************
    !-----------------------------------------------------
    !Initialization of matrices
    !-----------------------------------------------------
    do i=1,nt
        do j=1,nt
            tsm(i,j)=0.0D0
            tsg(i,j)=0.0D0
        end do
    end do
    !-----------------------------------------------------
    !To assemble matrices
    !-----------------------------------------------------
    do ne=1,NELT
        !Linear stiffness matrix
        call SM_FRM(ne,kakom,x,y,Em,AA,AI,sm,hen,NELT,NODT)
        call ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT)
        !Geometric stiffness matrix
        call GM_FRM(ne,kakom,x,y,Em,AA,AI,sresul,smg,hen,NELT,NODT)
        call ASMAT(ne,nod,nhen,kakom,smg,tsg,NELT,NODT)
    end do
    !-----------------------------------------------------
    !Treatment of fixed nodes
    !-----------------------------------------------------
    do i=1,nt
        mindex(i)=i
        reac(i)=ftvec(i)
    end do
    if(0<KOX)then
        do i=1,KOX
            n=3*nokx(i)-2
            mindex(n)=0
        end do
    end if
    if(0<KOY)then
        do i=1,KOY
            n=3*noky(i)-1
            mindex(n)=0
        end do
    end if
    if(0<KOZ)then
        do i=1,KOZ
            n=3*nokz(i)
            mindex(n)=0
        end do
    end if
    mm=0
    do i=1,nt
        if(mindex(i)>0)then
            mm=mm+1
            mindex(mm)=i
        end if
    end do
    do i=1,mm
        ia=mindex(i)
        do j=1,mm
            ja=mindex(j)
            tsm(i,j)=tsm(ia,ja)
            tsg(i,j)=tsg(ia,ja)
        end do
    end do
    !-----------------------------------------------------
    !Eigenvalue analysis
    !-----------------------------------------------------
    call GJACOBI(mm,tsm,tsg,ev,vec,nt)
    do k=1,mm
        do i=1,nt
            work(i)=0.0D0
        end do
        do i=1,mm
            ia=mindex(i)
            work(ia)=vec(i,k)
        end do
        do i=1,nt
            vec(i,k)=work(i)
        end do
    end do
    !*****************************************************
    !To print the results of eigenvalue analysis
    !*****************************************************
    open(12, file=fnameW, access='append')
    write(12,'(a)') '*eigen value & eigen vector'
    !Degree number
    kst=0
    do i=1,mm
        write(6,*) i,ev(i)
        if(0.0D0<ev(i))exit
    end do
    kst=i
    kk=kst+negv-1
    if(mm<kk)kk=mm
    linebuf='mode.'
    do i=kst,kk
        write(dummy,'(",",i0)') i
        linebuf=trim(linebuf)//dummy
    end do
    call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    !Eigenvalue
    linebuf='eigen value'
    do i=kst,kk
        write(dummy,'(",",e15.7)') ev(i)
        linebuf=trim(linebuf)//dummy
    end do
    call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    !Eigenvector
    linebuf='mode'
    do i=kst,kk
        write(dummy,'(",",i0)') i
        linebuf=trim(linebuf)//dummy
    end do
    call del_spaces(linebuf)
    write (12,'(a)') trim(linebuf)
    do i=1,NODT
        write(linebuf,'("dis-x_",i0)') i
        do k=kst,kk
            write(dummy,'(",",e15.7)') vec(3*i-2,k)
            linebuf=trim(linebuf)//dummy
        end do
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    do i=1,NODT
        write(linebuf,'("dis-y_",i0)') i
        do k=kst,kk
            write(dummy,'(",",e15.7)') vec(3*i-1,k)
            linebuf=trim(linebuf)//dummy
        end do
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    do i=1,NODT
        write(linebuf,'("dis-z_",i0)') i
        do k=kst,kk
            write(dummy,'(",",e15.7)') vec(3*i,k)
            linebuf=trim(linebuf)//dummy
        end do
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do


    !To finish the measurement of calculation time
    t2=DIFT()
    write(12,'("NODT=",i0," nt=",i0," mm=",i0," ib=",i0)') ,NODT,nt,mm,ib
    write(12,'("Calculation time=",f10.3,"(sec)")') t2-t1
    write(linebuf,*) 'Date_time='//trim(CALLDATE())//' (gfortran)'
    write (12,'(a)') trim(adjustl(linebuf))
    close(12)

stop

contains

    subroutine SM_FRM(ne,kakom,x,y,Em,AA,AI,sm,hen,NELT,NODT)
        !*************************************************************
        !Element stiffness matrix
        !*************************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,kakom(1:NELT,1:2)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::Em(1:NELT),AA(1:NELT),AI(1:NELT)
        real(8),intent(out)::sm(1:6,1:6),hen(1:6,1:6)
        integer::i,j
        real(8)::AL,EA,EI,cs,sn,w(1:6,1:6)

        !Initialization
        do i=1,6
            do j=1,6
                hen(i,j)=0.0D0
                sm(i,j)=0.0D0
            end do
        end do

        !Transform matrix
        i=kakom(ne,1)
        j=kakom(ne,2)
        AL=sqrt((x(j)-x(i))**2.0D0+(y(j)-y(i))**2.0D0)
        cs=(x(j)-x(i))/AL
        sn=(y(j)-y(i))/AL
        hen(1,1)= cs
        hen(1,2)= sn
        hen(2,1)=-sn
        hen(2,2)= cs
        hen(3,3)=1.0D0
        do i=1,3
            do j=1,3
                hen(i+3,j+3)=hen(i,j)
            end do
        end do

        !Element stiffness matrix
        EA=Em(ne)*AA(ne)
        EI=Em(ne)*AI(ne)
        sm(1,1)=EA/AL
        sm(2,1)=0.0D0
        sm(3,1)=0.0D0
        sm(4,1)=-EA/AL
        sm(5,1)=0.0D0
        sm(6,1)=0.0D0
        sm(2,2)=12.0D0*EI/AL/AL/AL
        sm(3,2)=6.0D0*EI/AL/AL
        sm(4,2)=0.0D0
        sm(5,2)=-12.0D0*EI/AL/AL/AL
        sm(6,2)=6.0D0*EI/AL/AL
        sm(3,3)=4.0D0*EI/AL
        sm(4,3)=0.0D0
        sm(5,3)=-6.0D0*EI/AL/AL
        sm(6,3)=2.0D0*EI/AL
        sm(4,4)=EA/AL
        sm(5,4)=0.0D0
        sm(6,4)=0.0D0
        sm(5,5)=12.0D0*EI/AL/AL/AL
        sm(6,5)=-6.0D0*EI/AL/AL
        sm(6,6)=4.0D0*EI/AL
        do i=1,5
            do j=i+1,6
                sm(i,j)=sm(j,i)
            end do
        end do

        !To convert local coordinates into global coordinates
        do i=1,6
            do j=1,6
                w(i,j)=0.0D0
                do k=1,6
                    w(i,j)=w(i,j)+hen(k,i)*sm(k,j)
                end do
            end do
        end do
        do i=1,6
            do j=1,6
                sm(i,j)=0.0D0
                do k=1,6
                    sm(i,j)=sm(i,j)+w(i,k)*hen(k,j)
                end do
            end do
        end do

    end subroutine SM_FRM


    subroutine GM_FRM(ne,kakom,x,y,Em,AA,AI,sresul,smg,hen,NELT,NODT)
        !*************************************************************
        !Element geometric stiffness matrix
        !*************************************************************
        integer,intent(in)::ne,NELT,NODT
        integer,intent(in)::kakom(1:NELT,1:2)
        real(8),intent(in)::x(1:NODT),y(1:NODT)
        real(8),intent(in)::Em(1:NELT),AA(1:NELT),AI(1:NELT)
        real(8),intent(in)::sresul(1:NELT,1:6)
        real(8),intent(in)::hen(1:6,1:6)
        real(8),intent(out)::smg(1:6,1:6)

        integer::i,j,k
        real(8)::smu(1:6,1:6)
        real(8)::smv(1:6,1:6)
        real(8)::w(1:6,1:6)
        real(8)::AL,P
        real(8)::xx1,yy1,xx2,yy2

        do i=1,6
            do j=1,6
                smv(i,j)=0.0D0
                smu(i,j)=0.0D0
                smg(i,j)=0.0D0
            end do
        end do

        !Axial force
        i=kakom(ne,1)
        j=kakom(ne,2)
        xx1=x(i)
        yy1=y(i)
        xx2=x(j)
        yy2=y(j)
        AL=sqrt((xx2-xx1)*(xx2-xx1)+(yy2-yy1)*(yy2-yy1))
        P=-0.5D0*(sresul(ne,4)-sresul(ne,1)) !圧縮軸力正

        !Geometric stiffness matrix
        smv(1,1)=0.0D0
        smv(2,1)=0.0D0
        smv(2,2)=6.0D0/5.0D0/AL
        smv(3,1)=0.0D0
        smv(3,2)=1.0D0/10.0D0
        smv(3,3)=2.0D0*AL/15.0D0
        smv(4,1)=0.0D0
        smv(4,2)=0.0D0
        smv(4,3)=0.0D0
        smv(4,4)=0.0D0
        smv(5,1)=0.0D0
        smv(5,2)=-6.0D0/5.0D0/AL
        smv(5,3)=-1.0D0/10.0D0
        smv(5,4)=0.0D0
        smv(5,5)=6.0D0/5.0D0/AL
        smv(6,1)=0.0D0
        smv(6,2)=1.0D0/10.0D0
        smv(6,3)=-AL/30.0D0
        smv(6,4)=0.0D0
        smv(6,5)=-1.0D0/10.0D0
        smv(6,6)=2.0D0*AL/15.0D0
        do i=1,6
            do j=i,6
                smv(i,j)=smv(j,i)
            end do
        end do
        smu(1,1)=1.0D0/AL
        smu(4,1)=-1.0D0/AL
        smu(4,4)=1.0D0/AL
        do i=1,6
            do j=i,6
                smu(i,j)=smu(j,i)
            end do
        end do
        do i=1,6
            do j=1,6
                smg(i,j)=smg(i,j)+P*smv(i,j)+P*smu(i,j)
            end do
        end do

        do i=1,6
            do j=1,6
                w(i,j)=0.0D0
                do k=1,6
                    w(i,j)=w(i,j)+hen(k,i)*smg(k,j)
                end do
            end do
        end do
        do i=1,6
            do j=1,6
                smg(i,j)=0.0D0
                do k=1,6
                    smg(i,j)=smg(i,j)+w(i,k)*hen(k,j)
                end do
            end do
        end do

    end subroutine GM_FRM


    subroutine CALST_FRM(ne,kakom,dis,sm,hen,reac,sresul,NELT,NODT,nt)
        !*************************************************************
        !Calculation of stress resultants
        !*************************************************************
        integer,intent(in)::NELT,NODT,nt
        integer,intent(in)::ne,kakom(1:NELT,1:2)
        real(8),intent(in)::dis(1:nt),sm(1:6,1:6),hen(1:6,1:6)
        real(8),intent(inout)::reac(1:nt)
        real(8),intent(out)::sresul(1:NELT,1:6)
        integer :: i,j,k,ia,ir,iw
        real(8) :: wd(1:6),ws(1:6),tempe

        !Nodal displacement
        do i=1,2
            ia=kakom(ne,i)
            wd(3*i-2)=dis(3*ia-2)
            wd(3*i-1)=dis(3*ia-1)
            wd(3*i)=dis(3*ia)
        end do
        !Nodal force in global coordinate system
        do i=1,6
            ws(i)=0.0D0
            do k=1,6
                ws(i)=ws(i)+sm(i,k)*wd(k)
            end do
        end do
        !Stress resultant in local coordinate system
        do i=1,6
            sresul(ne,i)=0.0D0
            do k=1,6
                sresul(ne,i)=sresul(ne,i)+hen(i,k)*ws(k)
            end do
        end do
        !Stress resultant in global coordinate system
        do i=1,6
            ws(i)=0.0D0
            do j=1,6
                ws(i)=ws(i)+hen(j,i)*sresul(ne,j)
            end do
        end do
        !To set nodal force in global coordinate system
        do i=1,2
            ia=kakom(ne,i)
            do j=1,3
                ir=3*(ia-1)+j
                iw=3*(i-1)+j
                reac(ir)=reac(ir)+ws(iw)
            end do
        end do
    end subroutine CALST_FRM


    subroutine ASVEC(ne,nod,nhen,kakom,fevec,ftvec,NELT,NODT)
        !*************************************************************
        !Assembly of vector
        !*************************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,nod,nhen
        integer,intent(in)::kakom(1:NELT,1:nod)
        real(8),intent(in)::fevec(1:nod*nhen)
        real(8),intent(inout)::ftvec(1:NODT*nhen)

        integer :: i,k,ki,ks,kik,ksk

        do i=1,nod
            ki=(kakom(ne,i)-1)*nhen
            ks=(i-1)*nhen
            do k=1,nhen
                kik=ki+k
                ksk=ks+k
                ftvec(kik)=ftvec(kik)+fevec(ksk)
            end do
        end do
    end subroutine ASVEC


    subroutine ASMAT(ne,nod,nhen,kakom,sm,tsm,NELT,NODT)
        !*************************************************************
        !Assembly of matrix
        !*************************************************************
        integer,intent(in)::NELT,NODT
        integer,intent(in)::ne,nod,nhen
        integer,intent(in)::kakom(1:NELT,1:nod)
        real(8),intent(in)::sm(1:nod*nhen,1:nod*nhen)
        real(8),intent(inout)::tsm(1:NODT*nhen,1:NODT*nhen)

        integer :: i,j,k,l,ki,kj,ks,js,kik,kjl,ksk,jsl

        do i=1,nod
            do j=1,nod
                ki=(kakom(ne,i)-1)*nhen
                kj=(kakom(ne,j)-1)*nhen
                ks=(i-1)*nhen
                js=(j-1)*nhen
                do k=1,nhen
                    do l=1,nhen
                        kik=ki+k
                        kjl=kj+l
                        ksk=ks+k
                        jsl=js+l
                        tsm(kik,kjl)=tsm(kik,kjl)+sm(ksk,jsl)
                    end do
                end do
            end do
        end do
    end subroutine ASMAT


    subroutine CHBAND10(n1,m1,array,nt)
        integer,intent(in)::n1,m1,nt
        real(8),intent(inout)::array(1:nt,1:nt)
        integer :: i,j,k,mm,mn
        real(8) :: z

        array(1,1)=sqrt(array(1,1))
        do j=2,m1
            array(1,j)=array(1,j)/array(1,1)
        end do
        do i=2,n1
            if(n1-i+1>m1)then
                mm=m1
            else
                mm=n1-i+1
            end if
            do j=1,mm
                z=array(i,j)
                if(j/=m1)then
                    if(m1-j+1<i)then
                        mn=m1-j+1
                    else
                        mn=i
                    end if
                    do k=2,mn
                        z=z-array(i-k+1,k)*array(i-k+1,j+k-1)
                    end do
                    if(j>1)then
                        array(i,j)=z/array(i,1)
                    else
                        array(i,1)=sqrt(z)
                    end if
                else
                    array(i,j)=z/array(i,1)
                end if
            end do
        end do
    end subroutine CHBAND10


    subroutine CHBAND11(n1,m1,array,vector,nt)
        integer,intent(in)::n1,m1,nt
        real(8),intent(in)::array(1:nt,1:nt)
        real(8),intent(out)::vector(nt)
        integer :: i,j,mm
        real(8) :: z

        vector(1)=vector(1)/array(1,1)
        do i=2,n1
            z=vector(i)
            if(i>m1)then
                mm=m1
            else
                mm=i
            end if
            do j=2,mm
                z=z-array(i-j+1,j)*vector(i-j+1)
            end do
            vector(i)=z/array(i,1)
        end do

        vector(n1)=vector(n1)/array(n1,1)
        do i=n1-1,1,-1
            z=vector(i)
            do j=2,m1
                if(i+j-1>n1)exit
                z=z-array(i,j)*vector(i+j-1)
            end do
            vector(i)=z/array(i,1)
        end do
    end subroutine CHBAND11


    integer function IBAND(mm,array,fact0,nt)
        integer,intent(in)::mm,nt
        real(8),intent(in)::array(1:nt,1:nt),fact0
        integer :: i,j

        IBAND=1
        do i=1,mm-1
            do j=mm,i+1,-1
                if(fact0<=abs(array(i,j)))exit
            end do
            if(IBAND<=j-i+1)IBAND=j-i+1
        end do
        if(mm<IBAND)IBAND=mm
    end function IBAND


    subroutine BAND(mm,ib,array,nt)
        integer,intent(in)::mm,ib,nt
        real(8),intent(inout)::array(1:nt,1:nt)
        integer :: i,j,mn
        real(8) :: work(1:nt)

        do i=2,mm
            if(i+ib-1<mm)then
                mn=ib
            else
                mn=mm-i+1
            end if
            do j=1,ib
                work(j)=0.0D0
            end do
            do k=1,mn
                work(k)=array(i,i+k-1)
            end do
            do j=1,ib
                array(i,j)=work(j)
            end do
        end do
    end subroutine BAND




    subroutine GJACOBI(nf,a,b,ev,vec,nt)
        !************************************
        !Generalized jacobi method
        !************************************
        integer,intent(in)::nf,nt
        real(8),intent(inout)::a(1:nt,1:nt),b(1:nt,1:nt)
        real(8),intent(out)::ev(1:nt),vec(1:nt,1:nt)

        integer::i,j,k
        real(8)::fmk,fmg,epsk,epsg
        real(8)::fact=1.0D-30
        integer::kmax,ip,iq,im
        real(8)::amax,dam,de1,de2,de3,dd,ss,x,gg,ww
        real(8)::amin

        fmk=0.0D0
        fmg=0.0D0
        do i=1,nf
            do j=1,nf
                if(abs(a(i,j))>fmk)fmk=abs(a(i,j))
                if(abs(b(i,j))>fmg)fmg=abs(b(i,j))
            end do
        end do
        epsk=fmk*fact
        epsg=fmg*fact
        do i=1,nf
            do j=1,nf
                vec(i,j)=0.0
            end do
            vec(i,i)=1.0
        end do
        kmax=5*nf*nf
        do k=1,kmax
            ip=1
            iq=2
            amax=a(1,2)**2.0D0/(a(1,1)**2.0D0+a(2,2)**2.0D0)
            do i=1,nf-1
                do j=i+1,nf
                    dam=a(i,j)**2.0D0/(a(i,i)**2.0D0+a(j,j)**2.0D0)
                    if(amax<dam)then
                        ip=i
                        iq=j
                        amax=dam
                    end if
                end do
            end do
            do i=1,nf-1
                do j=i+1,nf
                    dam=b(i,j)**2.0D0/(b(i,i)**2.0D0+b(j,j)**2.0D0)
                    if(amax<dam)then
                        ip=i
                        iq=j
                        amax=dam
                    end if
                end do
            end do
            de1=a(ip,ip)*b(ip,iq)-a(ip,iq)*b(ip,ip)
            de2=a(iq,iq)*b(ip,iq)-a(ip,iq)*b(iq,iq)
            de3=a(ip,ip)*b(iq,iq)-a(iq,iq)*b(ip,ip)
            dd=de3*de3+4.0D0*de1*de2
            ss=sqrt(dd)
            if(de3<0.0D0)then
                x=0.5D0*(de3-ss)
            else
                x=0.5D0*(de3+ss)
            end if
            if(x==0.0D0)then
                dam=0.0D0
                gg=-a(ip,iq)/a(iq,iq)
            else
                dam=de2/x
                gg=-de1/x
            end if
            do j=1,nf
                ww=a(ip,j)
                a(ip,j)=a(ip,j)+gg*a(iq,j)
                a(iq,j)=a(iq,j)+dam*ww
                ww=b(ip,j)
                b(ip,j)=b(ip,j)+gg*b(iq,j)
                b(iq,j)=b(iq,j)+dam*ww
            end do
            do i=1,nf
                ww=a(i,ip)
                a(i,ip)=a(i,ip)+gg*a(i,iq)
                a(i,iq)=a(i,iq)+dam*ww
                ww=b(i,ip)
                b(i,ip)=b(i,ip)+gg*b(i,iq)
                b(i,iq)=b(i,iq)+dam*ww
                ww=vec(i,ip)
                vec(i,ip)=vec(i,ip)+gg*vec(i,iq)
                vec(i,iq)=vec(i,iq)+dam*ww
            end do
            if(amax<epsk.and.amax<epsg)exit
        end do
        do i=1,nf
            ev(i)=a(i,i)/b(i,i)
        end do

        !Sorting
        !ev(k): Minimum eigenvalue: k=1, maximum eigenvalue: k=nf
        !vec(i,k): Eigenvector for (k)th eigenvalue
        do k=1,nf-1
            amin=1.0E+30
            do i=k,nf
                if(ev(i)<=amin)then
                    im=i
                    amin=ev(i)
                end if
            end do
            ww=ev(k)
            ev(k)=ev(im)
            ev(im)=ww
            do i=1,nf
                ww=vec(i,k)
                vec(i,k)=vec(i,im)
                vec(i,im)=ww
            end do
        end do
    end subroutine GJACOBI


    subroutine del_spaces(s)
        character (*), intent (inout) :: s
        character (len=len(s)) tmp
        integer i, j
        j = 1
        do i = 1, len(s)
            if (s(i:i)==' ') cycle
            tmp(j:j) = s(i:i)
            j = j + 1
        end do
        s = tmp(1:j-1)
    end subroutine del_spaces


    real(4) function DIFT()
        integer::count,count_rate,count_max
        call system_clock(count,count_rate,count_max)
        DIFT=real(count)/real(count_rate)
    end function DIFT


    character(len=50) function CALLDATE()
        character(len=8)::date
        character(len=10)::time
        character(len=5)::zone
        integer::values(1:8)
        integer::i
        character(len=4)::ye
        character(len=4)::da
        character(len=10)::mo
        character(len=5)::ti
        character(len=10)::str
        character(len=50)::datime

        call DATE_AND_TIME(date,time,zone,values)
        ye=date(1:4)
        mo=date(5:6)
        da=date(7:8)
        ti=time(1:2)//':'//time(3:4)

        if(trim(adjustl(mo))=='01')mo='January'
        if(trim(adjustl(mo))=='02')mo='February'
        if(trim(adjustl(mo))=='03')mo='March'
        if(trim(adjustl(mo))=='04')mo='April'
        if(trim(adjustl(mo))=='05')mo='May'
        if(trim(adjustl(mo))=='06')mo='June'
        if(trim(adjustl(mo))=='07')mo='July'
        if(trim(adjustl(mo))=='08')mo='August'
        if(trim(adjustl(mo))=='09')mo='September'
        if(trim(adjustl(mo))=='10')mo='October'
        if(trim(adjustl(mo))=='11')mo='November'
        if(trim(adjustl(mo))=='12')mo='December'

        CALLDATE=trim(adjustl(da))//' '//trim(adjustl(mo))//' '//trim(adjustl(ye))//' at '//trim(adjustl(ti))
    end function CALLDATE


end program f90_EGV_BFRM