..................の間がプログラム本体です
program 10
..................
C 行列演算の応用
C Gauss の消去法(ピボット選択なし)
C  a x = b の解 x を求める

      program main
      parameter(nn=3)
      real*8 a(nn,nn),b(nn),x(nn)
      real*8 awork(nn,nn),bwork(nn),check(nn)
      real*8 dd,xx

C  ファイルから行列aに値を読み込む
      do i=1,nn
        read(*,*)(a(i,j),j=1,nn)
      enddo

C  ファイルからベクトルbに行列を読み込む
      read(*,*)(b(i),i=1,nn)

C  作業用の行列aworkとベクトルbworkにaとbをコピー
      do i=1,nn
        do j=1,nn
          awork(i,j)=a(i,j)
        enddo
      enddo
      do i=1,nn
         bwork(i)=b(i)
      enddo

C  前進消去
      do i=2,nn
        do j=i,nn
          dd=awork(j,i-1)/awork(i-1,i-1)
          awork(j,i-1)=0
          bwork(j)=bwork(j)-dd*bwork(i-1)
          do k=i,nn
            awork(j,k)=awork(j,k)-dd*awork(i-1,k)
          enddo
        enddo
      enddo

C  消去結果の出力(右上三角行列になったことを確認)
      write(*,*)"result of forward elimination"
      write(*,*)"matrix awork ="
      do i=1,nn
        write(*,*)(awork(i,j),j=1,nn)
      enddo
      write(*,*)"vector bwork ="
      write(*,*)(bwork(i),i=1,nn)

C  後退代入
      do i=nn,1,-1
        xx=0
        do j=i+1,nn
            xx=xx+x(j)*awork(i,j)
        enddo
        x(i)=(bwork(i)-xx)/awork(i,i)
      enddo

C  代入結果の出力
      write(*,*)"result of backward substitution"
      write(*,*)"vector x ="
      write(*,*)(x(i),i=1,nn)

C  結果のチェック( a x  を計算)
      do i=1,nn
        check(i)=0
        do j=1,nn
          check(i)=check(i)+a(i,j)*x(j)
       enddo
      enddo        

C  チェック結果の出力(ax-b=0であることを確認)
      write(*,*)"result of consistency check"
      write(*,*)"a x - b ="
      write(*,*)(check(i)-b(i),i=1,nn)

      end
..................
program 11
..................
C 関数サブプログラムの利用
C 2分法による非線形方程式の解
      program main
      parameter(nmax=30)
      real*8 x1,x2,xtemp
      real*8 func

C 解探索の範囲を指定
      read(*,*)x1,x2

C 2分法を最大nmax回くりかえす
      do irepeat = 1,nmax

C  x1とx2の差が与えた精度以内なら終了
        if(abs((x1-x2)/x1).le. 1e-10)then
          write(*,*)"solution is between ",x1," and ",x2
          stop

C  func(x1)とfunc(x2)の符号が逆なら範囲を二分する
        else if(func(x1)*func(x2).lt.0)then
            xtemp = (x1+x2)/2
            if(func(x1)*func(xtemp).lt.0)then
              x2=xtemp
            else
              x1=xtemp
            endif
       write(*,*)irepeat," step  x1=",x1,"  x2=",x2

C  func(x1)とfunc(x2)が同符号なら、この範囲に解がないと判断して終了
        else
          write(*,*)"no solution is found between ",x1," and ",x2
          stop
        endif
      enddo

C  nmax回のくりかえしで解に収束しなかったときは終了
      write(*,*)"not converged after ",nmax," steps"
      write(*,*)"current values of x are ",x1," and ",x2
      end

C  解を得たい関数
      function func(x)
      real*8 func
      real*8 x
      func= x**3 + 3*x**2 + 5*x + 1
      end