  !両端固定の梁の問題をとく

       !長さell=１ｍ　荷重位置１ｍ　荷重P=１N　 E=9d9 b1=0.1 h=0.1

       implicit real*8(a-h, o-z)
       dimension a(1000,1000),x(1000),b(1000)
!
      n=102 !節点数*2(自由度、マトリックスの次元)
	  n=6
      call keisan(n,a,x,b)
!メインプログラムはここまで。
      end
!
      subroutine keisan(n,a,x,b)
     implicit real*8(a-h, o-z)
       !dimension a(102,102),x(102),b(102)
       dimension a(n,n),x(n),b(n)
       !
!マトリクス(n,n)
! x(n)は最後に変位を入れるベクトルだけど、
!最初は、境界条件を入れるのに使う。
! 拘束しないところは、1.が入っている。
       do i=1,n
	   !2012/7/26 下のx(i)がx(1)となっていたのを修正
       x(i)=1.d0
       end do
! 境界条件を拘束したい変位成分を変位ベクトルの番号で入れる
! 例えば両端固定の場合
      x(1)=0.d0
      x(2)=0.d0
      x(n-1)=0.d0
      x(n)=0.d0
!	  print*,'x='
!     print*, (x(i),i=1,n)
!       
       n1=n/2-1
       n8=n/2
       n2=2
       n3=-1 
       n4=n-n2
       n5=n1-1
       n6=n1+1
       n7=n6+1
       !
!断面
      b1=0.1
      h=0.1
      h3=h*h*h
      ell=1.
      el=ell/n1
!材料諸元 ヤング率e マトリクスの係数
! print*,'h3,el,n1=',h3,el,n1
      e=9.00d9
      f1=(1.d0/el/el/el)*e*h3*b1 *12.d0
      f2=(1.d0/3/el)*e*h3*b1 *12.d0
      f3=(1.d0/2/el/el)*e*h3*b1 *12.d0
      f4=(1.d0/6/el)*e*h3*b1 *12.d0
!	  print*,'f1,f2,f3,f4=',f1,f2,f3,f4
!荷重ｐ
      p=1.0/n1
      p1=1.0
!理論解
      riron=p1*ell*ell*ell/e/h3/b1*3.0/2.0
!初期化
       do i=1,n
       do j=1,n
          a(i,j)=0.d0
       end do
       end do
!マトリクスの入力
!マトリクス番号：a(i,j)
      !1,2行目まで 
       a(1,1)=f1 !ここは適切に定義
       a(1,2)=n3*f3 !今は0を入れておくけど。
       a(1,3)=n3*f1
       a(1,4)=n3*f3
       a(2,1)=n3*f3
       a(2,2)=f2
       a(2,3)=f3
       a(2,4)=f4
     !3行目からの部分
       do m=0,n8-3
       a(m*2+3,m*2+1)=n3*f1
       a(m*2+3,m*2+2)=f3
       a(m*2+3,m*2+3)=n2*f1
       a(m*2+3,m*2+4)=0.0
       a(m*2+3,m*2+5)=n3*f1
       a(m*2+3,m*2+6)=n3*f3
       a(m*2+4,m*2+1)=n3*f3
       a(m*2+4,m*2+2)=f4
       a(m*2+4,m*2+3)=0.0
       a(m*2+4,m*2+4)=n2*f2
       a(m*2+4,m*2+5)=f3
       a(m*2+4,m*2+6)=f4
       end do
! おしりの2行ぶん
      a(n-1,n-3)=n3*f1
      a(n-1,n-2)=f3
      a(n-1,n-1)=f1
      a(n-1,n)=f3
      a(n,n-3)=n3*f3
      a(n,n-2)=f4
      a(n,n-1)=f3
      a(n,n)=f2
!
!       
!       print*,'a='
!       do i=1,n
!       print'(1p6d11.3)', (a(i,j),j=1,n)
!       print*, (a(i,j),j=1,n)
!       end do
!　境界条件を入れる
 do i=1,n
 do j=1,n
 a(i,j)=x(i)*a(i,j)
 a(j,i)=x(i)*a(j,i)
end do
end do
do i=1,n
 if(x(i)<1.d-3) then
 a(i,i)=1.d0
 end if
end do
!      
! 荷重変位ベクトル
! 両端固定で等分布の場合
      do i=3,n-3,2 !1個とばしに鉛直荷重を入れる
      b(i)=p !p*L/(節点数-2)
       end do
!       
!       print*,'a='
!       do i=1,n
!       print'(1p6d11.3)', (a(i,j),j=1,n)
!       end do
!      !
!     print*,'b='
!     print*, (b(i),i=1,n)
!念の為xを初期化
     do i=1,n
     x(i)=0.d0
     end do
!
    ! 
     call gausu(n,a,x,b)
     print*,'x='
     print*, (x(i),i=1,n)
     print*, (x(i),i=n-1,n)

     print*,riron
     gosa=(x(n-1)-riron)/riron
     print*,gosa

return  !subroutin keisanの終わり
      end
     
     subroutine gausu(n,a,x,b)
     implicit real*8(a-h, o-z)
     dimension a(n,n),x(n),b(n)
     !
     !ガウスの消去法の参考としたのは、
     !名取亮「すうがくぶっくす12 線形計算」(朝倉書店)p.10-15
     !
     do k=1,n-1  !a(k,k)を消去
     do i=k+1,n  !k+1行からn行まで
     do j=k+1,n    !k+1列からn列まで
        a(i,j)=a(i,j)-a(k,j)*a(i,k)/a(k,k)
     end do
        b(i)=b(i)-b(k)*a(i,k)/a(k,k)
     end do
     end do
     !
     ! 後退代入
     x(n)=b(n)/a(n,n)
     do k=n,1,-1
        akjxj=0.d0
     do j=k+1,n
        akjxj=akjxj+a(k,j)*x(j)
     end do
        x(k)=(b(k)-akjxj)/a(k,k)
     end do
     
     return
     end

