線形演算ライブラリ

    数値解析の授業でC言語を用いてプログラミングの演習を行っていますが、行列やベクトルを用いる必要があります。線形演算のライブラリには、Lapack、uBlas、Eigen等々がありますが、授業の中でインストールには時間をかけたくないですし、特に最近のオンライン授業では生徒への個別対応が難しいということから、様々な環境で確実に動くものが欲しいということで、簡単なものを自作してみました。
    ただし、あくまで簡易的なものですので、機能は限定的ですし、計算効率も良くありません。本格的な数値計算には用いないようにして下さい。とはいえ、ちょっとした計算をお試しで実行してみるには便利ですので、紹介することにしました。
    使い方は、Linear.h を cpp ファイルと同じ場所に置いて、プログラムの最初に #include "Linear.h" でインクルードするだけです。
    以下、特徴です。
    基本的な機能しか用いていないので、大抵のバージョンのC++で動くはずです。ただ、書き方が古すぎるということで警告は出るかもしれません。
    すべてを1つのヘッダファイルに記述しているので、使いたいプログラムからインクルードするだけで済みます。C言語のファイル管理の慣習からは外れますが、使い易さを優先するため、この形にしました。

前バージョンとの違い

    2025/4/14にバージョンアップを行いました。以前のものとの主な違いは以下の通りです。
    行列やベクトルの代入やコピーコンストラクタ実行の際に、新規のオブジェクトを作成するのを止め、参照を渡すことにしました。これにより、関数呼び出しの際のオーバーヘッドを低減させることができます。
    その影響により、ある変数に格納されている行列やベクトルを他の変数に代入したとしても、実体としては同じものになります。すなわち、一方を変更すると他方も変化してしまいます。Python の list 等と同じ動きになります。
    実体が別のコピーを必要とする場合もあるかと思いますので、その目的のためにメソッド Copy() を実装しました。
    行列やベクトルを新規に作成したときに、ゼロクリアされるのを止めました(ゼロクリアが必要無い場合には無駄な動きになるため)。
    その代わり、ゼロクリアを行うメソッド Clear() を実装しました。

具体的な機能

    型(クラス)としては、ベクトルを表す Vector、整数値を取るベクトルを表す IntVector、行列を表す Matrixが使用できます。
    
        Vector b(100); // ベクトルを宣言
        IntVector p(30); // 整数型のベクトルを宣言
        Matrix A(10,20); // (行数,列数)で行列を宣言
    
    のように、行列やベクトルのサイズを指定して宣言することもできますし、
    
        Matrix A;
        A = B;
        Vector b = c;
    
    のように、とりあえず空の状態で宣言して後から代入したり、宣言と同時に別の行列やベクトルで初期化したりしても構いません。なお、A = B とすると行列 A と行列 B は同じ実体を持つことになるので、一方を変更すると他方も変化してしまいます。B とは別の実体を持つコピーを得たい場合は、A = B.Copy() とします。
    サイズを指定して宣言した場合、中身は初期化されないことに注意が必要です。添え字の範囲は、0から指定した値までで、普通の配列のように使用することができます。例えば、Matrix A(20,30); と宣言した場合には A[0][0]~A[20][30] が使えます。添え字の範囲についてはチェックが入り、範囲を外れた添え字でアクセスするとエラーが出ます。
    Vector および Matrix については、演算子オーバーロードにより、"+", "-", "+=", "-=", "*", "*=", "/", "/=" がオーバーロードされているため、自然な演算は自由に行えます。ただし、VectorとVectorの積は内積となります。また、Matrix同士の積、VectorとMatrixの積、Vector同士の積(内積)については、計算される添え字の範囲は1からになります。それ以外の演算については添え字は0から計算されます。
    IntVector については、添え字を指定してのアクセス、IntVectorからIntVectorへの代入、以外の演算はできません。
    その他の機能については以下を見て下さい。
    
        int n = 5;
    
        Vector b(n); // ベクトルを宣言、b[0]~b[n]を使用可能。
        b.Shift(2);  // ベクトルの中身を添え字2つ分右にシフトしたベクトルを返す。
                     // 添え字の範囲は0から。b自体は変化しないことに注意。
                     // 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自体は変化しないことに注意。
                     // 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.Clear(); // ベクトルbの要素を全てゼロでクリアする。
        b.Copy();  // ベクトルbと同じ値を持つが、実体としては異なるベクトルを返す。
        b.Show(k); // ベクトルの中身を有効数字k桁で表示する。kを省略するとk=16とみなされる。表示する添え字の範囲は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.Clear(); // 行列Aの要素を全てゼロでクリアする。
        A.Copy();  // 行列Aと同じ値を持つが、実体としては異なる行列を返す。
        A.Show(k); // 行列Aの中身を有効数字k桁で表示する。kを省略するとk=16とみなされる。添え字の範囲は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()
    {
        Vector x(3);
    
        x[1] = 0.5; x[2] = 0.5; x[3] = 0.5; // 初期値
        for(int k=1; k<=10; k++){
            Matrix 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);
        
        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++){
            Vector k1 = h * f(t[n],       x       );
            Vector k2 = h * f(t[n] + h/2, x + k1/2);
            Vector k3 = h * f(t[n] + h/2, x + k2/2);
            Vector 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;
    }