8－FEM-VA-beam

　　VA(Vibration  Analysis)(振動解析）

　　8-1-梁の振動解析例

         8-2-  solver(FORTRAN)

         8-3- 全体資料

         8-4- FEM-VA-理論

注（謝辞）

8－3ー全体資料では、梁の固有値、固有モードの記述に際し、

S.Yamauchi　2014年1月28日　機械力学講義ノート（連続体の振動）　Table5　種々のはり曲げ振動を参照させて頂きました。S.Yamauchi先生ありがとうございました。

謝辞

有限要素法による振動解析では、戸川隼人先生の本を参考に

作らせて頂きました。大変に、わかりやすい解説で、理解の

助けになりました。ここにお礼を申し上げます。

また、固有値、固有モードの計算プログラムは、文献（2）を

参照させて頂きました。

文献

（1）戸川隼人著　有限要素法による振動解析　サイエンス社

（2）D.J.ハター著　機械振動解析とプログラミング　

　　　　理工学海外名著シリーズ27　ブレイン図書出版（株）

　　　　東京大学教授　鷲津久一郎　砂川惠　監訳

　　　　日本　I.B.M（株）　高野夏樹　訳

8-1-梁の振動解析例

　梁ー固定ー支持の解析例

ここに出ているFORTRANプログラムをコピー＆ペイストすると、1行おきにブランク行が出来るが、

構わず、******.forに保存して、コンパイル、

実行出来ることを確認しました。

8－1－2．solverをコンパイルして実行

対話例

C:\Users\junsk>cd C:\08_FEM_VA

C:\08_FEM_VA>gfortran FEM_VA_beam_solver.for -o FEM_VA_beam_solver.exe

C:\08_FEM_VA>FEM_VA_beam_solver

 input file name ?

1)FEM_VA_Beam_50ELEM_L0.3_h=0.001_fix_siji.txt

 file read end

 Output file name   OUT

1)FEM_VA_Beam_50ELEM_L0.3_h=0.001_fix_siji-out.txt

 Output file name   Excel

1)FEM_VA_Beam_50ELEM_L0.3_h=0.001_fix_siji-Excel.txt

8－1－3．入力データ

1)FEM_VA_Beam_50ELEM_L0.3_h=0.001_fix_siji.txt

ファイルの中

0.3,0.01,0.001,205.8E+9,7870.,50,3,1.0D-15,1.0D-10

1,2,101

このデータの説明

0．3＝L（梁の全長「ｍ」）

0．01＝ｂ（梁の幅「ｍ」）

0．001＝ｈ（梁の厚み「ｍ」）

205．8E+9＝E（梁のヤング率）

7870．＝ρ（梁の密度「ｋｇ/m^3」）

50＝要素数（梁を50分割）

3＝固定自由度番号数（2行目の数）

1．0D-15＝ERROR収束判定1

1．0D-10＝ERROR収束判定2

1＝X1

2＝θ1

101＝X50

等のように設定する。

8－2．solver

　（以下、Cの書いてある所から下をコピーしてTeraPadに

　　　貼り付け、FEM-VA-beam-solver.forという名称で保存）

C---------------------------------------------------------------

C         FEM   vibration analysis

C

C   solver file name : FEM-VA-beam-solver.for

C---------------------------------------------------------------

C  C:\Users\junsk>cd C:\08_FEM_VA

C  >gfortran FEM_VA_beam_solver.for -o FEM_VA_beam_solver.exe

C  >FEM_VA_beam_solver

C input file name ?

C 1)FEM_VA_Beam_50ELEM_L0.3_h=0.001_fix_siji.txt

C file read end

C

C Output file name   OUT

C 1)FEM_VA_Beam_50ELEM_L0.3_h=0.001_fix_siji-out.txt

C

C Output file name   Excel

C 1)FEM_VA_Beam_50ELEM_L0.3_h=0.001_fix_siji-Excel.txt

C

C

C

C---------------------------------------------------------------

C     variable      mean

C

C     SL(i)         element No.i  L  Beam's Length

C     b(i)          element No.i  B  Beam's width

C     h(i)          element No.i  h  Beam's Hight

C     E(i)          element No.i  E  Beam's Elastic modulus

