精度保証数値計算レポート(演習1,2)

岩村誠

1 演習1

Matlab(Scilab,Octave)により連立一次方程式の精度保証をするプログラムを作れ。

1.1 環境

実行環境は以下の通り。

CPU Mobile Celeron 300MHz

OS Linux 2.2.14(Vine Linux 2.0)

数値計算ソフト GNU Octave, version 2.1.31

OctaveにはMatlibのsetround相当の関数がないため、まずそれと同等の関数を作った。以下のプログラムをsetround.ccという名前で保存し、 "mkoctfile setround.cc"にてsetround.octを作成した。

#include <fpu_control.h>
#include <octave/oct.h>
#include <octave/pager.h>
int _RoundUp    = _FPU_DEFAULT | _FPU_RC_UP;
int _RoundDown  = _FPU_DEFAULT | _FPU_RC_DOWN;
int _RoundNear  = _FPU_DEFAULT | _FPU_RC_NEAREST;
DEFUN_DLD (setround, args, , "The `setround'.")
{
    octave_value    retval = "OK";
    if(args.length() != 1){
        octave_stdout << "usage: setround([\"up\"|\"down\"|\"near\"])\n";
        retval = "NG";
    }else if(args(0).string_value() == "up"){
        asm("fldcw _RoundUp");
        octave_stdout << "setround ( up )\n";
    }else if(args(0).string_value() == "down"){
        octave_stdout << "setround ( down )\n";
        asm("fldcw _RoundDown");
    }else if(args(0).string_value() == "near"){
        asm("fldcw _RoundNear");
        octave_stdout << "setround ( near )\n";
    }else{
        octave_stdout << "usage: setround([\"up\"|\"down\"|\"near\"])\n";
        retval = "NG";
    }
    return retval;
}

1.2 プログラム

連立一次方程式の精度保証をするプログラムは以下のとおり。

n = 3
A = rand(n,n)
b = rand(n,1)
R = inv(A)
c = R * b
setround("down")
X1 = R * A - eye(n,n)
setround("up")
X2 = R * A - eye(n,n)
X = max(abs(X1),abs(X2))
D = norm(X,inf)
N = norm(R,inf)
r2 = A * c - b
setround("down")
r1 = A * c - b
k = 1 - D
r = max(abs(r1),abs(r2))
p = norm(r,inf)
setround("up")
e = (N * p) / k

1.3 実行結果

実行結果は以下のとおり。

n = 3
A =

  0.038796  0.613068  0.599751
  0.903270  0.823821  0.162296
  0.418286  0.911795  0.482365

b =

  0.055407
  0.012832
  0.792917

R =

   3.4901   3.5142  -5.5218
  -5.1472  -3.2487   7.4929
   6.7031   3.0935  -7.3021

c =

  -4.1399
   5.6144
  -5.3788

setround ( down )
ans = OK
X1 =

   -2.2204e-16   -9.9747e-18   -1.1970e-16
    8.4568e-18    0.0000e+00    2.8233e-16
   -1.3639e-16    2.4720e-16    0.0000e+00

setround ( up )
ans = OK
X2 =

   -1.1102e-16   -9.1073e-18   -1.1905e-16
    9.1073e-18    2.2204e-16    2.8298e-16
   -1.3574e-16    2.4850e-16    2.2204e-16

X =

   2.2204e-16   9.9747e-18   1.1970e-16
   9.1073e-18   2.2204e-16   2.8298e-16
   1.3639e-16   2.4850e-16   2.2204e-16

D =  6.0694e-16
N = 17.099
r2 =

   -1.5266e-16
   -1.7347e-18
   -2.2204e-16

setround ( down )
ans = OK
r1 =

   -1.5959e-16
   -3.4694e-18
   -3.3307e-16

k = 1.00000
r =

   1.5959e-16
   3.4694e-18
   3.3307e-16

p =  3.3307e-16
setround ( up )
ans = OK
e =  5.6950e-15

2 演習2

Octave,Scilabをインストール、動作確認をして感想を述べよ。

2.1 感想

GNU Octave version 2.1.31をVine Linux2.0にインストールすることにした。まずはoctave-2_1_31_tar.bz2を手に入れ展開、octave-2.1.31に移動し

./configure
make
su
make install

で簡単にインストールすることが出来た。

その後、最適化BLAS or ATLASで高速化出来るのを思いだし、最適化BLASを探す。しかしCeleron(L2Cache=128KB)のBLASが見つからないので、 ATLASをコンパイルすることにした。寝る前に仕込んでおいたら、さすがに朝にはコンパイル完了。出来上がったライブラリを/usr/local/libに放り込み、もう一度octaveを作り直した。その際、config.cache中に

ac_cv_lib_atlas_ATL_xerbla=${ac_cv_lib_atlas_ATL_xerbla='yes'}

という記述があるのを確認。勝手に認識された模様。再度makeし実行してみると、行列同士のかけ算は10倍程度速くなったが、 A \bは変わらなかった。

File translated from T_EX by T_TH, version 2.80 modified by Oishi
On 1 Dec 2000, 09:19.

URI: http://www.oishi.info.waseda.ac.jp/~oishi/accuracy/main.html

CPU	Mobile Celeron 300MHz
OS	Linux 2.2.14(Vine Linux 2.0)
数値計算ソフト	GNU Octave, version 2.1.31