2.10 ライブラリ関数の実装  (1999/05/03 初版)

2.10.1 ベクトル操作関数群  行列操作ライブラリの関数群を実装し、ライブラリを完成させましょう。以下では、 C++言語の関数機能である、 ・関数のデフォルト引数 ・インライン関数 ・関数多重定義(オーバーロード) についても触れます。  行列の計算ではベクトルと共に利用することが一般的です。まずはベクトルの型とそ の関数群を定義しておきましょう。以下ではヘッダファイルを matrix10.h 、関数定 義部をmatrix10.cpp とします。ヘッダファイルは、始めと終わりにプリプロセッサ命 令 #ifndef MATRIX10_H #define MATRIX10_H  ...... #endif でコードを囲んでいるものとします。まずベクトルに関する部分を見て行きます。
// ヘッダファイル matrix10.h の一部(前半)

#include <math.h>

const double NEARLY_ZERO = 1.E-20;   // ゼロと見なせる
const double ZERO_TOLERANCE = 1.E-8; // 許容誤差ゼロ
const int NUM_OF_ITEMS  = 5;      // 画面に表示する1行の要素数


typedef double *Vector, **Matrix;

// ベクトルの動的作成 **************************************************
Vector vect_new(int n);

/*  機能:次元数nのベクトルを作成
    引数:(入力) n ベクトルの次元

    戻り値:ベクトルを返す  */


// ベクトルの解放 ******************************************************
void vect_del( Vector x);

/*  機能:ベクトルを解放する
    引数:(入力) x ベクトル */


// ベクトル要素の読み込み **********************************************
void vect_read(Vector x, int n);

/*  機能:ベクトル要素を入力する
    引数:(入力)x ベクトル
                  n ベクトルの次元
          (出力)x ベクトルを上書きして返す */



// ベクトルの表示 ******************************************************
void vect_print(const Vector x, int n,
              const char *s = 0, int items = NUM_OF_ITEMS);

/*  機能:ベクトル要素を出力する
    引数:(入力)x  ベクトル
                  n ベクトルの次元
                  *s タイトル表示のための文字列ポインタ(省略可)
                  items 1行に表示する要素数(省略時は既定値が使われる)*/


//  ベクトルのコピー ***************************************************
void vect_copy(const Vector x, Vector y, int n);

/*  機能:ベクトルxをベクトルyにコピーする
    引数:(入力) x  ベクトル
                   n ベクトルの次元
          (出力) y ベクトル  */


// ベクトルの内積 ******************************************************
double vect_prod(const Vector x, const Vector y, int n);

/*  機能:ベクトルxとベクトルyの内積計算
    引数:(入力) x  ベクトル
                   y ベクトル
                   n ベクトルの次元
    戻り値  内積の値 */


// ベクトルのノルム ****************************************************
inline double vect_norm(Vector x, int n)
{ return sqrt(vect_prod(x, x, n)); }

/*  機能:ベクトルxのノルム(大きさ)を計算
    引数:(入力) x  ベクトル
                   n ベクトルの次元
    戻り値: ノルムの値     */


// ベクトルの単位ベクトル化 ********************************************
void vect_unit(Vector x, Vector y, int n, double zero = NEARLY_ZERO);

/*  機能:ベクトルxを単位ベクトル化して
    引数:(入力) x  ベクトル
                   n ベクトルの次元
                  zero ノルムが0と見なせる判定値(省略時は既定値が使われる)   
          (出力) yベクトルに上書きして返す */

ベクトルは1次元配列として表現できますから、その型を typedef double *Vector; で定義します。型 Vector はポインタ「double *」のことです。  次に関数の定義部を示します。
// matrix10.cpp の一部
// 行列操作ライブラリ 関数定義部
#include <iostream.h>
#include <iomanip.h>
#include <stdlib.h>
#include "matrix10.h" // 行列操作ライブラリのヘッダファイル

//  エラー関数
void error(const char *message)
{
    cout << message << endl;
    abort();
}