C     ganma         element No.i  γ Beam's density

C

C     sm            Mass Matrix

C     sk            Stiffness Matrix

C     systemm       System Matrix

C

C--------------------------------------------------------------

C

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

      dimension sm(202,202),sk(202,202),systemm(202,202)

      dimension slamda(202),x(202,202),ww(202,202),R(203,203)

      dimension SL(100),b(100),h(100),ganma(100),E(100),SI(100),A(100)

      DIMENSION V1(202),V2(202),BB(202,202),AA(203,203) 

      DIMENSION CC(202,202),EIG(203),X2(203,203)       

      DIMENSION SUBMX(203,203) 

      dimension sm1(203,203),sk1(203,203)

      dimension U(202,101),V(202,101),IUV(203)

      dimension IRIGD(203),itrace(203)

      character*200 file_name_out,file_name_Excel

      character*200 file_name_in

      write(6,*)' '

      write(6,*)'input file name ?'

      read(5,9000)file_name_in

      OPEN(UNIT=8,FILE=file_name_in,STATUS='OLD')

      read(8,*)QL,SSB,SSH,E1,ganma1,NELEM,NRIGD,ERR1,ERR2

      if(NRIGD.ne.0)then

      read(8,*)(IRIGD(i),I=1,NRIGD)

      else

      endif

C      NELEM=50           ! max 50

      N100=200

      N101=201

      N102=202

      N103=203

      write(6,*)'file read end'

      close(8)

      NELEM=NELEM

C      ERR1=1.0D-15

C      ERR2=1.0D-5

C      NRIGD=2

C      IRIGD(1)=1

C      IRIGD(2)=2

C      IRIGD(2)=101

C      IRIGD(3)=102

      N=NELEM*2+2

      N=N-NRIGD

C      write(6,*)'input L'

C      read(5,*)QL

C      write(6,*)'input b'

C      read(5,*)SSB

C      write(6,*)'input h'

C      read(5,*)SSh

      dal=QL/float(NELEM)

C

C     set up  input data

C

      do 10 i=1,NELEM

      ganma(i)=ganma1

      b(i)=SSB

      h(i)=SSh

      E(i)=E1

      SI(i)=b(i)*h(i)*h(i)*h(i)/12.0

      SL(i)=daL

      A(i)=b(i)*h(i)

   10 continue

C

C     write   output file

C

      write(6,*)' '

      write(6,*)'Output file name   OUT'

      read(5,9000)file_name_out

 9000 format(A)

      write(6,*)' '

      write(6,*)'Output file name   Excel'

      read(5,9000)file_name_Excel

      OPEN(UNIT=3,FILE=file_name_out,STATUS='UNKNOWN')

      write(3,*)' '

      write(3,*)'CANTILEVER vibration analysis'

      write(3,*)' '

      write(3,*)'INPUT DATA'

      write(3,*)' '

      write(3,410)(i,SL(i),i=1,NELEM)

  410 format(5(' L(',i3,')=',g14.6))

      write(3,*)' '

      write(3,420)(i,b(i),i=1,NELEM)

  420 format(5(' B(',i3,')=',g14.6))

      write(3,*)' '

      write(3,430)(i,h(i),i=1,NELEM)

  430 format(5(' H(',i3,')=',g14.6))

      write(3,*)' '

      write(3,440)(i,E(i),i=1,NELEM)

  440 format(5(' E(',i3,')=',g14.6))

      write(3,*)' '

      write(3,450)(i,A(i),i=1,NELEM)

  450 format(5(' A(',i3,')=',g14.6))

      write(3,*)' '

      write(3,460)(i,SI(i),i=1,NELEM)

  460 format(5(' I(',i3,')=',g14.6))

      write(3,*)' '

      write(3,*)' '

      write(3,470)(i,ganma(i),i=1,NELEM)

  470 format(5(' GAM(',i3,')=',g14.6))

      write(3,*)' '

C

C

      call make_m(sm,N102,N102,sm1,N103,N103,

     &ganma,NELEM,A,NELEM,SL,NELEM,NRIGD,IRIGD,N103,

     &N,NELEM)

C      do 5100 i=1,N

