[Matrix] 矩陣與向量乘積

最近對線性代數應用感興趣，開始研究線性代數運算的底層程式實作。
先從矩陣乘積開始。
在此用 C++ 做矩陣與向量的乘積。

我們定義以下的矩陣 Class：

    template <class T, size_t, size_t j>
    class mtx
    {
    public:
       mtx(){}
       ~mtx(){}

       inline T& operator()(size_t i, size_t j)
       {
            return arr_[i][j];
       }
       void info()
       {
            cout << endl;
            for(size_t i = 0; i < ROWS; i++)
            {
                for(size_t j = 0; j < COLS; j++)
                    cout << " " << arr_[i][j];
                cout << endl;
            }
            cout << endl;
       }

    private: 
       T arr_[ROWS][COLS];
    };

接著我們就可以使用此 Class 來做矩陣與向量乘積：

    mtx a;
    mtx x;
    mtx b;

    a(0, 0) = 2;
    a(0, 1) = 1;

    a(1, 0) = 3;
    a(1, 1) = 5;

    a(2, 0) = 4;
    a(2, 1) = 7;

    x(0, 0) = 10;
    x(1, 0) = 12;
    
    // matrix mul:
    for(size_t r = 0; r < 3; r++)    b(r, 0) = 0;

    for(size_t i = 0; i < 3; i++)
    {
        for(size_t j = 0; j < 2; j++)
            b(i, 0) = b(i, 0) + a(i, j) * x(j, 0);
    }

    cout << " a:" << endl;
    a.info();
    cout << " x:" << endl;
    x.info();
    cout << " ax = b, b:" << endl;
    b.info();

執行結果
 a:

 2 1
 3 5
 4 7

 x:

 10
 12

 ax = b, b:

 32
 90
 124

做乘法的部份，重點在這段 Code：

    for(size_t r = 0; r < 3; r++)    b(r, 0) = 0;

    for(size_t i = 0; i < 3; i++)
    {
        for(size_t j = 0; j < 2; j++) // row 水平方向掃描
            b(i, 0) = b(i, 0) + a(i, j) * x(j, 0);
    }

這段 Code 在執行 ax = b 的矩陣乘積。這種寫法視 a 矩陣跟向量 x 做內積。
若把兩個 for 迴圈對調，變成 x 內的每個元素，先對 a 矩陣的 column vector 內每個分量先相乘，我們觀看一下寫法：

    for(size_t r = 0; r < 3; r++)    b(r, 0) = 0;

    for(size_t j = 0; j < 2; j++)
    {
        for(size_t i = 0; i < 3; i++) // column 垂直方向掃描
            b(i, 0) = b(i, 0) + a(i, j) * x(j, 0);
    }

兩種寫法的執行結果相同。但執行速度根據矩陣儲存方式會有差異，因為二維陣列在記憶體內儲存方式，是同 row 的值存在同一區域。
因此我推論：
上述第二種寫法在 a 矩陣有很多個 row 時，執行速度較慢，因為內層迴圈讓 a 矩陣取值一次跳一個 row。
而第一種寫法，內層迴圈只讓 a 取值一次跳一個 column，在相鄰記憶體位置取值。而對於 x 和 b 向量而言，都是在相鄰記憶體位置讀寫。

我之後會補充矩陣向量乘積執行速度的實驗在此。