// ベクトルの動的作成
Vector vect_new(int size)
{
    Vector x = new double[size];

    return x;
}

// ベクトルの解放
void vect_del( Vector x)
{
    delete [] x;
}
//  ベクトルの読み込み
void vect_read(Vector x, int n)
{
    cout << n << " 次元ベクトルの要素を入力してください" << endl;     
    for(int i = 0; i < n; i++){
        cout << "第 " << (i+1) << " 成分: ";
        cin  >> x[i];
    }
}

// ベクトルの表示
void vect_print(const Vector x, int n, const char *s, int items)
{
    if(s) cout << s << endl;
    cout.setf(ios::scientific); // 科学表記法
    for(int i = 0; i < n; i++){
        cout << setw(15) << x[i];
        if( !( (i+1) % items ) ) cout << endl;
    }
    if( n % items ) cout << endl;
}

//  ベクトルのコピー
void vect_copy(const Vector x, Vector y, int n)
{
    for(int i = 0; i < n; i++)
        y[i] = x[i];
}

// ベクトルの内積
double vect_prod(const Vector x, const Vector y, int n)
{
    double prod = 0.0;
    for(int i = 0; i < n; i++)
        prod += x[i] * y[i];
    return prod;
}

// ベクトルの単位ベクトル化
void vect_unit(Vector x, Vector y, int n, double zero)
{
    double norm = vect_norm(x, n);
    if(norm <= zero)
        error("0ノルムのベクトルです。プログラムを終了します。");

    for(int i = 0; i < n; i++)
        y[i] = x[i]/norm;
}
 次の例はベクトル操作関数のテストプログラム lst02_16.cpp です。
// lst02_16.cpp
// 行列操作ライブラリーのテストプログラム
// ベクトルのテスト
#include <iostream.h>
#include "matrix10.h"

int main()
{
    int dim;
    cout << "ベクトルの次元数を入力してください: ";
    cin  >> dim;

    Vector x = vect_new(dim);
    Vector y = vect_new(dim);

    vect_read(x, dim);
    vect_read(y, dim);
    vect_print(x, dim, "ベクトルX");
    vect_print(y, dim, "ベクトルY");

    cout << "ベクトルXのノルム: " << vect_norm(x, dim) << endl;
    cout << "ベクトルYのノルム: " << vect_norm(y, dim) << endl;
    cout << "ベクトルXYの内積: " << vect_prod(x, y, dim) << endl;    

    vect_unit(x, x, dim);
    vect_print(x, dim, "ベクトルXを単位ベクトル化");

    vect_del(y);
    vect_del(x);

    return 0;
}

ソースプログラム lst02_16.cpp
[実行結果]
ベクトルの次元数を入力してください: 3 
3 次元ベクトルの要素を入力してください 
第 1 成分: 1 
第 2 成分: 2 
第 3 成分: 2 
3 次元ベクトルの要素を入力してください 
第 1 成分: -1 
第 2 成分: 2 
第 3 成分: -2 
ベクトルX 
   1.000000e+00   2.000000e+00   2.000000e+00     
ベクトルY 
  -1.000000e+00   2.000000e+00 -2.000000e+00      
ベクトルXのノルム: 3.000000e+00 
ベクトルYのノルム: 3.000000e+00 
ベクトルXYの内積: -1.000000e+00 
ベクトルXを単位ベクトル化 
   3.333333e-01   6.666667e-01   6.666667e-01
テストプログラムにあるように、ベクトルを作成するには Vector x = vect_new(10); で次元数10のベクトルxが生成されます。このベクトルを削除するには関数 vect_del(x); を呼び出します。  また、ベクトル操作関数で、引数で渡すベクトルを変更しないものには修飾子const が付いていることに注目しましょう。例えば、ベクトルxをベクトルyにコピーする関 数は void vect_copy(const Vector x, Vector y, int n); という関数プロトタイプになっています。この場合yはxの要素がコピーされて変更さ れますが、コピー元のxは変更されません。


