cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c 付録 O の主な例題のソースコード
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program area1
c                            ex1: 四角形の面積を求める
c                            縦と横の長さを入力し，面積を出力
c ← １桁目が``c''の行はコメント行（プログラムとは無関係な文）
c    | ← 1〜5桁が行番号の桁
c     |←　7桁目から72桁目までがプログラム本体（73桁以降は無視）     →|
      read(5,100) a,b
      write(6,200) a,b
      a=a*b
      write(6,300) a
      stop
c
  100 format(2f10.0)
  200 format(' ','a =',1pe15.7,'  b =',e15.7)
  300 format(' ','A =',1pe15.7)
c                            5, 6は標準入出力装置番号
c                            キーボードと画面が5と6，あるいは
c                            area1  <inarea.dat   >outarea.txt
c                            ファイル`inarea.dat'の例
c                            ----+----1----+----2
c                                   1.2       2.7
c                            ----+----1----+----2
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c input-data file `inarea.dat'
       1.2       2.7
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program area11
c                            ex1_1: 四角形の面積と
c                                             断面2次モーメントを求める
c                            縦と横の長さを入力
      read(5,100) a,b
      write(6,200) a,b
c                            A= a*b
c                            I= a*a*a*b/12.
      b=a*b
      a=a*a*b/12.
      write(6,300) b,a
      stop
c
  100 format(2f10.0)
  200 format(' ','a =',1pe15.7,'  b =',e15.7)
  300 format(' ','A =',1pe15.7,'  I =',e15.7)
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program area12
c                            ex1_2: 四角形の面積と
c                                             断面2次モーメントを求める
c                            縦と横の長さを入力
      read(5,100) a,b
      write(6,200) a,b
c                            A= a*b
c                            I= a*a*a*b/12.= a**3*b/12. べきには ** を使う
      area=a*b
      I=a**3*b/12.
c                       ↓実行するとエラー
      write(6,300) area,I
      stop
c
  100 format(2f10.0)
  200 format(' ','a =',1pe15.7,'  b =',e15.7)
  300 format(' ','A =',1pe15.7,'  I =',e15.7)
c                                       ↑頭文字が`i'から`n'までの変数は
c                                         通常は整数として解釈される
c                一つの解決策: program文の次に`real i'の１行を入れる
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program area13
c                            ex1_3: 長い計算
      read(5,100) a,b
      write(6,200) a,b
c
      x = a * b**3 + b * sin(a)
      y = x * ( a * cos(b) - exp(a) / 2.4 ) / ( 120.9 * tan( a/b ) -
     &    23.03 ) - x
c    ↑
c    6桁目に `&' を入れると，継続行とみなされる。
      write(6,300) x,y
      stop
c
  100 format(2f10.0)
  200 format(' ','a =',1pe15.7,'  b =',e15.7)
  300 format(' ','x =',1pe15.7,'  y =',e15.7)
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program area2
      open(5,file='inarea.dat')
      open(6,file='outarea.txt')
c
      read(5,100) a,b
      write(6,200) a,b
      a=a*b
      write(6,300) a
c
      close (5)
      close (6)
      stop
c
  100 format(2f10.0)
  200 format(' ','a =',1pe15.7,'  b =',e15.7)
  300 format(' ','A =',1pe15.7)
c                            5, 6を明確にファイルに割り当てる
c                            入力: inarea.dat   出力: outarea.txt
c                            ファイル`inarea.dat'の例
c                            ----+----1----+----2
c                                   1.2       2.7
c                            ----+----1----+----2
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c input-data file `inarea.dat'
       1.2       2.7
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program area3
      open(5,file='inarea.dat')
      open(6,file='outarea.txt')
      read(5,100) a,b
c ----------------- ここから
c                            連続入力を，あり得ない入力値で制御する
      do while(a.gt.0.)
          write(6,200) a,b
          a=a*b
          write(6,300) a
c                            次の入力をして
          read(5,100) a,b
      end do
c ----------------- ここまでがループ
      close (5)
      close (6)
      stop
c
  100 format(2f10.0)
  200 format(' ','a =',1pe15.7,'  b =',e15.7)
  300 format(' ',15x,'A =',1pe15.7)
c                            5, 6を明確にファイルに割り当てる
c                            入力: inarea.dat   出力: outarea.txt
c                            ファイル`inarea.dat'の例
c                            ----+----1----+----2
c                                   1.2       2.7
c                                   2.3       8.1
c                                   -1.       8.1
c                            ----+----1----+----2
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c input-data file `inarea.dat'
       1.2       2.7
       2.3       8.1
       -1.       8.1
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program area4
c                            ex4: 円の面積を求める
c                            半径を入力し，面積を出力
      data pi / 3.14159265 /
c
      open(5,file='inarea.dat')
      open(6,file='outarea.txt')
      read(5,100) r
c ----------------- ここから
      do while(r.gt.0.)
          write(6,200) r
          a=pi*r*r
          write(6,300) a
c                            次の入力をして
          read(5,100) r
      end do
c ----------------- ここまで
      close (5)
      close (6)
      stop
c
  100 format(f10.0)
  200 format(' ','radius =',1pe15.7)
  300 format(' ',15x,'A =',1pe15.7)