C      write(3,5110)(i,j,sm(i,j),j=1,N)

C 5110 format(5(' m(',i3,',',i3,')=',e14.6))

C 5100 continue

C

      call make_k(sk,N102,N102,sk1,N103,N103,

     &E,NELEM,SI,NELEM,SL,NELEM,NRIGD,IRIGD,N103,

     &N,NELEM)

C      do 5102 i=1,N

C      write(3,5112)(i,j,sk(i,j),j=1,N)

C 5112 format(5(' k(',i3,',',i3,')=',e14.6))

C 5102 continue

C

C     [ww]=inverse matrix[sm]-1     R(N102,N102) is work area  

C

      call GINV(R,N103,N103,N,N102,ww,sm,ERR1,ILL) 

C

C     [systemm]=[ww]*[sk]

C

      call MMUL2(systemm,ww,sk,N102,N102,N)

C

      write(6,*)'EIGEN START'

C

C

C

C

      call EIGEN(systemm,N102,N102,slamda,N102,x,N102,N102,N,

     &ERR2,ILL,V1,V2,BB,CC,AA,X2,SUBMX,EIG,N102,N103)

C

      if(ILL.eq.1)then

      write(6,*)'eigen err'

      go to 30

      else

      endif

C

C　　　固有値と固有ベクトルが次数にたいして順番が逆になっているので

C　　　入れ替えを行う。

C                                        

      DO 1070 I=1,NELEM   

      J=NELEM*2+2-NRIGD-I+1            

      AAA=slamda(I)      

      slamda(I)=slamda(J)   

      slamda(J)=AAA  

      DO 1080 K=1,N   

      AAA=x(K,I)        

      x(K,I)=x(K,J)    

      x(K,J)=AAA      

 1080 CONTINUE

 1070 CONTINUE

C

C     RIGD elimination  hukkatsu

C

      write(6,*)'RIGD   hukkatsu'

      Iskip=0

      do 2200 INEW=1,NELEM*2+2

      do 2210 IOLD=1,NRIGD

      if(IRIGD(IOLD).eq.INEW)then

      IUV(INEW)=0

      Iskip=Iskip+1

      else

      IUV(INEW)=INEW-Iskip

      endif

 2210 continue

 2200 continue

C

      do 2120 K=1,N

      ii=0

      do 2100 INEW=1,NELEM*2+2

      if(IUV(INEW).eq.0)then

      AA(K,INEW)=0.0

      else

      AA(K,INEW)=x(IUV(INEW),K)

      endif

 2100 continue

 2120 continue

C

C     奇数番号：U

C     偶数番号：θ

C

      write(6,*)' U   V   huruiwake'

      do 2150 K=1,N

      ii=0

      do 2130 i=1,NELEM*2+2,2

      ii=ii+1

      U(K,ii)=AA(K,i)

 2130 continue

      ii=0

      do 2110 i=2,NELEM*2+2,2

      ii=ii+1     

      V(K,ii)=AA(K,i)

 2110 continue

 2150 continue

C

C

      write(6,*)'eigen vevtor   normalies'

      do 2060 K=1,N

      umax=0.0

      vmax=0.0

      do 2040 i=1,NELEM+1

      if(DABS(U(k,i)).gt.umax)then

      umax=DABS(U(k,i))

      else

      endif

      if(DABS(V(k,i)).gt.vmax)then

      vmax=DABS(V(K,i))

      else

      endif

 2040 continue

      do 2050 i=1,NELEM+1

      U(k,i)=U(K,i)/umax

      V(K,i)=V(K,i)/vmax

 2050 continue

 2060 continue

C

C     write   result

C

      write(3,*)' '

      write(3,*)'cantilever analysis  result'

      write(3,*)' '

      write(3,100)(i,SQRT(slamda(i))/2.0/3.14159265,i=1,N)

  100 format('natural freq.',i3,'=',G16.8)

      do 20 k=1,N

      write(3,*)' '

      if(k.eq.1)then

      write(3,*)'primary vibration n=',k,'  mode'

      else

      write(3,*)'higher-order vibration n=',K,'  mode'

      endif

      WRITE(3,200)(i,U(k,i),i=1,NELEM+1)

  200 FORMAT(5(' U ',I3,'=',E16.8))

      WRITE(3,300)(i,V(k,i),i=1,NELEM+1)

  300 FORMAT(5(' V ',I3,'=',E16.8))

   20 continue

      CLOSE(3)