2.10.2 デフォルト引数  C++では、関数に渡す引数にデフォルト値(既定値)を与えることができます。ヘッ ダファイルにあるベクトルの要素をディスプレイに表示する関数のプロトタイプ宣言、 // ベクトルの表示 void vect_print(const Vector x, int n, const char *s = 0, int items = NUM_OF_ITEMS); /* 機能:ベクトル要素を出力する 引数:(入力)x ベクトル n ベクトルの次元 *s タイトル表示のための文字列定数ポインタ(省略可) items 1行に表示する要素数(省略時は既定値が使われる)*/ では2つの引数にデフォルト値を与えています。デフォルト引数を持つ関数はそれらを 省略して呼び出すことができ、その場合は既定値が使われます。この場合次の3通りの 呼び出しが可能です。 vect_print(x, 10); vect_print(y, 8, "ベクトルyの要素"); vect_print(z, 4, "ベクトルzの要素", 4);  デフォルト引数は、パラメータリストの後ろ(右端)から順に与えなければなりませ ん。またデフォルト引数は2重に定義することはできません。上の例では、関数定義部 のパラメータリストにはデフォルト引数がないことに注目しましょう。  変数名を省略した関数プロトタイプ宣言は次のようになります。 void vect_print(const Vector, int, const char * = 0, int = NUM_OF_ITEMS); つまり、 データ型 = デフォルト値 で与えます。


2.10.3 インライン関数  関数の呼び出しでは、実行時間のオーバーヘッドを伴います。プログラムの中で、関 数を呼び出すたびに制御が関数コードへ渡されるからです。C++では関数呼び出しによ るオーバーヘッドを少なくするために、インライン関数が用意されています。  インライン関数は、関数定義の先頭(戻り型の前)にキーワード inline を付けま す。コンパイラは、その関数を呼び出しているすべての場所に関数のコピーを埋め込み ます。matrix10.h にあるベクトルのノルム(大きさ)を計算する関数 // ベクトルのノルム inline double vect_norm(Vector x, int n) { return sqrt(vect_prod(x, x, n)); } はインライン関数として定義しています。しかしインライン関数化されるかどうかはコ ンパイラ(の能力)に依存します。インライン関数は小さな関数に対して使います。ル ープを含むような複雑な関数(for ループなど)では、コンパイラは inline を無視 して、普通の関数としてコンパイルします。インライン関数を使うと実行時間は減りま すが、実行ファイルのサイズは大きくなります。 [注]C言語では、これと似た機能にマクロ展開があります。マクロは#defineプリプロ セッサ命令でマクロ式を定義し、コンパイルする前にプログラムの中にマクロ式を展開 します。しかしC++ではインライン関数が使われます。マクロ展開は、C++の関数受け 渡しメカニズムと違って型の安全性を保証しないからです。

2.10.4 行列操作関数群  関数ライブラリの残りの部分を完成させます。
// matrix10.h の後半部分

// 行列の動的作成 ******************************************************
Matrix mat_new(int m, int n);

/*  機能:m×nの行列を作成
    引数:(入力) m 行数
                   n 列数
    戻り値: 行列を返す */

// 行列の解放 **********************************************************
void mat_del( Matrix a);

/*  機能:m×nの行列を解放
    引数:(入力) a 行列  */


// 行列要素の読み込み **************************************************
void mat_read(Matrix a, int m, int n);

/*  機能:m×nの行列aの要素を読み込む
    引数:(入力) a 行列
                   m 行数
                   n 列数 */


// 行列要素の表示 ******************************************************
void mat_print( const Matrix a, int m, int n,
                 const char *s = 0, int items = NUM_OF_ITEMS);

/*  機能:m×nの行列要素を出力する
    引数:(入力) a 行列
                   m 行数
                   n 列数
                  *s タイトル表示のための文字列ポインタ(省略可)
                  items 1行に表示する要素数(省略時は既定値が使われる)*/ 


//  行列のコピー *******************************************************
void mat_copy( const Matrix a, Matrix b, int m, int n);

