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


岩村 誠

2001年02月05日

1  問題

Ax = bを精度保証するJavaプログラムを書け。

2  環境

利用した環境は以下の通り。

CPU Mobile Celeron 300MHz
OS Linux 2.2.14(Vine Linux 2.0)
JAVAパッケージ Jama Version 1.0.1
JDK IBM Developer Kit for Linux,Java Technology Edition,Version 1.1.8

3  答え

3.1  プログラム

行列演算にはJamaを用いた。 以下にAx = bを精度保証するJavaプログラム示す(このプログラムの名前をJNItest.javaとする)。 
import java.util.*;
import Jama.*;
public class JNItest
{
    public static void main (String[] args)
    {
        Matrix  A,b,R,c,E,X,X1,X2,r,r1,r2;
        double  D,N,k,p,e;
        int     n = 3;
        System.loadLibrary("setround");
        E = Matrix.identity(n,n);
        A = Matrix.random(n,n);
        b = Matrix.random(n,1);
        R = A.inverse();
        c = R.times(b);
        down();
        X1 = R.times(A).minus(E);
        up();                    
        X2 = R.times(A).minus(E);
        X = max(abs(X1),abs(X2));
        D = X.normInf();
        N = R.normInf();
        r2 = A.times(c).minus(b);
        down();
        r1 = A.times(c).minus(b);
        k = 1 - D;
        r = max(abs(r1),abs(r2));
        p = r.normInf();
        up();
        e = (N * p) / k;
        printMatrix("A",A);
        printMatrix("b",b);
        printMatrix("R",R);
        printMatrix("c",c);
        printMatrix("X1",X1);
        printMatrix("X2",X2);
        printMatrix("X",X);
        System.out.println("D=\n" + D + "\n");
        System.out.println("N=\n" + N + "\n");
        printMatrix("r1",r1);
        printMatrix("r2",r2);
        System.out.println("k=\n" + k + "\n");
        printMatrix("r",r);
        System.out.println("p=\n" + p + "\n");
        System.out.println("e=\n" + e + "\n");
    }
    public static Matrix abs(Matrix m)
    {
        int     row,col;                   
        double  elm;
        Matrix  ret;
        ret = m.copy();
        for(row = 0;row < ret.getRowDimension();row++){
            for(col = 0;col < ret.getColumnDimension();col++){
                if((elm = ret.get(row,col)) < 0){
                    ret.set(row,col,-elm);
                }
            }
        }    
        return ret;
    }
    public static Matrix max(Matrix m1,Matrix m2)
    {
        int     row,col;
        double  elm;
        Matrix  ret;
        ret = m1.copy();
        for(row = 0;row < ret.getRowDimension();row++){
            for(col = 0;col < ret.getColumnDimension();col++){
                if(ret.get(row,col) < (elm = m2.get(row,col))){
                    ret.set(row,col,elm);
                }
            }
        }
        return ret;
    }
    public static void printMatrix(String name,Matrix m)
    {
        int     row,col;
        System.out.println(name + "=");
        for(row = 0;row < m.getRowDimension();row++){
            for(col = 0;col < m.getColumnDimension();col++){
                System.out.print("" + m.get(row,col) + "\t");
            }
            System.out.println("");
        }
        System.out.println("");
    }                       
    public static native void up();
    public static native void down(); 
}
メンバ関数up(),down()はそれぞれ上への丸め、下への丸めを行う関数である。 JNIを用いて実装している。 以下はそのCによる実装である(次のプログラムの名前をJNItest.cとする)。 
#include    "JNItest.h"
#include    <fpu_control.h>
int _RoundUp    = _FPU_DEFAULT | _FPU_RC_UP;
int _RoundDown  = _FPU_DEFAULT | _FPU_RC_DOWN;
JNIEXPORT void JNICALL Java_JNItest_up(JNIEnv *env, jclass obj)
{
    asm volatile("fldcw _RoundUp");
}
 
JNIEXPORT void JNICALL Java_JNItest_down(JNIEnv *env, jclass obj)
{
    asm volatile("fldcw _RoundDown");
}
JNItest.hは(上記のJavaプログラムを JNItest.javaという名前で保存し、) 
javac JNItest.java
javah JNItest
で生成する。 最後に 
gcc -shared -I${JDKDIR}/include/ -I${JDKDIR}/include/linux/\
JNItest.c -o libsetround.so

としてShared Objectを生成しておく(Vine Linux 2.0の場合には，libsetround.soを${JDKDIR}/include/linux/i386/にコピーしておく)。

java JNItest

とすれば4.の実行結果が得られる。

注　ピンクの文字部分は大石による追補。

4  実行結果

実行結果は以下のとおり。 
A=
0.22159941403993644     0.10522301765697573     0.060393508599793755    
0.18360899847992895     0.2638448087126104      0.9763161167000047      
0.26051010362435856     0.2635488894146636      0.8587558384680678      

b=
0.4986128679241 
0.8051000126710501      
0.5059417590655145      

R=
-14.40476314619811      -34.89733015509862      40.687676848290806      
45.3137469791732        81.83186339851287       -96.22109051490244      
-9.536821742356382      -14.527516511267628     18.351445108776062      

c=
-14.69264641723548      
39.79458380239908       
-7.16652334740929       

X1=
-1.7763568394002505E-15 -1.9506965487359196E-15 -6.782768791069316E-15  
4.430483757644765E-15   4.218847493575595E-15   1.504352198367087E-14   
-7.125376677574735E-16  -6.782768791069316E-16  -2.3314683517128287E-15 

X2=
-1.6653345369377348E-15 -1.9480944635219544E-15 -6.772360450213455E-15  
4.435687928072696E-15   4.440892098500626E-15   1.5064338665382593E-14  
-7.116703060194851E-16  -6.774095173689432E-16  -2.220446049250313E-15  

X=
1.7763568394002505E-15  1.9506965487359196E-15  6.782768791069316E-15   
4.435687928072696E-15   4.440892098500626E-15   1.5064338665382593E-14  
7.125376677574735E-16   6.782768791069316E-16   2.3314683517128287E-15  

D=
2.3940918691955915E-14

N=
223.3667008925885

r1=
-7.216449660063518E-16  
-3.3306690738754696E-16 
-5.551115123125783E-16  

r2=
-6.661338147750939E-16  
-2.220446049250313E-16  
-4.440892098500626E-16  

k=
0.999999999999976

r=
7.216449660063518E-16   
3.3306690738754696E-16  
5.551115123125783E-16   

p=
7.216449660063518E-16

e=
1.6119145527258685E-13

5  最後に

Javaにはjava.math.BigDecimalというクラスが存在する。 これは任意精度のスケールなし整数値と、小数点以下の桁数を表す負でない 32bit整数のスケールから構成される。 このクラスには丸めの制御も実装されているため、 精度保証数値計算にも使える。 ただ任意精度という仕様を見る限りソフトウェアで、 各演算を実装していると考えられる。 またJavaは演算子のオーバロードをサポートしていないため、 JamaをこのBigDecimalに対応させるには、基本演算部分を全て該当する関数で 置き換える必要がある。



File translated from TEX by TTH, version 2.80. Modified by Oishi Shin'ichi
On 5 Feb 2001, 09:04.

©waseda university

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