C

C      write Excel Data

C

      OPEN(UNIT=2,FILE=file_name_Excel,STATUS='UNKNOWN')

      SX=0.0

      write(2,510)(i,SQRT(slamda(i))/2.0/3.14159265,i=1,NELEM*2)

  510 format('X',202(',',i3,'mode FEM ',f12.3,'Hz'))

      SX=0.0

      do 530 i=1,NELEM+1

      write(2,540)SX,(U(K,i),k=1,N)

      SX=SX+SL(i)

  540 format(e14.6,202(',',E14.6))

  530 continue

      CLOSE(2)

      write(6,*)'calculation  stop end'

   30 stop

      end

C

C------------------------------------------------------------------

C

      subroutine make_k(sk,n1,n2,sk1,n6,n7,E,n3,SI,n4,SL,n5,

     &NRIGD,IRIGD,n8,NN,NELEM)

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

      dimension sk(n1,n2),E(n3),SI(n4),SL(n5)

      dimension Q(4,4)

      dimension sk1(n6,n7)

      dimension IRIGD(n8)

      write(3,*)'make_k   in'

      do 80 i=1,n6

      do 90 j=1,n7

      sk1(i,j)=0.0

   90 continue

   80 continue

      do 82 i=1,n1

      do 92 j=1,n2

      sk(i,j)=0.0

   92 continue

   82 continue

C      write(3,110)(i,SL(i),i=1,NELEM)

C  110 format(9(' SL(',i2,')=',E16.6))

      do 10 N=1,NELEM

      SL3=SL(N)*SL(N)*SL(N)

      SL2=SL(N)*SL(N)

      S=E(N)*SI(N)        

      Q(1,1)=12.*S/SL3    

      Q(1,2)=-6.*S/SL2    

      Q(1,3)=-Q(1,1)      

      Q(1,4)=Q(1,2)       

      Q(2,1)=Q(1,2)       

      Q(2,2)=4.*S/SL(N)  

      Q(2,3)=6.*S/SL2    

      Q(2,4)=2.*S/SL(N) 

      Q(3,1)=Q(1,3)     

      Q(3,2)=Q(2,3)     

      Q(3,3)=Q(1,1)     

      Q(3,4)=Q(2,3)     

      Q(4,1)=Q(1,4)     

      Q(4,2)=Q(2,4)     

      Q(4,3)=Q(3,4)     

      Q(4,4)=Q(2,2)

C

C    

      if(N.eq.1)then

      DO 30 I=1,4               

      DO 40 J=1,4             

      sk1(I,J)=Q(I,J)  

   40 CONTINUE         

   30 CONTINUE      

      else

      L=N*2-2   

      do 32 I=1,4

      IL=I+L

      do 42 J=1,4

      JL=J+L

      sk1(IL,JL)=sk1(IL,JL)+Q(I,J)

   42 continue

   32 continue

      endif

   10 continue

C

C     RIGD elemination

C

      if(NRIGD.eq.0)go to 350

      do 34 m=1,NRIGD

      mb=NRIGD-m+1      !  NRIGD  from back trace No. is  mb

      write(3,*)'mb=',mb,' IRIGD(mb)=',IRIGD(mb)

      if(IRIGD(mb).eq.1)then

      do 36 i=1,NELEM*2+1

      IL=1

      do 46 j=1,NELEM*2+1

      JL=1

      sk1(i,j)=sk1(i+IL,j+JL)

   46 continue

   36 continue

      else

      write(3,*)'else  mb=',mb,' IRIGD(mb)=',IRIGD(mb)

      NUKI=IRIGD(mb)

      do 262 j=NUKI,NELEM*2+2

      do 120 i=1,NUKI-1

      sk1(i,j)=sk1(i,j+1)

  120 continue

  262 continue

      do 132 i=NUKI,NELEM*2+2

      do 130 j=1,NUKI-1

      sk1(i,j)=sk1(i+1,j)

  130 continue

  132 continue

      do 140 i=NUKI,NELEM*2+2

      do 142 j=NUKI,NELEM*2+2

      sk1(i,j)=sk1(i+1,j+1)

  142 continue

  140 continue

      endif

   34 continue