/*  機能:行列aを行列bにコピー
    引数:(入力)a 行列
                  m 行数
                  n 列数
          (出力)b 行列に上書きして返す */


// 行列とベクトルの積 **************************************************
void mat_vect( const Matrix a, const Vector x, Vector y,
                int m, int n);

/*  機能:行列a(m×n)とベクトルx(n次元縦ベクトル)の積
    引数:(入力)a 行列
                  x ベクトル
                  m 行数
                  n 列数
          (出力)y ベクトル(m次元縦ベクトル)に上書きして返す  */


// ベクトルと行列の積 **************************************************
void vect_mat( const Vector x, const Matrix a, Vector y,
                 int m, int n);

/*  機能:ベクトルx(m次元横ベクトル)と行列a(m×n)の積
    引数:(入力)x ベクトル
                  a 行列
                  m ベクトルの次元(行数)
                  n aの列数
          (出力)y ベクトル(n次元横ベクトル)に上書きして返す  */


// 行列と行列の積 ******************************************************
void mat_prod(const Matrix a, const Matrix b, Matrix c,   
                 int m, int n, int k);

/* 機能:m×nの行列aとn×kの行列bの積
    引数:(入力)a 行列
                  b 行列
                  m aの行数
                  n aの列数(bの行数)
                  k bの列数
          (出力)c 行列に上書きして返す  */
及び行列操作の関数定義部分は次のようになります。
// matrix10.cpp の後半部分
// 行列の動的作成
Matrix mat_new(int row, int col)
{
    Matrix a = new double *[row + 1]; // 行を作る
    for(int i = 0; i < row; i++)      // 列を作る
        a[i] = new double[col];
    a[row] = 0;

    return a;
}

//行列の解放
void mat_del(Matrix a)
{
    Matrix b = a;
    while(*b != 0)
        delete [] *b++;  // 列を解放
    delete [] a;         // 行を解放
}

// 行列要素の読み込み
void mat_read(Matrix a, int m, int n)
{
    cout << m << " 行" << n << " 列の行列要素を入力してください"     
    << endl;
    for(int i = 0; i < m; i++)
        for(int j = 0; j < n; j++){
            cout << (i+1) << " 行" << (j+1) << " 列: ";
            cin  >> a[i][j];
        }
}

// 行列の表示
void mat_print(const Matrix a, int m, int n,
                 const char *s, int items)
{
    if(s) cout << s << endl;
    for(int i = 0; i < m; i++)
        vect_print( a[i], n, 0, items);
}

//  行列のコピー
void mat_copy(const Matrix a, Matrix b, int m, int n)
{
    for(int i = 0; i < m; i++)
        vect_copy( a[i], b[i], n);
}

// 行列とベクトルの積
void mat_vect(const Matrix a, const Vector x, Vector y,
                int m, int n)
{
    for(int i = 0; i < m; i++){
        y[i] = 0.0;
        for(int j = 0; j < n; j++)
            y[i] += a[i][j] * x[j];
    }
}

// ベクトルと行列の積
void vect_mat(const Vector x, const Matrix a, Vector y,
                 int m, int n)
{
    for(int i = 0; i < n; i++){
        y[i] = 0.0;
        for(int j = 0; j < m; j++)
            y[i] += x[j] * a[j][i];
    }
}

// 行列と行列の積
void mat_prod(const Matrix a, const Matrix b, Matrix c,
                int m, int n, int k)
{
    for(int i = 0; i < m; i++)
        vect_mat( a[i], b, c[i], n, k);
}
 行列の表示関数 mat_print()、行列のコピー関数 mat_copy()及び行列と行列の積 を計算する関数 mat_prod() ではベクトル操作関数を使っています。例えば、行列と 行列の積は次のコードと同じです。 // m×n行列とn×k行列の積 void mat_prod(const Matrix a, const Matrix b, Matrix c, int m, int n, int k) { for(int i = 0; i < m; i++) for(int j = 0; j < k; j++){ c[i][j] = 0.0; for(int l = 0; l < n; l++) c[i][j] += a[i][l]*b[l][j]; } }  テストプログラム lst02_17.cpp を示します。
