東京理科大学 infoserv[更新日]2000.10.23


問13  次の Fortran プログラムの説明及びプログラムを読んで,設問 1 〜 4 に答えよ。

〔プログラムの説明〕

ガウスの消去法を用いて n 元連立 1 次方程式

を解くプログラムである。

(1)

n 元連立 1 次方程式の n(2≦n≦100)を変数 na に, 係数 aij(i=1, …, n; j=1, …, n),定数 bi(i=1, …, n)を 配列 a,b に入力する。

(2)

ガウスの消去法は,次に示す前進消去と後退代入の二つの処理を行うことによって, n 元連立 1 次方程式の解を求める。

 前進消去は,次に示す計算を k=1, …, n−1 について行う。ここで,akk をピボットと呼ぶ。

 後退代入は,次に示す計算を行い, n 元連立 1 次方程式の解 xi(i=1, …, n)を求める。

プログラムでは,内積 akj・xj を計算するために組込み関数 dot_product を 用いている。 dot_product (s,t) は,同じ大きさの実数型などの 1 次元配列 s,t を。

ベクトルとみなして,s と t の内積を返す。

(3) 図 1 に示す入力データに対する出力結果を図 2 に示す。

  3
  2.0   3.0  1.0  6.0
  1.0   4.0  2.0  4.0
  5.0  -1.0  6.0  3.0

図1 入力データ

  coefficients of simultaneous equations constant vector
  0.20000E+01 0.30000E+01 0.10000E+01 0.60000E+01
  0.10000E+01 0.40000E+01 0.20000E+01 0.40000E+01
  0.50000E+01 -0.10000E+01 0.60000E+01 0.30000E+01
solution =  0.20000E+01  0.10000E+01  -0.10000E+01

図2 出力結果

〔プログラム〕

(行番号)

01  real,parameter:: eps=1.0e-6
02  real,dimension(100,100):: a
03  real,dimension(100):: b,x
04  real:: piv,work,sum
05  integer:: na,row,column,k
06  read(*,*) na
07  read(*,*)(((a(row,column),column=1,na),b(row)),row=1,na)
08  write(*,'("       coefficients of simultaneous equations", &
   &            "   constant vector  ")')
09  do row=1,na
10    write(*,'("   ",101e14.5)')(a(row,column),column=1,na),b(row)
11  enddo
12  do k=1,na-1
13    piv = a(k,k)
14    if( abs(piv) < eps ) then
15      write(*,'(a,e12.5,a,e12.5)') &
   & ' The absolute value of pivot ',piv,' is less than ',eps
16      stop
17    endif
18    a(k,k:na) = a(k,k:na)/piv
19    b(k) = b(k)/piv
20    do row=k+1,na
21      work = a(row,k)
22      a(row,k+1:na) = a(row,k+1:na)-work*a(k,k+1:na)
23      b(row) = b(row)-work*b(k)
24    enddo
25  enddo
26  if( abs(a(na,na)) < eps ) then
27    write(*,'(a,e12.5,a,e12.5)') &
   & ' The absolute value of pivot ',a(na,na),' is less than ',eps
28    stop
29  endif
30  x(na) = b(na)/a(na,na)
31  do k=na-1,1,-1
32    sum = dot_product(a(k,k+1:na),x(k+1:na))
33    x(k) = b(k)-sum
34  enddo
35  write(*,'(/"  solution = ",100e14.5)')(x(row),row=1,na)
36  end
設問1   図 3 に示すデータを入力して実行したとき,行番号 22 の計算の結果,

The absolute value of pivot 0.00000E+00 is less than 0.10000E-05

というメッセージが表示され,プログラムが停止した。 絶対値が 1.0×10-6 より小さくなった配列要素とメッセージの出力場所として 正しい答えを,解答群の中から選べ。

   3
   2.0   2.0   4.0   6.0
   1.0   1.0   3.0   4.0
   3.0   5.0   6.0 5.0

図3 入力データ

解答群

ア a(1, 1) の絶対値が小さくなり,行番号 15 でメッセージを出力した。

イ a(1, 1) の絶対値が小さくなり,行番号 27 でメッセージを出力した。

ウ a(2, 2) の絶対値が小さくなり,行番号 15 でメッセージを出力した。

エ a(2, 2) の絶対値が小さくなり,行番号 27 でメッセージを出力した。

オ a(3, 3) の絶対値が小さくなり,行番号 15 でメッセージを出力した。

カ a(3, 3) の絶対値が小さくなり,行番号 27 でメッセージを出力した。

設問2   設問 1 で起こった問題を防ぐ方法としてピボット選択処理がある。 これは,方程式の順番を入れ替えても解は変わらないという性質に基づき, 絶対値が最も大きい値がピボットになるように,式の順序を入れ替える方式である。 ピボット選択処理を記述した次に示すプログラムの挿入場所として正しい答えを, 解答群の中から選べ。ここで,プログラムの行番号 4 と 5 の間に
real:: mpivot
integer:: mindex

を挿入する。

〔ピボット選択処理〕

mindex = k
mpivot = abs(a(k,k))
do row=k+1,na
  if( abs(a(row,k)) > mpivot ) then
    mindex = row
    mpivot = abs(a(row,k))
  endif
enddo
if( mindex /= k ) then
  do column=k,na
    work = a(k,column)
    a(k,column) = a(mindex,column)
    a(mindex,column) = work
  enddo
  work = b(k)
  b(k) = b(mindex)
  b(mindex) = work
endif

解答群

ア 行番号 12 と行番号 13 の間    イ 行番号 13 と行番号 14 の間

ウ 行番号 17 と行番号 18 の間    エ 行番号 19 と行番号 20 の間

設問3 図 3 に示す入力データを使い,設問 2 で示したピボット選択処理を行った場合, 1 回目のピボット選択処理が終了したときの配列 a の内容として正しい答えを, 解答群の中から選べ。

解答群

ア 4.0   2.0   2.0       イ 6.0   5.0   3.0
   3.0   1.0   1.0         4.0   2.0   2.0
   6.0   5.0   3.0         3.0   1.0   1.0

ウ 3.0   5.0   6.0    エ 3.0   5.0   6.0
   2.0   2.0   4.0         1.0   1.0   3.0
   1.0   1.0   3.0         2.0   2.0   4.0
設問4 プログラムの行番号 18,19 を削除し,行番号 21,33 を変更することによって, 同等の結果を得るプログラムを作成することができる。 行番号 21,33 の変更内容として正しい答えを,解答群の中から選べ。

解答群

ア 行番号 21 work = a(row,k)*piv

  行番号 33 x(k) = (b(k)-sum)/a(k,k)

イ 行番号 21 work = a(row,k)*piv

  行番号 33 x(k) = b(k)-sum/a(k,k)

ウ 行番号 21 work = a(row,k)/piv

  行番号 33 x(k) = (b(k)-sum)/a(k,k)

エ 行番号 21 work = a(row,k)/piv

  行番号 33 x(k) = b(k)-sum/a(k,k)


東京理科大学 infoserv 戻る 次頁:問14