C

C

  350 do 270 i=1,NN

      do 280 j=1,NN

      sk(i,j)=sk1(i,j)

  280 continue

  270 continue

      write(6,*)'make_k   return'

      return

      end

C

C-----------------------------------------------------------------

C

      subroutine make_m(sm,n1,n2,sm1,n6,n7,ganma,n3,A,n4,SL,n5,

     &NRIGD,IRIGD,n8,NN,NELEM)

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

      dimension sm(n1,n2),ganma(n3),A(n4),SL(n5)

      dimension Q(4,4)

      dimension sm1(n6,n7)

      dimension IRIGD(n8)

      write(3,*)'make_m     in'

      do 80 i=1,n6

      do 90 j=1,n7

      sm1(i,j)=0.0

   90 continue

   80 continue

      do 82 i=1,n1

      do 92 j=1,n2

      sm(i,j)=0.0

   92 continue

   82 continue

      do 10 N=1,NELEM

      W=A(N)*SL(N)*ganma(N) 

      Q(1,1)=13./35.*W      

      Q(1,2)=-11.*SL(N)/210.*W

      Q(1,3)=9./70.*W                    

      Q(1,4)=13./420.*SL(N)*W

      Q(2,1)=-11.*SL(N)/210.*W  

      Q(2,2)=SL(N)*SL(N)/105.*W   

      Q(2,3)=-13.*SL(N)/420.*W   

      Q(2,4)=-SL(N)*SL(N)/140.*W 

      Q(3,1)=Q(1,3)      

      Q(3,2)=Q(2,3)  

      Q(3,3)=Q(1,1)    

      Q(3,4)=-Q(2,1)     

      Q(4,1)=Q(1,4)      

      Q(4,2)=Q(2,4)

      Q(4,3)=Q(3,4)      

      Q(4,4)=Q(2,2)

C

C    

      if(N.eq.1)then

      DO 30 I=1,4               

      DO 40 J=1,4             

      sm1(I,J)=Q(I,J)  

   40 CONTINUE         

   30 CONTINUE      

      else

      L=N*2-2   

      do 32 I=1,4

      IL=I+L

      do 42 J=1,4

      JL=J+L

      sm1(IL,JL)=sm1(IL,JL)+Q(I,J)

   42 continue

   32 continue

      endif

   10 continue

C

C     RIGD elemination

C

C

      write(6,*)'make_m   RIGD  elimination'

      do 34 m=1,NRIGD

      write(6,*)'RIGD   m=',m

      mb=NRIGD-m+1

      if(IRIGD(mb).eq.1)then

      do 36 i=1,NELEM*2+1

      IL=1

      do 46 j=1,NELEM*2+1

      JL=1

      sm1(i,j)=sm1(i+IL,j+JL)

   46 continue

   36 continue

      else

      NUKI=IRIGD(mb)

      do 262 j=NUKI,NELEM*2+2

      do 120 i=1,NUKI-1

      sm1(i,j)=sm1(i,j+1)

  120 continue

  262 continue

      do 132 i=NUKI,NELEM*2+2

      do 130 j=1,NUKI-1

       sm1(i,j)=sm1(i+1,j)

  130 continue

  132 continue

      do 140 i=NUKI,NELEM*2+2

      do 142 j=NUKI,NELEM*2+2

      sm1(i,j)=sm1(i+1,j+1)

  142 continue

  140 continue

      endif

   34 continue

C

C

      do 270 i=1,NN

      do 280 j=1,NN

      sm(i,j)=sm1(i,j)

  280 continue

  270 continue

      write(6,*)'make_m    return'

      return

      end

C

C

C------------------------------------------------------------------

C

C　　　マトリックスの掛け算　[A] * [B] = [C]