// lst02_17.cpp
// 行列操作ライブラリーのテストプログラム
// 行列、ベクトルの積
#include 
#include "matrix10.h"

int main()
{
    Vector x,y;
    Matrix a,b,c;

    // 領域確保
    x = vect_new(2);
    y = vect_new(3);
    a = mat_new(2,3);
    b = mat_new(3,2);
    c = mat_new(2,2);

    vect_read(x, 2);
    mat_read(a, 2, 3);
    vect_print(x, 2, "ベクトルX");
    mat_print(a, 2, 3, "行列A");
    vect_mat(x, a, y, 2, 3);
    vect_print(y, 3, "ベクトルY=XA");

    mat_vect(a, y, x, 2, 3);
    vect_print(x, 2, "ベクトルX=AY");    

    mat_read(b, 3, 2);
    mat_print(b, 3, 2, "行列B");

    mat_prod(a, b, c, 2, 3, 2);
    mat_print(c, 2, 2, "行列C=AB");

    // 領域の解放
    vect_del(y);
    vect_del(x);
    mat_del(a);
    mat_del(b);
    mat_del(c);

    return 0;
}

ソースプログラム lst02_17.cpp
[実行結果]
2 次元ベクトルの要素を入力してください
第 1 成分: 1 
第 2 成分: -1 
2 行3 列の行列要素を入力してください 
1 行1 列: 1 
1 行2 列: 2 
1 行3 列: 3 
2 行1 列: 4 
2 行2 列: 5 
2 行3 列: 6 
ベクトルX 
   1.000000e+00  -1.000000e+00 
行列A 
   1.000000e+00   2.000000e+00   3.000000e+00          
   4.000000e+00   5.000000e+00   6.000000e+00          
ベクトルY=XA 
  -3.000000e+00  -3.000000e+00  -3.000000e+00          
ベクトルX=AY 
  -1.800000e+01  -4.500000e+01 
3 行2 列の行列要素を入力してください 
1 行1 列: -1 
1 行2 列: 0 
2 行1 列: 1 
2 行2 列: 0 
3 行1 列: 1 
3 行2 列: 2 
行列B 
  -1.000000e+00   0.000000e+00 
   1.000000e+00   0.000000e+00 
   1.000000e+00   2.000000e+00 
行列C=AB 
   4.000000e+00   6.000000e+00 
   7.000000e+00   1.200000e+01
【演習問題2.13】連立一次方程式を解く関数をつくりなさい。 以上のプログラムコードはこちらです。 ヘッダファイル: matrix10.h 関数定義プログラム:matrix10.cpp このプログラムには、ガウスの消去法を使った連立一次方程式を解く関数を含んでいま す。 int gauss(Matrix a, Vector x, int n, double eps = ZERO_TOLERANCE); /* 機能:n元連立1次方程式をガウスの消去法で解く 引数:(入力)a 係数行列(正方行列) x 定数項ベクトル n 次元数 eps 解の判定基準値(省略可) ピボットの絶対値がこれより小さいと解無しと判定 (出力)x 定数項ベクトルに上書きして返す 戻り値:1(成功)、0(解は不定または不能) */


2.10.5 関数の多重定義  C++では、引数の型が違う同名の関数をいくつも定義することができます。これを 関数多重定義(オーバーロード)といいます。  ここで作成した行列を表現するデータ Matrix a; はポインタへのポインタ「double **」として実現しています。これは次のような2 次元配列のデータ double b[5][10]; とは互換性がありません。ですからこの様な2次元配列をこれらの関数に渡すことは できません。なぜならここで作成した関数群は行列は「double **」型であるのに対 して、2次元配列bは「double [][10]」型(あるいは「double (*)[10]」型)で あるために型が一致しないからです。  一方ベクトル Vector の方は、単なるポインタ「double *」なので、1次元配列 をそのまま関数に渡すことができます。  しかし、2次元配列を使いたい場合もあります。次のプログラムは連立方程式 | 1 -2 3 1 2 | |x| | 13 | | 2 1 -1 2 4 | |y| | 7 | | 3 -1 -2 1 -1 | |z| = | -5 | | 1 3 1 -4 -2 | |u| | 2 | | 4 -2 1 3 -3 | |v| | -7 |, の解(x,y,z,u,v)を求めます。
// lst02_18.cpp
// 連立方程式の解

