簡易線形演算ライブラリ

線形演算ライブラリ

数値解析の授業でC言語を用いてプログラミングの演習を行っていますが、行列やベクトルを用いる必要があります。線形演算のライブラリには、Lapack、uBlas、Eigen等々がありますが、授業の中でインストールには時間をかけたくないですし、特に最近のオンライン授業では生徒への個別対応が難しいということから、様々な環境で確実に動くものが欲しいということで、簡単なものを自作してみました。

ただし、あくまで簡易的なものですので、機能は限定的ですし、計算効率も良くありません。本格的な数値計算には用いないようにして下さい。とはいえ、ちょっとした計算をお試しで実行してみるには便利ですので、紹介することにしました。

注意：2023/5/1に重大なバグを発見して修正しました。以前にダウンロードされた方は差し替えをお願いします。 行列ベクトル積を計算するときにメモリの初期化を忘れるという初歩的なミスです。PCの起動直後などでメモリが0クリアされている状態だと不具合が起きないので、長らく気が付きませんでした。

こちらにLinear.hというファイルがあります。普通にクリックするとソースが表示されてしまいます。右クリックして「リンク先を保存」等でダウンロードして下さい。

使い方は、Linear.h を cpp ファイルと同じ場所に置いて、プログラムの最初に #include "Linear.h" でインクルードするだけです。

以下、特徴です。

基本的な機能しか用いていないので、大抵のバージョンのC++で動くはずです。ただ、書き方が古すぎるということで警告は出るかもしれません。

すべてを1つのヘッダファイルに記述しているので、使いたいプログラムからインクルードするだけで済みます。C言語のファイル管理の慣習からは外れますが、使い易さを優先するため、この形にしました。

関数は値渡しで戻り値を返していますのでオーバーヘッドが大きいです。ここは効率化できないこともないですが、そもそも簡易的なものですので、効率よりシンプルさを優先としました。

具体的な機能

型（クラス）としては、ベクトルを表す Vector、整数値を取るベクトルを表す IntVector、行列を表す Matrixが使用できます。

    Vector b(100); // ベクトルを宣言
    IntVector p(30); // 整数型のベクトルを宣言
    Matrix A(10,10); // 行列を宣言

のように、行列やベクトルのサイズを指定して宣言することもできますし、

    Matrix A;
    A = B;
    Vector b = c;

のように、とりあえず空の状態で宣言して後から代入したり、宣言と同時に別の行列やベクトルで初期化したりしても構いません。

サイズを指定して宣言した場合は、中身は全て0で初期化されます。添え字の範囲は、0から指定した値までで、普通の配列のように使用することができます。例えば、Matrix A(20,30); と宣言した場合には A[0][0]～A[20][30] が使えます。添え字の範囲についてはチェックが入り、範囲を外れた添え字でアクセスするとエラーが出ます。

Vector および Matrix については、演算子オーバーロードにより、"+", "-", "+=", "-=", "*", "*=", "/", "/=" がオーバーロードされているため、自然な演算は自由に行えます。ただし、VectorとVectorの積は内積となります。また、Matrix同士の積、VectorとMatrixの積、Vector同士の積（内積）については、計算される添え字の範囲は1からになります。それ以外の演算については添え字は0から計算されます。

IntVector については、添え字を指定しての代入、IntVector同士の代入以外の演算はできません。

その他の機能については以下を見て下さい。

    int n = 5;

    Vector b(n); // ベクトルを宣言、b[0]～b[n]を使用可能。
    b.Shift(2);  // ベクトルの中身を添え字2つ分右にシフトする。添え字の範囲は0から。
                 // b[0]=8; b[1]=7; b[2]=6; b[3]=5; b[4]=4; b[5]=3; のときには c=b.Shift(2)は
                 // c[0]=0; c[1]=0; c[2]=8; c[3]=7; c[4]=6; c[5]=5; となる。
    b.Shift(-1); // ベクトルの中身を添え字1つ分左にシフトする。添え字の範囲は0から。
                 // b[0]=8; b[1]=7; b[2]=6; b[3]=5; b[4]=4; b[5]=3; のときには c=b.Shift(-1)は
                 // c[0]=7; c[1]=6; c[2]=5; c[3]=4; c[4]=3; c[5]=0; となる。
    b.Norm1();     // ベクトルbの1ノルムを返す。
    b.Norm2();     // ベクトルbの2ノルムを返す。
    b.NormInfty(); // ベクトルbの無限大ノルムを返す。
    b.Size();    // ベクトルbのサイズ（添え字の最大値）を返す。この場合は5が返る。
    b.Show(k);   // ベクトルの中身を有効数字k桁で表示する。添え字の範囲は1から。

    Matrix A(n,n);   // 行列を宣言、A[0][0]～A[n][n]を使用可能。
    A.Transpose();   // 行列Aの転置行列を返す。添え字の範囲は0から。
    A.RowVector(k);  // 行列Aのk行目をベクトルとして返す。添え字の範囲は0から。
    A.ColVector(k);  // 行列Aのk列目をベクトルとして返す。添え字の範囲は0から。
    A.SetRowVector(k, v); // 行列Aのk行目にベクトルvをセットする。添え字の範囲は0から。
    A.SetColVector(k, v); // 行列Aのk列目をベクトルvをセットする。添え字の範囲は0から。
    A.Gauss(b);           // Gaussの消去法を用いて Ax=b を解いた結果を返す。計算に用いる添え字の範囲は1から。
                          // スケーリングおよび部分ピボット選択を使用。
    A.RowSize();    // 行列Aの行数（添え字の最大値）を返す。この場合は5が返る。
    A.ColSize();    // 行列Aの列数（添え字の最大値）を返す。この場合は5が返る。
    A.Show(k);      // 行列Aの中身を有効数字k桁で表示する。添え字の範囲は1から。