C

      SUBROUTINE MMUL2(C,A,B,L1,L2,N)

      REAL*8 C,A,B

      DIMENSION C(L1,L1),A(L1,L2),B(L2,L1)

      DO 10 I=1,L1

      DO 12 J=1,L1

      C(I,J)=0.

   12 CONTINUE

   10 CONTINUE

      DO 20 I=1,N

      DO 22 J=1,N

      DO 24 K=1,N

      C(I,J)=C(I,J)+A(I,K)*B(K,J)

   24 CONTINUE

   22 CONTINUE

   20 CONTINUE

      RETURN

      END 

C------------------------------------------------------------------

C　　　固有値、固有ベクトルを求めるサブルーチン

C　　　　　X(I,J)　　　　 ：対象となるマトリックス

C　　　　　VALU(I)　　　　：固有値

C　　　　　VEX(I,J)　　　 ：固有ベクトル

C　　　　　ILL　　　　　  ：ILL=1の時エラー

C

      SUBROUTINE EIGEN(X,NA,NB,VALU,NC,VEX,ND,NE,N,ERR,ILL,

     &V1,V2,BB,CC,AA,X2,SUBMX,EIG,N102,N103)

      REAL*8 X,VALU,VEX,ERR,V1,V2,BB,AA,CC,EIG,X2,SUBMX    

      REAL*8 AAA,ANORM,B        

      DIMENSION X(NA,NB),VALU(NC),VEX(ND,NE)      

      DIMENSION V1(N102),V2(N102),BB(N102,N102),AA(N103,N103) 

      DIMENSION CC(N102,N102),EIG(N103),X2(N103,N103)       

      DIMENSION SUBMX(N103,N103) 

      write(6,*)'EIGEN IN'  

      NN=N

      DO 2206 I=1,ND

      DO 2206 J=1,NE

      VEX(I,J)=0.0

 2206 CONTINUE

      DO 2208 I=1,NC

      VALU(I)=0.0

 2208 CONTINUE

      DO 2100 I=1,N103

      EIG(I)=0.

      DO 2100 J=1,N103

      X2(I,J)=0.

      AA(I,J)=0.

      SUBMX(I,J)=0.

 2100 CONTINUE

      DO 2200 I=1,N102

      V1(I)=0

      V2(I)=0.

      DO 2200 J=1,N102

      BB(I,J)=0.

      CC(I,J)=0.

 2200 CONTINUE        

      DO 10 I=1,N  

      DO 10 J=1,N  

   10 X2(I,J)=X(I,J)

      LL=N-1

      DO 1000 L=1,LL   

      CALL ITER(X2,NN,ERR,B,EIG,MAQ,KK,N103)

C

      IF(MAQ)2000,2000,2001     

 2001 WRITE(6,2002)L          

 2002 FORMAT(1H ,"NO CONV. AFTER 10000 LOOPS",   

     &" CRASH ON VECTOR NO.",I5)           

      GO TO 500                

 2000 CONTINUE            

      IF(L-1)3000,3000,3001   

 3000 DO 3002 K=1,N         

 3002 VEX(K,1)=EIG(K)       

 3001 VALU(L)=B

C

      DO 4000 I=1,NN        

      DO 4000 J=1,NN     

 4000 SUBMX(I,J)=EIG(I)*X2(KK,J) 

      DO 5000 I=1,NN        

      DO 5000 J=1,NN      

 5000 SUBMX(I,J)=X2(I,J)-SUBMX(I,J) 

      NN=NN-1                       

      II=0                       

      DO 6000 I=1,NN        

   25 II=II+1             

      IF(II-KK)26,25,26    

   26 CONTINUE           

      JJ=0               

      DO 6000 J=1,NN     

   28 JJ=JJ+1            

      IF(JJ-KK)27,28,27  

   27 X2(I,J)=SUBMX(II,JJ)  

 6000 CONTINUE          

 1000 CONTINUE         

      VALU(N)=X2(1,1)

      write(6,*)'VALU(',N,')=',N,'f=',sqrt(VALU(N))/2./3.14159

      write(3,*)'VALU(',N,')=',N,'f=',sqrt(VALU(N))/2./3.14159