#include <iostream.h>
#include <iomanip.h>
#include "matrix10.h"

const int DIM = 5; // 行列の次元

// ユーザ定義関数
void mat_copy(const double a[][DIM], Matrix b,
        int row, int col)
{
    for(int i = 0; i < row; i++)
        for(int j = 0; j < col; j++)
            b[i][j] = a[i][j];
}

int main()
{
    int i,j;
    double a[][DIM] = { { 1.0,-2.0, 3.0, 1.0, 2.0},
                         { 2.0, 1.0,-1.0, 2.0, 4.0},
                         { 3.0,-1.0,-2.0, 1.0,-1.0},
                         { 1.0, 3.0, 1.0,-4.0,-2.0},
                         { 4.0,-2.0, 1.0, 3.0,-3.0}};
    double b[] = { 13.0, 7.0, -5.0, 2.0, -7};

    cout << "係数行列と定数項ベクトル" << endl;
    for(i = 0; i < DIM; i++){
        cout << "|";
        for(j = 0; j < DIM; j++)
            cout << setw(5) << a[i][j];
        cout << " |  |" << setw(5) << b[i] << " |" << endl;
    }

    Matrix mat = mat_new(DIM, DIM);
    Vector vec = vect_new(DIM);

    // 読み込み
    mat_copy(a, mat, DIM, DIM);
    vect_copy(b, vec, DIM);

    int status = gauss( mat, vec, DIM);

    if(status)
        vect_print(vec, DIM, "連立方程式の解");
    else cout << "解けませんでした。" << endl;

    mat_del(mat);
    vect_del(vec);

    return 0;
}

ソースプログラム lst02_18.cpp
[実行結果]
係数行列と定数項ベクトル
|    1   -2    3    1    2 |  |   13 | 
|    2    1   -1    2    4 |  |    7 | 
|    3   -1   -2    1   -1 |  |   -5 | 
|    1    3    1   -4   -2 |  |    2 | 
|    4   -2    1    3   -3 |  |   -7 | 
連立方程式の解 
   1.000000e+00  -1.000000e+00   2.000000e+00  -2.000000e+00   3.000000e+00  
この連立方程式を解くプログラムでは、はじめに5×5の係数行列を2次元配列とし て用意しておいて、これを行列 Matrix型に変換(コピー)してライブラリ関数を利 用します。  コピーするためのユーザ定義関数 void mat_copy(const double a[][DIM], Matrix b, int row, int col); を使って、2次元配列をMatrix型(つまり「double **」型)のデータにコピーしま す。  この関数は行列ライブラリ関数 void mat_copy( const Matrix a, Matrix b, int m, int n); と同じ名前です。コンパイラは同名の関数をパラメータリスト(引数の数、型、順序) の違いから判別して関数を区別します。この機能を使うことで、似通った仕事をする 複数の関数を同名の関数にして異なる働き(または呼び出し方法)を持たせることが できます。この様に関数の多重定義は、デフォルト引数と共に柔軟性のある関数を定 義することができます。  注意すべきことは、多重定義する関数のパラメータリスト(引数の数、型、順序) が同じにならないようにすることです。そのような関数をつくると構文エラーになり ます。特にデフォルト引数を持つ関数で、引数を省略した場合の関数名と衝突しない ようにしなければなりません。 (*) 上のユーザー定義の関数で、古いコンパイラにはこの様に const を使えないも のがあります。 その場合は、関数のパラメータリストの宣言からキーワード const を削除してください。

| 目次 | 前のページ | 次のページ |

Copyright(c) 1999 Yamada,K