c                            ファイル`inarea.dat'の例
c                            ----+----1----+----2
c                                   1.2
c                                   2.3
c                                   -1.
c                            ----+----1----+----2
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c input-data file `inarea.dat'
       1.2
       2.3
       -1.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program area5
      data pi / 3.14159265 /
c
      open(5,file='inarea.dat')
      open(6,file='outarea.txt')
c
      read(5,100) a,b
c
      do while(a.gt.0.)
c                            入力の数で円か四角か判断
          if(b.gt.0.)then
c                            四角
              write(6,400) a,b
              a=a*b
              write(6,300) a
          else
c                            円
              write(6,200) a
              a=pi*a*a
              write(6,300) a
c
          end if
          read(5,100) a,b
      end do
c
      close (5)
      close (6)
      stop
c
  100 format(2f10.0)
  200 format(' ','radius =',1pe15.7)
  300 format(' ',15x,'A =',1pe15.7)
  400 format(' ','a =',1pe15.7,'  b =',e15.7)
c                            ファイル`inarea.dat'の例
c                            ----+----1----+----2
c                                   1.2        0.
c                                   2.3       4.5
c                                   3.6       .45
c                                    .2        0.
c                                   -1.       .45
c                            ----+----1----+----2
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c input-data file `inarea.dat'
       1.2        0.
       2.3       4.5
       3.6       .45
        .2        0.
       -1.       .45
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program area6
      common / const / pi
      pi=4.*atan(1.)
c
      open(5,file='inarea.dat')
      open(6,file='outarea.txt')
c
      read(5,100) a,b
      do while(a.gt.0.)
          if(b.gt.0.)then
              call rect(a,b)
              write(6,200) a
c             ^^^^^^^^^^^^^^^ rect を call したあとには a の値が変更されている!!!
          else
              call circle(a)
          end if
c
          read(5,100) a,b
      end do
c
      close (5)
      close (6)
      stop
c
  100 format(2f10.0)
  200 format(' ','> call されたあとのメインプログラムの a の値は',
     &       1pe15.7,'になっている <')
      end
c -------------------------------------------------------------------
      subroutine circle(r)
c                            円
      common / const / pi
      write(6,100) r
      a=pi*r*r
      write(6,200) a
      return
c
  100 format(' ','radius =',1pe15.7)
  200 format(' ',15x,'A(circle) =',1pe15.7)
      end
c -------------------------------------------------------------------
      subroutine rect(a,b)
c                            四角
      write(6,100) a,b
      a=a*b
      write(6,200) a
      return
c
  100 format(' ','a =',1pe15.7,'  b =',e15.7)
  200 format(' ',15x,'A(rectangle) =',1pe15.7)
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program area7
      implicit real*8(a-h,o-z)
c                            倍精度（8バイト実数）
      character*4 title(19),aend
      common / const / pi
      data aend / 'end ' /
      pi=4.d0*atan(1.d0)
c
      open(5,file='inarea.dat')
      open(6,file='outarea.txt')
c                            タイトルで終了制御
      read(5,100) title
      do while(title(1).ne.aend)
          write(6,200) title
          read(5,300) a,b
          if(b.gt.0.)then
              call rect(a,b)
          else
              call circle(a)
          end if
          read(5,100) title
      end do
c
      close (5)
      close (6)
      stop