C

      NK=N-1         

      DO 888 L=2,N   

      DO 37 I=2,N    

   37 V1(I-1)=-X(I,1) 

      DO 1500 I=1,NK  

      DO 1500 J=1,NK  

 1500 SUBMX(I,J)=X(I+1,J+1)  

      DO 750 M=2,N          

  750 SUBMX(M-1,M-1)=X(M,M)-VALU(L) 

      do 761 M=1,N-1      

      do 760 MM=1,N-1     

      CC(M,MM)=SUBMX(M,MM)

  760 continue

  761 continue

      CALL GINV(AA,N103,N103,N-1,N102,BB,CC,ERR,ILL)    

      CALL GMPRD(BB,V1,V2,N102,N102,N-1)       

      IFLG=0                 

      B=1.0                  

      BX=1.0                 

      VEX(1,L)=1.0           

      DO 850 K=1,NK        

      IF(BX-DABS(V2(K)))900,850,850   

  900 B=V2(K)               

      BX=ABS(V2(K))         

      IFLG=1              

  850 CONTINUE            

      IF(B.LT.0.0)B=-B    

      IF(IFLG)875,875,876 

  876 VEX(1,L)=VEX(1,L)/B 

      DO 69 K=1,NK        

   69 V2(K)=V2(K)/B       

  875 DO 877 I=2,N      

  877 VEX(I,L)=V2(I-1)

C

C     

      f1=SQRT(VALU(N-L+1))/2./3.14159265358979

      write(6,*)'L',L,'MODE VECTOR ',N-L+1,' freq=',f1,'Hz'

      write(3,*)'L',L,'MODE VECTOR ',N-L+1,' freq=',f1,'Hz'

  888 CONTINUE           

C

  500 CONTINUE

C

      DO 1272 K=1,NN

      K1=K+1

      DO 1270 J=K,NN

      AAA=0.0

      DO 1242 I=1,NN

      AAA=AAA+VEX(K,I)*VEX(J,I)

 1242 CONTINUE

 1270 CONTINUE

 1272 CONTINUE

      DO 1274 K=1,NN

      K1=K+1

      DO 1276 J=K,NN

      AAA=0.0

      DO 1243 I=1,NN

      AAA=AAA+VEX(I,K)*VEX(I,J)

 1243 CONTINUE

      IF(K.EQ.J)ANORM=AAA

      AAA=AAA/ANORM

 1276 CONTINUE

 1274 CONTINUE            

      write(6,*)'EIGEN  OUT' 

      RETURN                

      END

C

C------------------------------------------------------------------

C                                                             

      SUBROUTINE ITER(X,N,ERR,B,EIG,MAQ,KK,N103)

      REAL*8 X,ERR,B,EIG,ANUD  

C-------THIS SUBROUTINE CARRIES OUT THE ITERATION         

      DIMENSION X(N103,N103),EIG(N103),ANUD(N103)

      DO 1000 I=1,N103

      ANUD(I)=0.

 1000 CONTINUE       

      L=0            

      MAQ=0         

C-------SET UP VECTOR OF UNITY AS FIRST TEST       

      DO 1 K=1,N    

      EIG(K)=1.0          

    1 CONTINUE            

C-------CALCULATE PRODUCT OF MATRIX AND TEST VECTOR

    7 DO 2 J=1,N         

      ANUD(J)=0.0       

      DO 2 K=1,N            

      ANUD(J)=ANUD(J)+X(J,K)*EIG(K)     

    2 CONTINUE  

      B=ANUD(1)          

      KK=1               

C-------FIND LARGEST ELEMENT AND NORMALISE    

      DO 30 J=2,N        

      IF(DABS(B)-DABS(ANUD(J)))31,30,30       

   31 B=ANUD(J)      

      KK=J        

   30 CONTINUE

      IF(DABS(B).LT.1.D-300)GO TO 1020

       GO TO  1010

 1020 WRITE(3,100)B,L

      WRITE(6,100)B,L

  100 FORMAT(1H ,"ITER B.LT.1.D-300 B=",E18.8,3X,"L=",I4)

      GO TO 20

 1010 L=L             

      DO 3 J=1,N         

      ANUD(J)=ANUD(J)/B  

    3 CONTINUE           

C-------COMPARE NEW VECTOR WITH TEST VECTOR    

      DO 5 J=1,N                

      IF(DABS(EIG(J)-ANUD(J))-ERR)5,5,4    

    5 CONTINUE      

      GO TO 9       