固有値関係、LU分解、疎行列なんかも実装したいところですが、今のところそこまでは必要を感じていませんので未実装です。

使用例

3次元ニュートン法

#include <stdio.h>
#include <math.h>
#include "Linear.h"

// x^2+y^2+z^2=1
// y=sin(x)
// z=x+y
// を解く。プログラムではxをx[1]、yをx[2]、zをx[3]とする。
// F(x) = ( x[1]^2 + x[2]^2 + x[3]^2 - 1, x[2] - sin(x[1]), x[3]-x[1]-x[2] )

// 関数値
Vector f(Vector x)
{
    Vector y(3);
    y[1] = x[1]*x[1] + x[2]*x[2] + x[3]*x[3] - 1;
    y[2] = x[2] - sin(x[1]);
    y[3] = x[3] - x[1] - x[2];
    return y;
}
// ヤコビ行列
Matrix df(Vector x)
{
    Matrix y(3, 3);
    y[1][1] = 2*x[1];     y[1][2] = 2*x[2]; y[1][3] = 2*x[3];
    y[2][1] = -cos(x[1]); y[2][2] = 1;      y[2][3] = 0;
    y[3][1] = -1;         y[3][2] = -1;     y[3][3] = 1;
    return y;
}

int main()
{
    int n = 3;
    Matrix J(n, n);
    Vector x(n);

    // 初期値
    x[1]=0.5; x[2]=0.5; x[3] = 0.5;
    for(int k=1; k<=10; k++){
        J = df(x);
        x = x - J.Gauss(f(x));
        printf("x^(%d)=(%.15f,%.15f,%.15f)\n",k,x[1],x[2],x[3]);
    }

    return 0;
}

ロトカ・ボルテラ方程式

#include <stdio.h>
#include "Linear.h"

// x_1'=(a-b*x_2)*x_1  を解く
// x_2'=(-c+d*x_1)*x_2
Vector f(double t, Vector x)
{
    Vector y(2);
    y[1] = ( 1.0 - 1.0*x[2])*x[1];
    y[2] = (-1.0 + 0.5*x[1])*x[2];
    return y;
}

int main()
{
    double a = 0.0, b = 30.0; // 計算範囲
    int N=1000; // 分割数
    int d = 2; // 従属変数の次元
    Matrix S(d, N);
    Vector t(N);
    Vector x(d), k1(d), k2(d), k3(d), k4(d);
    
    double h = (b-a)/N; // 区間幅
    x[1] = 0.5; x[2] = 0.5; // 初期条件
    S.SetColVector(0, x);
    for(int n=0; n<=N; n++){
        t[n] = a + n*h;
    }

    for(int n=0; n<N; n++){
        k1 = h * f(t[n],       x       );
        k2 = h * f(t[n] + h/2, x + k1/2);
        k3 = h * f(t[n] + h/2, x + k2/2);
        k4 = h * f(t[n] + h/2, x + k3  );
        x += (k1 + 2*k2 + 2*k3 + k4)/6;
        S.SetColVector(n+1, x);
    }
    FILE *fp;
    fp = fopen("l-v.dat", "w");
    for(int n=0; n<=N; n++){
        fprintf(fp, "%f", t[n]);
        for(int i=1; i<=d; i++){
            fprintf(fp, " %f", S[i][n]);
        }
        fprintf(fp, "\n");
    }
    fclose(fp);

    return 0;
}