c
  100 format(19a4)
  200 format(//' ',19a4)
  300 format(2f10.0)
c                            ファイル`inarea.dat'の例
c                            ----+----1----+----2
c                            Problem No.1 --> title は全部 76 文字より長くする
c                                   1.2        0.
c                            Next Problem
c                                   3.6       .45
c                            Last Problem
c                                    .2        0.
c                            end_                 <-- 4 文字目はスペース!
c                            ----+----1----+----2
      end
c -------------------------------------------------------------------
      subroutine circle(r)
c                            円
      implicit real*8(a-h,o-z)
      common / const / pi
      write(6,100) r
      a=pi*r*r
      write(6,200) a
      return
c
  100 format(' ','radius =',1pd15.7)
  200 format(' ',15x,'A(circle) =',1pd15.7)
      end
c -------------------------------------------------------------------
      subroutine rect(a,b)
      implicit real*8(a-h,o-z)
c                            四角
      write(6,100) a,b
      a=a*b
      write(6,200) a
      return
c
  100 format(' ','a =',1pd15.7,'  b =',d15.7)
  200 format(' ',15x,'A(rectangle) =',1pd15.7)
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c input-data file `inarea.dat'
Problem No.1                                                                
       1.2        0.
Next Problem                                                                
       3.6       .45
Last Problem                                                                
        .2        0.
end                                                                         
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program inprod
c                            ex8: ２つのベクトルの内積を求める
      character*4 title(19),aend
      dimension a(3),b(3)
      data aend / 'end ' /
c
      open(5,file='invec.dat')
      open(6,file='outvec.txt')
      read(5,100) title
      do while(title(1).ne.aend)
          write(6,200) title
c
          read(5,300) (a(i), i=1,3)
          read(5,300) (b(i), i=1,3)
          write(6,400) (a(i), i=1,3)
          write(6,400) (b(i), i=1,3)
c                            内積の計算
          prod=0.
          do i=1,3
              prod=prod + a(i)*b(i)
          end do
          write(6,500) prod
c
          read(5,100) title
      end do
c
      close (5)
      close (6)
      stop
c
  100 format(19a4)
  200 format(//' ',19a4)
  300 format(3f10.0)
  400 format(' ','1:',1pe15.7,'  2:',e15.7,'  3:',e15.7)
  500 format(' ',15x,'Inner product =',1pe15.7)
c                            ファイル`invec.dat'の例
c                            ----+----1----+----2----+----3
c                            Problem No.1 (3-D) --> title は全部 76 文字より長くする
c                                   1.2       2.2       1.1
c                                   -.3       3.9      -4.5
c                            Next Problem (2-D)
c                                   2.3      -4.5        0.
c                                   8.1       .34        0.
c                            Last Problem (3-D)
c                                    .2       2.1       -.2
c                                  -1.1      -3.2       9.1
c                            end_                 <-- 4 文字目はスペース!
c                            ----+----1----+----2----+----3
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c input-data file `invec.dat'
Problem No.1 (3-D)                                                          
       1.2       2.2       1.1
       -.3       3.9      -4.5
Next Problem (2-D)                                                          
       2.3      -4.5        0.
       8.1       .34        0.
Last Problem (3-D)                                                          
        .2       2.1       -.2
      -1.1      -3.2       9.1
end                                                                         
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program param1
c                            ex9: データのやりとりと宣言について
      data n / 5 /
      write(6,100) n
      n = 2
      write(6,100) n
      data n / 8 /
      write(6,100) n
      stop
c
  100 format(' ',i5)
c                            このプログラムの出力結果は
c                                8
c                                2
c                                2
c                            です。なぜでしょう?
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program param2
c                            ex10: データのやりとりと宣言について
      parameter (n=5)
      write(6,100) n
      n = 2
      write(6,100) n
      stop
c
  100 format(' ',i5)
c                            このプログラムは
c                            第5行の`n=2'が文法違反です。
c                            なぜでしょう?
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program param3
c                            ex11: データのやりとりと宣言について
      parameter (n=5)
      common n
      write(6,100) n
      call suba
      stop
c
  100 format(' ',i5)
c                            このプログラムは
c                            第4行の`common n'が文法違反です。
c                            なぜでしょう?
      end
c -------------------------------------------------------------------
      subroutine suba
      common n
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program param4
c                            ex12: データのやりとりと宣言について
      common n
      data n / 5 /
      write(6,100) n
      call suba
      stop
c
  100 format(' ',i5)
c                            このプログラムも
c                            第3行の`common n'が文法違反です。
c                            なぜでしょう?
      end
c -------------------------------------------------------------------
      subroutine suba
      common k
      write(6,100) k
      return
c
  100 format(' ',i5)
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program param5
c                            ex13: データのやりとりと宣言について
      data n / 5 /
      write(6,100) n
      call suba(n)
      write(6,100) n
      stop
c
  100 format(' ',i5)
c                            このプログラムも subroutine suba の
c                            第2行の`common k'が文法違反です。
c                            なぜでしょう?
      end
c -------------------------------------------------------------------
      subroutine suba(k)
      common k
      k=2*k
      call subb
      return
      end
c -------------------------------------------------------------------
      subroutine subb
      common j
      j=j/2
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program param5
c                            ex13_1: データのやりとりと宣言について
      common n
      n=5
      write(6,100) n
      call suba(n)
      write(6,100) n
      stop
c
  100 format(' ',i5)
      end
c -------------------------------------------------------------------
      subroutine suba(k)
      k=2*k
      call subb
      return
      end
c -------------------------------------------------------------------
      subroutine subb
      common j
      j=j/2
      write(6,100) j
      return
  100 format(' ',i5)
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program inv22
c                            ex14: 2x2の逆行列
      character*4 title(19),aend
      dimension a(2,2),b(2,2),c(2,2)
      data aend / 'end ' /
c
      open(5,file='inmat.dat')
      open(6,file='outmat.txt')
c
      read(5,100) title
      do while(title(1).ne.aend)
          write(6,200) title
c
          do i=1,2
              read(5,300) (a(i,j), j=1,2)
              write(6,400) (a(i,j), j=1,2)
          end do
c
          det=a(1,1)*a(2,2)-a(1,2)*a(2,1)
          write(6,500) det
c
          b(1,1)=a(2,2)/det
          b(1,2)=-a(1,2)/det
          b(2,1)=-a(2,1)/det
          b(2,2)=a(1,1)/det
c
          do i=1,2
              write(6,400) (b(i,j), j=1,2)
          end do
c                            チェック
          do i=1,2
              do j=1,2
                  c(i,j)=0.
                  do k=1,2
                      c(i,j)=c(i,j)+a(i,k)*b(k,j)
                  end do
              end do
          end do
c                             きちんと単位行列になるか?
c                             どのくらいおかしいか?　なぜか?
          write(6,600)
          do i=1,2
              write(6,400) (c(i,j), j=1,2)
          end do
c
          read(5,100) title
      end do
c
      close (5)
      close (6)
      stop
c
  100 format(19a4)
  200 format(//' ',19a4/' ','Given matrix:')
  300 format(2f10.0)
  400 format(' ',' |',1pe15.7,'  |',e15.7,' |')
  500 format(' ',15x,'Determinant =',1pe15.7/' ','Inverse matrix:')
  600 format(' ','Check.....')
c                            ファイル`inmat.dat'の例
c                            ----+----1----+----2
c                            Problem No.1 --> title は全部 76 文字より長くする
c                                   1.2       2.2
c                                   -.3       3.9
c                            Next Problem
c                                   2.3      -4.5
c                                   8.1       .34
c                            Last Problem
c                                    .2       2.1
c                                  -1.1      -3.2
c                            end_                 <-- 4 文字目はスペース!
c                            ----+----1----+----2
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c input-data file `inmat.dat'
Problem No.1                                                                
       1.2       2.2
       -.3       3.9
Next Problem                                                                
       2.3      -4.5
       8.1       .34
Last Problem                                                                
        .2       2.1
      -1.1      -3.2
end                                                                         
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program matprd
c                            ex15: 行列の積 {d} = [A(T)] [B] [A] {c}
      implicit real*8 (a-h,o-z)
      character*4 title(19),aend
      dimension a(8,8),b(8,8),c(8),d(8),e(8,8)
      data aend / 'end ' /
c
      open(5,file='inmatp.dat')
      open(6,file='outmatp.txt')
c
      read(5,100) title
      do while(title(1).ne.aend)
          write(6,200) title
c
          read(5,600) mat
          if(mat.le.8.and.mat.gt.0)then
c
              write(6,500) mat
              do i=1,mat
                  read(5,300) (a(i,j), j=1,mat)
                  write(6,400) (a(i,j), j=1,mat)
              end do
              write(6,350)
              do i=1,mat
                  read(5,300) (b(i,j), j=1,mat)
                  write(6,400) (b(i,j), j=1,mat)
              end do
              write(6,450)
              read(5,300) (c(i), i=1,mat)
              write(6,400) (c(i), i=1,mat)
c
              do i=1,mat
                  d(i)=0.d0
                  do j=1,mat
                     d(i)=d(i)+a(i,j)*c(j)
                     e(i,j)=0.d0
                     do k=1,mat
                        e(i,j)=e(i,j)+a(k,i)*b(k,j)
                     end do
                  end do
              end do
c
              do i=1,mat
                  c(i)=0.d0
                  do j=1,mat
                      c(i)=c(i)+e(i,j)*d(j)
                  end do
              end do
c
              write(6,550)
              write(6,400) (c(i), i=1,mat)
c
          else
              write(6,700)
              close (5)
              close (6)
              stop
          end if
c
          read(5,100) title
      end do
c
      close (5)
      close (6)
      stop
c
  100 format(19a4)
  200 format(//' ',19a4)
  300 format(8f10.0)
  350 format(' ',' [B]')
  400 format(' ',8(' |',1pd15.7),' |')
  450 format(' ',' {c}')
  500 format(' ','Size of matrices =',i5/' ','Given matrices:'/
     &       ' ',' [A]')
  550 format(' ','Obtained vector:'/' ',' {d}')
  600 format(i5)
  700 format(' ','***** Data error : 0<mat<9 *****')
c                            ファイル`inmatp.dat'の例
c                            ----+----1----+----2----+----3
c                            Problem No.1 --> title は全部 76 文字より長くする
c                                2
c                                   1.2       2.2
c                                   -.3       3.9
c                                   2.3      -4.5
c                                   8.1       .34
c                                   .55       -.2
c                            Next Problem
c                                3
c                                    .2       2.1       -.4
c                                  -1.1      -3.2       2.7
c                                   1.2       2.2       -.9
c                                   -.3       3.9       5.7
c                                   7.1       8.8       -.8
c                                   3.2      -1.2        1.
c                                    1.       -.5        2.
c                            end_                 <-- 4 文字目はスペース!
c                            ----+----1----+----2----+----3
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c input-data file `inmatp.dat'
Problem No.1                                                                
    2
       1.2       2.2
       -.3       3.9
       2.3      -4.5
       8.1       .34
       .55       -.2
Next Problem                                                                
    3
        .2       2.1       -.4
      -1.1      -3.2       2.7
       1.2       2.2       -.9
       -.3       3.9       5.7
       7.1       8.8       -.8
       3.2      -1.2        1.
        1.       -.5        2.
end                                                                         
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program matprd
      implicit real*8 (a-h,o-z)
      character*4 title(19),aend
      dimension a(8,8),b(8,8),c(8),d(8),e(8,8)
      data aend / 'end ' /
c
      open(5,file='inmatp.dat')
      open(6,file='outmatp.txt')
c
      read(5,100) title
      do while(title(1).ne.aend)
          write(6,200) title
c
          read(5,600) mat
          if(mat.le.8.and.mat.gt.0)then
              write(6,500) mat
              do i=1,mat
                  read(5,300) (a(i,j), j=1,mat)
                  write(6,400) (a(i,j), j=1,mat)
              end do
              write(6,350)
              do i=1,mat
                  read(5,300) (b(i,j), j=1,mat)
                  write(6,400) (b(i,j), j=1,mat)
              end do
              write(6,450)
              read(5,300) (c(i), i=1,mat)
              write(6,400) (c(i), i=1,mat)
c
              call prod(a,b,c,d,e,mat)
c
              write(6,550)
              write(6,400) (c(i), i=1,mat)
          else
              write(6,700)
              close (5)
              close (6)
              stop
          end if
c
          read(5,100) title
      end do
c
      close (5)
      close (6)
      stop
c
  100 format(19a4)
  200 format(//' ',19a4)
  300 format(8f10.0)
  350 format(' ',' [B]')
  400 format(' ',8(' |',1pd15.7),' |')
  450 format(' ',' {c}')
  500 format(' ','Size of matrices =',i5/' ','Given matrices:'/
     &       ' ',' [A]')
  550 format(' ','Obtained vector:'/' ',' {d}')
  600 format(i5)
  700 format(' ','***** Data error : 0<mat<9 *****')
      end
c -------------------------------------------------------------------
      subroutine prod(a,b,c,d,e,mat)
      implicit real*8 (a-h,o-z)
      dimension a(8,1),b(8,1),c(1),d(1),e(8,1)
c
      do i=1,mat
          d(i)=0.d0
          do j=1,mat
              d(i)=d(i)+a(i,j)*c(j)
              e(i,j)=0.d0
              do k=1,mat
                  e(i,j)=e(i,j)+a(k,i)*b(k,j)
              end do
          end do
      end do
      do i=1,mat
          c(i)=0.d0
          do j=1,mat
              c(i)=c(i)+e(i,j)*d(j)
          end do
      end do
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program matprd
      implicit real*8 (a-h,o-z)
      character*4 title(19),aend
      dimension a(8,8),b(8,8),c(8),d(8),e(8,8)
      data aend / 'end ' /
c
      open(5,file='inmatp.dat')
      open(6,file='outmatp.txt')
c
      read(5,100) title
      do while(title(1).ne.aend)
          write(6,200) title
c
          read(5,600) mat
          if(mat.le.8.and.mat.gt.0)then
              write(6,500) mat
c
              call matin(a,b,c,mat)
c
              call prod(a,b,c,d,e,mat)
c
              call matout(c,mat)
          else
              write(6,700)
              close (5)
              close (6)
              stop
          end if
c
          read(5,100) title
      end do
c
      close (5)
      close (6)
      stop
c
  100 format(19a4)
  200 format(//' ',19a4)
  500 format(' ','Size of matrices =',i5/' ','Given matrices:'/
     &       ' ',' [A]')
  600 format(i5)
  700 format(' ','***** Data error : 0<mat<9 *****')
      end
c -------------------------------------------------------------------
      subroutine matin(a,b,c,mat)
      implicit real*8 (a-h,o-z)
      dimension a(mat,1),b(mat,1),c(1)
c
      do i=1,mat
          read(5,300) (a(i,j), j=1,mat)
          write(6,400) (a(i,j), j=1,mat)
      end do
      write(6,350)
      do i=1,mat
          read(5,300) (b(i,j), j=1,mat)
          write(6,400) (b(i,j), j=1,mat)
      end do
      write(6,450)
      read(5,300) (c(i), i=1,mat)
      write(6,400) (c(i), i=1,mat)
      return
  300 format(8f10.0)
  350 format(' ',' [B]')
  400 format(' ',8(' |',1pd15.7),' |')
  450 format(' ',' {c}')
      end
c -------------------------------------------------------------------
      subroutine prod(a,b,c,d,e,mat)
      implicit real*8 (a-h,o-z)
      dimension a(mat,1),b(mat,1),c(mat),d(mat),e(mat,1)
c
      do i=1,mat
          d(i)=0.d0
          do j=1,mat
              d(i)=d(i)+a(i,j)*c(j)
              e(i,j)=0.d0
              do k=1,mat
                  e(i,j)=e(i,j)+a(k,i)*b(k,j)
              end do
          end do
      end do
      do i=1,mat
          c(i)=0.d0
          do j=1,mat
              c(i)=c(i)+e(i,j)*d(j)
          end do
      end do
      return
      end
c -------------------------------------------------------------------
      subroutine matout(c,mat)
      implicit real*8 (a-h,o-z)
      dimension c(mat)
      write(6,550)
      write(6,400) (c(i), i=1,mat)
      return
  400 format(' ',8(' |',1pd15.7),' |')
  550 format(' ','Obtained vector:'/' ',' {d}')
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program matprd
      implicit real*8 (a-h,o-z)
      character*4 title(19),aend
      common / dtarea / a(8000)
      data aend / 'end ' /
c                            データ用メモリの最大 = 8000 とした
      lasta=8000
      open(5,file='inmatp.dat')
      open(6,file='outmatp.txt')
c
      read(5,100) title
      do while(title(1).ne.aend)
          write(6,200) title
c
          read(5,300) mat
          mmat=mat*mat
c                            各配列の番地を設定する
          addrs1=1
          addrs2=addrs1+mmat
          addrs3=addrs2+mmat
          addrs4=addrs3+mat
          addrs5=addrs4+mat
          lastad=addrs5+mmat-1
c                                      最後の番地が許容できるか?
          if(lastad.gt.lasta)then
              write(6,400)
              close (5)
              close (6)
              stop
          end if
c                                                  本当のメインへ
          call matcal(a(addrs1),a(addrs2),a(addrs3),
     &                                    a(addrs4),a(addrs5),mat)
c
          read(5,100) title
      end do
c
      close (5)
      close (6)
      stop
c
  100 format(19a4)
  200 format(//' ',19a4)
  300 format(i5)
  400 format(' ','***** Data error : mat is too large! *****')
c
      end
c -------------------------------------------------------------------
      subroutine matcal(a,b,c,d,e,mat)
      implicit real*8 (a-h,o-z)
      dimension a(mat,1),b(mat,1),c(mat),d(mat),e(mat,1)
c
      write(6,100) mat
      do i=1,mat
          read(5,200) (a(i,j), j=1,mat)
          write(6,300) (a(i,j), j=1,mat)
      end do
      write(6,400)
      do i=1,mat
          read(5,200) (b(i,j), j=1,mat)
          write(6,300) (b(i,j), j=1,mat)
      end do
      write(6,500)
      read(5,200) (c(i), i=1,mat)
      write(6,300) (c(i), i=1,mat)
c
      do i=1,mat
          d(i)=0.d0
          do j=1,mat
              d(i)=d(i)+a(i,j)*c(j)
              e(i,j)=0.d0
              do k=1,mat
                  e(i,j)=e(i,j)+a(k,i)*b(k,j)
              end do
          end do
      end do
c
      do i=1,mat
          c(i)=0.d0
          do j=1,mat
              c(i)=c(i)+e(i,j)*d(j)
          end do
      end do
c
      write(6,600)
      write(6,300) (c(i), i=1,mat)
c
      return
c
  100 format(' ','Size of matrices =',i5/' ','Given matrices:'/
     &       ' ',' [A]')
  200 format(8f10.0)
  300 format(' ',8(' |',1pd15.7),' |')
  400 format(' ',' [B]')
  500 format(' ',' {c}')
  600 format(' ','Obtained vector:'/' ',' {d}')
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program c_arry
c                            ex16_1: 整合配列
      implicit real*8 (a-h,o-z)
      parameter (dtmaxa=8000)
      common / dtarea / a(dtmaxa)
c
      read(5,100) n,m
c
      n1=1
      n2=n1+n
      n3=n2+m
      n4=n3+n*n
      n5=n4+m*n
      n6=n5+n
      nlast=n6+n*m-1
      if(nlast.gt.dtmaxa) stop
c
      call sub1(a(n3),a(n1),a(n4),a(n5),a(n6),n,m)
c
      n6=n5+m*m
      nlast=n6+m-1
      if(nlast.gt.dtmaxa) stop
c
      call sub2(a(n3),a(n4),a(n2),a(n5),a(n6),n,m)
c
      nlast=n3+m*m*n-1
      if(nlast.gt.dgmaxa) stop
c
      call sub3(a(n1),a(n2),a(n3),n,m)
c
      stop
  100 format(2i5)
      end
c -------------------------------------------------------------------
      subroutine sub1(a,b,c,e,f,n,m)
      implicit real*8 (a-h,o-z)
      dimension a(n,n),b(n),c(m,n),e(n),f(n,m)
c
      return
      end
c -------------------------------------------------------------------
      subroutine sub2(a,c,d,e,f,n,m)
      implicit real*8 (a-h,o-z)
      dimension a(n,n),c(m,n),d(m),e(m,m),f(m)
c
      return
      end
c -------------------------------------------------------------------
      subroutine sub3(b,d,e,n,m)
      implicit real*8 (a-h,o-z)
      dimension b(n),d(m),e(m,m,n)
c
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program l_regr
c            ex17: ソートと最小2乗法
      implicit real*8 (a-h,o-z)
      character*4 title(19),aend
      common / dtarea / a(8000)
      data aend / 'end ' /
c
      lasta=8000
      open(5,file='regrin.dat')
      open(6,file='regrout.txt')
c                                  グラフィックス用ファイル
      open(1,file='regrph.dat')
c
      read(5,100) title
      do while(title(1).ne.aend)
          write(6,200) title
c                                  タイトルもグラフィックス用ファイルに
          write(1,500) title
          read(5,300) nn
          addrs1=1
          addrs2=addrs1+nn
          lastad=addrs2+nn-1
          if(lastad.gt.lasta)then
              write(6,400)
              close (1)
              close (5)
              close (6)
              stop
          end if
c
          call regr(a(addrs1),a(addrs2),nn)
          read(5,100) title
      end do
c
      close (1)
      close (5)
      close (6)
      stop
c
  100 format(19a4)
  200 format(//' ',19a4)
  300 format(i5)
  400 format(' ','***** Data error : nn is too large! *****')
  500 format(' ',19a4)
c
      end
c
      subroutine regr(x,y,nn)
      implicit real*8 (a-h,o-z)
      dimension x(nn),y(nn)
c
      write(6,100) nn
      do i=1,nn
          read(5,200) x(i),y(i)
      end do
c                            小さい順に並べる
      do i=1,nn-1
          do j=i+1,nn
              if(x(i).gt.x(j))then
                  xy=x(i)
                  x(i)=x(j)
                  x(j)=xy
                  xy=y(i)
                  y(i)=y(j)
                  y(j)=xy
              end if
          end do
      end do
c                            グラフィックスデータ用出力も
      do i=1,nn
          write(6,300) x(i),y(i)
          write(1,600) x(i),y(i)
      end do
c                            最小2乗法の係数計算
      sumxx=0.d0
      sumxy=0.d0
      sumx =0.d0
      sumy =0.d0
      do i=1,nn
          sumxx=sumxx+x(i)*x(i)
          sumxy=sumxy+x(i)*y(i)
          sumx =sumx +x(i)
          sumy =sumy +y(i)
      end do
c
      a=sumxy/sumxx
      write(6,400) a
c
      det=sumxx*float(nn)-sumx*sumx
      a=(float(nn)*sumxy-sumx*sumy)/det
      b=(-sumx*sumxy+sumxx*sumy)/det
      write(6,500) a,b
c
      return
c
  100 format(' ','Number of data =',i5/' ','Sorted data (x,y):')
  200 format(2f10.0)
  300 format(' ',5x,'x=',1pd15.7,': y=',d15.7)
  400 format(' ',2x,'Linear regression : y = (',1pd15.7,' ) * x')
  500 format(' ',2x,'Linear regression : y = (',1pd15.7,' ) * x + (',
     &           d15.7,' )')
c                            グラフィックス用は例えばE形式出力
  600 format(' ',1p2e15.7)
c
c     data example
c     ----+----1----+----2
c     Data No.1 --> title は全部 76 文字より長くする
c        11
c         23.186    23.754
c         81.354    81.192
c         73.363    74.974
c         44.265    45.186
c         61.761    61.290
c         33.535    34.550
c         16.067    16.713
c         60.598    59.715
c         79.709    80.166
c         90.358    92.894
c         53.447    54.974
c     end_                 <-- 4 文字目はスペース!
c     ----+----1----+----2
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c input-data file `regrin.dat'
Data No.1                                                                   
   11
    23.186    23.754
    81.354    81.192
    73.363    74.974
    44.265    45.186
    61.761    61.290
    33.535    34.550
    16.067    16.713
    60.598    59.715
    79.709    80.166
    90.358    92.894
    53.447    54.974
end                                                                         
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program bisect
c                            ex18: f(x)=tan(x)-x=0 の解は?
      implicit real*8 (a-h,o-z)
c                                  キーボード入力・画面出力
      do
c
      write(*,510)
      read(*,*) eps
      if(eps.le.0.d0.or.eps.ge.1.d0) stop
      write(*,520) eps
c
      diff=1.0d1*eps
      do while(diff.gt.eps)
          write(*,530)
          read(*,*) xl
          write(*,540) xl
          write(*,550)
          read(*,*) xh
          write(*,560) xh
          dx=(xh-xl)/1.d1
c
          do i=0,10
              x=xl+float(i)*dx
              write(*,600) x,f(x)
          end do
c
          fl=f(xl)
          fh=f(xh)
          if(fl*fh.le.0.d0)then
              do while(diff.gt.eps)
                  x=(xl+xh)/2.d0
                  fc=f(x)
                  diff=abs((xh-xl)/x)
                  if(fc*fl.le.0.d0)then
                      xh=x
                  else
                      xl=x
                      fl=fc
                  end if
              end do
          else
              write(*,570)
          end if
c
      end do
      write(*,590) x,fc
c
      end do
c
  510 format(//' ','*** Bisection method to solve f(x)=0 ***'/
     &         ' ',' Tolerance (Enter 0.0 to stop) ?')
  520 format(' ',5x,'Tolerance =',1pd15.7)
  530 format(' ',' Lower  trial value ?')
  540 format(' ',5x,'Lower  trial value =',1pd15.7)
  550 format(' ',' Higher trial value ?')
  560 format(' ',5x,'Higher trial value =',1pd15.7)
  570 format(/' ','Both of f(xh) and f(xl) have the same sign.',
     &            'Try again !!!'/)
  590 format(' ',' >>> Possible solution =',1pd15.7,' <<<',
     &       5x,'(f(x)=',d15.7,')')
  600 format(' ',10x,'x=',1pd15.7,':  f(x)=',d15.7)
      end
c -------------------------------------------------------------------
      real*8 function f(x)
      implicit real*8 (a-h,o-z)
      f=tan(x)-x
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program bisect
      implicit real*8 (a-h,o-z)
      f(z)=tan(z)-z
c                            これ以下がメインプログラム
      do
      .....
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program nwrphs
c                            ex19: f(x)=tan(x)-x=0 の解は?
c                                  by Newton-Raphson method
      implicit real*8 (a-h,o-z)
      nmax=20
c
      do
c
      write(*,510)
      read(*,*) eps
      if(eps.le.0.d0.or.eps.ge.1.d0) stop
      write(*,520) eps
      error=1.d1*eps
c
      do while(error.gt.eps)
          iter=0
          write(*,530)
          read(*,*) x
          write(*,540) x
c
          do while(error.gt.eps.and.iter.le.nmax)
              iter=iter+1
              dx=-f(x)/df(x)
              x =x+dx
              error=abs(dx/x)
              write(*,570) iter,x,f(x),error
          end do
c
          if(error.le.eps)then
              write(*,550) x,f(x)
          else
              write(*,560) nmax
          end if
      end do
c
      end do
c
  510 format(//' ','*** Newton-Raphson method to solve f(x)=0 ***'/
     &         ' ',' Tolerance (Enter 0.0 to stop) ?')
  520 format(' ','Tolerance =',1pd15.7)
  530 format(' ',' First trial value for x ?')
  540 format(' ','First trial value =',1pd15.7)
  550 format(' ',' >>> Possible solution =',1pd15.7,' <<<',
     &       5x,'(f(x)=',d15.7,')')
  560 format(' ','No convergence within iteration ',i5,
     &           ' !!! Try again !!!'/)
  570 format(' ',' iter=',i5,':  x=',1pd15.7,':  f(x)=',d15.7,
     &           ':  error=',d15.7)
      end
c -------------------------------------------------------------------
      real*8 function f(x)
      implicit real*8 (a-h,o-z)
      f=tan(x)-x
      return
      end
c -------------------------------------------------------------------
      real*8 function df(x)
      implicit real*8 (a-h,o-z)
      df=1.d0/cos(x)/cos(x)-1.d0
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program nwrph2
c                            ex20: 連立 f(x)=0 の解は?
c                                  by Newton-Raphson method
      implicit real*8 (a-h,o-z)
      dimension a(2,2),x(2),f(2),dx(2)
c
      nmax=20
c
      do
c
      write(*,510)
      read(*,*) eps
      if(eps.le.0.d0.or.eps.ge.1.d0) stop
      write(*,520) eps
      error=1.d1*eps
c
      do while(error.gt.eps)
          iter=0
          write(*,530)
          read(*,*) x(1),x(2)
          write(*,540) x(1),x(2)
c
          do while(error.gt.eps.and.iter.le.nmax)
              iter=iter+1
              f(1)=f1(x)
              f(2)=f2(x)
              call tangnt(a,x)
              do i=1,2
                  dx(i)=0.d0
                  do j=1,2
                      dx(i)=dx(i)-a(i,j)*f(j)
                  end do
              end do
              do i=1,2
                  x(i)=x(i)+dx(i)
              end do
              error=sqrt((dx(1)**2+dx(2)**2)/(x(1)**2+x(2)**2))
              write(*,570) iter,x(1),x(2),error
          end do
c
          if(error.le.eps)then
              write(*,550) x(1),x(2)
              write(*,580) f1(x),f2(x)
          else
              write(*,560) nmax
          end if
      end do
c
      end do
c
  510 format(//' ','*** Newton-Raphson method to solve f(x,y)=0 ***'/
     &         ' ',' Tolerance (Enter 0.0 to stop) ?')
  520 format(' ','Tolerance =',1pd15.7)
  530 format(' ',' First trial value for x & y ?')
  540 format(' ','First trial value (x,y) =',1p2d15.7)
  550 format(' ',' >>> Possible solution (x,y) =',1p2d15.7,' <<<')
  580 format(' ',' >>>  f1(x) and f2(x) @(x,y) =',1p2d15.7,' <<<')
  560 format(' ','No convergence within iteration ',i5,
     &           ' !!! Try again !!!'/)
  570 format(' ',' iter=',i5,':  x=',1pd15.7,':  y=',d15.7,
     &           ':  error=',d15.7)
      end
c -------------------------------------------------------------------
      subroutine tangnt(a,x)
      implicit real*8 (a-h,o-z)
      dimension a(2,2),x(2)
      a(1,1)=2.d0*x(1)-3.d0*x(2)
      a(1,2)=-3.d0*x(1)
      a(2,1)=-4.d0*x(1)*x(2)-5.d0
      a(2,2)=3.d0*x(2)*x(2)-2.d0*x(1)*x(1)
c
      det=a(1,1)*a(2,2)-a(1,2)*a(2,1)
      b=a(2,2)/det
      a(2,2)=a(1,1)/det
      a(1,1)=b
      a(1,2)=-a(1,2)/det
      a(2,1)=-a(2,1)/det
      return
      end
c -------------------------------------------------------------------
      real*8 function f1(x)
      implicit real*8 (a-h,o-z)
      dimension x(2)
      f1=x(1)*x(1)-3.d0*x(1)*x(2)+8.1d0
      return
      end
c -------------------------------------------------------------------
      real*8 function f2(x)
      implicit real*8 (a-h,o-z)
      dimension x(2)
      f2=x(2)*x(2)*x(2)-2.d0*x(1)*x(1)*x(2)-5.d0*x(1)-1.58d1
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      dimension what(1000),is(1000),neces(1000),sary(1000)
c
c       データ入力や最初の荷重ステップの直前までの準備をしておく
c
c       そして，初めての計算か，中断したところからの継続か
      read(5,100) icont
      open(2,file='temp.dat',form='unformatted')
c       icont=0: 初めて     =1: 中断継続
      if(icont.eq.1)then
c       継続の場合，前回の最終結果の必要なデータを読んで巻戻し
          read(2) what,is,neces,sary
          rewind 2
      else
c       初めての計算ならここでデータ what, is, neces, sary を初期化する
      end if
c
c       実際各荷重ステップの計算はここから開始
c       ここで各荷重ステップの長い計算をする
c       そして，結果が出たとする　　必要なデータだけとりあえず記憶する
c
      write(2) what,is,neces,sary
      rewind 2
c
c       まだ続けるなら・・・次のステップの計算をする
c
      close(2)
      stop
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c