C-------REPLACE TESR VECTOR BY NEW VECTOR 

    4 DO 6 J=1,N     

      EIG(J)=ANUD(J) 

    6 CONTINUE       

      L=L+1  

C-------TEST FOR 10000 ITERATIONS   

      IF(10000-L)20,20,21     

   21 GO TO 7    

   20 MAQ=1      

    9 CONTINUE   

      RETURN     

      END      

C

C------------------------------------------------------------------

C                                                       

      SUBROUTINE GMPRD(BB,V1,V2,K1,K2,NK)

      REAL*8 BB,V1,V2    

      DIMENSION BB(K1,K1),V1(K2),V2(K2)  

      DO 10 I=1,K2         

      V2(I)=0.             

      DO 20 J=1,NK         

      V2(I)=V2(I)+BB(I,J)*V1(J)   

   20 CONTINUE         

   10 CONTINUE         

      RETURN          

      END        

C    *******************************************************

C    ***  THIS PROGRAM COMPUTES                          ***

C    ***    THE INVERS OF A MATRIX                       ***

C    ***    USING THE GAUSS ELIMINATION METHOD           ***

C    *******************************************************

C                   

      SUBROUTINE GINV(A,LL,KK,IB,IC,B,C,ESP,ILL)  

      REAL*8 A,B,C,ESP,BMAX,AA,CC            

      DIMENSION A(LL,KK),B(IC,IC),C(IC,IC),JB(1000)

      ILL=0      

      IB1=IB+1  

      DO 90 III=1,IB 

      DO 95 JJJ=1,IB 

      A(JJJ,IB1)=0.  

   95 CONTINUE       

      A(III,IB1)=1.  

      DO 45 IIJ=1,IB 

      DO 45 JJI=1,IB 

   45 A(IIJ,JJI)=C(IIJ,JJI) 

      DO 15 I=1,IB   

   15 JB(I)=I        

      DO 10 K=1,IB-1 

      BMAX=DABS(A(K,K)) 

      IMAX=K    

      JMAX=K    

      DO 20 I=K,IB 

      DO 20 J=K,IB 

      AA=DABS(A(I,J))  

      IF(BMAX.GE.AA)GO TO 20  

      BMAX=AA  

      IMAX=I   

      JMAX=J   

   20 CONTINUE 

      DO 30 J=1,IB1 

      AA=A(IMAX,J)  

      A(IMAX,J)=A(K,J)

      A(K,J)=AA 

   30 CONTINUE  

      KB=JB(K) 

      JB(K)=JB(JMAX) 

      JB(JMAX)=KB   

      DO 35 I=1,IB  

      AA=A(I,JMAX)  

      A(I,JMAX)=A(I,K) 

      A(I,K)=AA    

   35 CONTINUE     

      DO 40 I=K+1,IB 

      CC=A(I,K)/A(K,K) 

      DO 40 J=1,IB1    

      A(I,J)=A(I,J)-CC*A(K,J)

   40 CONTINUE  

   10 CONTINUE  

      DO 51 KKK=1,IB-1

      K=IB1-KKK 

      DO 50 II=1,K-1 

      I=K-II     

      CC=A(I,K)/A(K,K)   

      DO 50 J=1,IB1      

      A(I,J)=A(I,J)-CC*A(K,J)

   50 CONTINUE       

   51 CONTINUE       

      DO 60 K=1,IB   

      AA=DABS(A(K,K)) 

      IF(AA.LT.ESP) GO TO 70  

   60 A(K,IB1)=A(K,IB1)/A(K,K) 

      DO 65 I=1,IB 

      DO 66 J=I,IB 

      KB=JB(J)   

      IF(KB.EQ.I) GO TO 67  

   66 CONTINUE            

   67 JB(J)=JB(I)         

      AA=A(KB,IB1)        

      A(KB,IB1)=A(J,IB1)  

      A(J,IB1)=AA         

   65 CONTINUE           

      DO 25 I=1,IB       

      B(I,III)=A(I,IB1)  

   25 CONTINUE     

   90 CONTINUE    

      GO TO 80    

   70 ILL=1       

   80 RETURN      

      END