用 C 語言實作矩陣運算的 Function

業界上常用的 Matlab 軟體，因為本身就是為了做數學運算而設計，所以在數值分析和處理方面非常快速好用，但是在 C 語言就中沒有那些好用的指令和數學函式庫了，得自己打造所有的運算流程和架構。舉例來說，使用 C 語言來做矩陣運算，程式開發者都是先『指定行列數量』建立一個陣列，再『指定行列數量』去對內部的值，一一做處理或高斯消去解方程式，這比起 Matlab 中用兩行就可以搞定一切運算，可說易用性相差甚遠。

所以就想到，要是我們能設計出一套好用的 C 語言矩陣函式庫，是否就可以省下這些重覆設計矩陣運算流程的時間？於是可以做一些實作，來完成這個任務。

基本的矩陣資料結構：
typedef struct {
   int rows;
   int cols;
   int total;
   double *table;
} Matrix;


其中定義了幾個基本的 functions 做矩陣的處理：
建立矩陣：
Matrix *matrix_create(int rows, int cols)
矩陣最佳化：
Matrix *matrix_optimize(Matrix *matrix)
解矩陣的方程式：
Matrix *matrix_solution(Matrix *matrix)
顯示並列出矩陣內容：
void matrix_print(Matrix *matrix)
清空矩陣：
void matrix_clear(Matrix *matrix)


函式的程式碼：
void matrix_clear(Matrix *matrix)
{
   int i, j;
   double *ptr;

   /* initializing matrix */
   ptr = matrix->table;
   for (i=0;i<matrix->total;i++)
       *(ptr++) = 0;
}

Matrix *matrix_create(int rows, int cols)
{
   Matrix *matrix;

   /* allocate */
   matrix = (Matrix *)malloc(sizeof(Matrix));
   matrix->rows = rows;
   matrix->cols = cols;
   matrix->total = rows * cols;
   matrix->table = (double *)malloc(sizeof(double)*rows*cols);
   matrix_clear(matrix);

   return matrix;
}

Matrix *matrix_optimize(Matrix *matrix)
{
   Matrix *new_matrix;
   int i, j;
   double *ptr, *o_ptr;

   /* create a new matrix */
   new_matrix = matrix_create(matrix->rows, matrix->cols);

   /* copy old matrix to the new one */
   ptr = new_matrix->table;
   o_ptr = matrix->table;

   for (i=0;i<matrix->total;i++)
       *(ptr++) = *(o_ptr++);

   return new_matrix;
}

Matrix *matrix_solution(Matrix *matrix)
{
   Matrix *result;
   Matrix *source;
   int i, j, k;
   double rate;
   double *ptr;

   /* 對矩陣做最佳化 */
   source = matrix_optimize(matrix);

   ptr = source->table;
   for (i=0;i<source->rows-1;i++) {
       if (*(ptr+i*source->cols+i)==0) {
           /* ignore the column */
           continue;
       }

       for (j=i+1;j<source->rows;j++) {
           rate = *(ptr+j*source->cols+i) / *(ptr+i*source->cols+i);
           *(ptr+j*(source->cols)+i) = 0;

           for (k=i+1;k<source->cols;k++) {
               *(ptr+j*source->cols+k) -= *(ptr+i*source->cols+k) * rate;
           }
       }
   }

   for (i=source->rows-1;i>0;i--) {
       if (*(ptr+i*source->cols+i)==0)
           continue;

       for (j=i-1;j>=0;j--) {
           rate = *(ptr+j*(source->cols)+i) / *(ptr+i*source->cols+i);
           *(ptr+j*(source->cols)+i) = 0;
           *(ptr+j*(source->cols)+source->cols-1) -= *(ptr+i*source->cols+source->cols-1) * rate;
       }
   }

   /* 建立空白的矩陣 */
   result = matrix_create(source->rows, 1);

   ptr = source->table;
   for (i=0;i<source->rows;i++) {
       if (*(source->table+i*source->cols+i)==0)
           continue;

           *(result->table+i) = *(ptr+i*source->cols+source->cols-1) / *(source->table+i*source->cols+i);
   }

   /* release */
   free(source->table);
   free(source);

   return result;
}

另外再設計兩個 Quadratic Splines 的初始化並填值和計算加總函式，以測試矩陣函式：
void matrix_init_qs(Matrix *matrix, Matrix *xftable)
double quadratic_splines(Matrix *matrix, Matrix *xftable, double x)

Quadratic Splines 函式的程式碼：
void matrix_init_qs(Matrix *matrix, Matrix *xftable)
{
    int i;
    int cols = xftable->rows-1;
    int rows = xftable->rows-1;
    double *ptr = matrix->table;
    double *xf = xftable->table;

    for (i=0;i<rows;i++) {
        *(ptr+i) = *(xf+(i+1)*xftable->cols) - *(xf+i*xftable->cols);  /* bi */

        if (i!=0) /* c1 = 0 */
            *(ptr+i+cols) = pow(*(xf+(i+1)*xftable->cols) - *(xf+i*xftable->cols), 2); /* ci */

        *(ptr+matrix->cols-1) = *(xf+(i+1)*xftable->cols+1) - *(xf+i*xftable->cols+1); /* f */
        ptr += matrix->cols;
    }

    /* insert a empty row into the matrix for c1 = 0 */
    ptr += matrix->cols;
    for (i=0;i<cols-1;i++) {
        *(ptr+i) = 1;     /* bi */
        *(ptr+i+1) = -1;    /* b(i+1) */

        if (i!=0) /* c1 = 0 */
            *(ptr+i+cols) = 2 * (*(xf+(i+1)*xftable->cols) - *(xf+i*xftable->cols)); /* ci */

        ptr += matrix->cols;
    }
}

double quadratic_splines(Matrix *matrix, Matrix *xftable, double x)
{
   int i;

   /* 找到 X 所在的範圍 */
   for (i=0;i<xftable->rows;i++) {

       if (x>*(xftable->table+i*xftable->cols))
           continue;

       i--;
       break;
   }

   /* 代入 QS 的公式 */
   return *(xftable->table+i*xftable->cols+1) \
           + *(matrix->table+i) * (x - *(xftable->table+i*xftable->cols)) \
           + *(matrix->table+xftable->rows-1+i) * pow(x - *(xftable->table+i*xftable->cols), 2);

}

使用實例：

#define FNUM 4

void main(void)
{
   Matrix *xftable;
   Matrix *matrix;
   Matrix *solution;
   int i, j;
   double fx[FNUM][2] = {{3.0, 2.5}, {4.5, 1.0}, {7.0, 2.5}, {9.0, 0.5}};
   double *ptr;

   /* 建立一個空矩陣 */
   xftable = matrix_create(FNUM, 2);

   /* 將陣列內的數值填入空矩陣 */
   ptr = xftable->table;
   for (i=0;i<xftable->rows;i++) {
       for (j=0;j<xftable->cols;j++) {
           *(ptr++) = fx[i][j];
       }
   }

   /* 列出目前矩陣內容 */
   matrix_print(xftable);

   /* 針對 Quadratic Splines 建立一個空矩陣 */
   matrix = matrix_create((FNUM-1)*2, (FNUM-1)*2+1);

   /* 針對 Quadratic Splines 初始化並填入數值 */
   matrix_init_qs(matrix, xftable);

   /* 列出目前 QS 矩陣內容 */
   matrix_print(matrix);

   /* 解方程式 */
   solution = matrix_solution(matrix);

   /* 列出目前所求得的解 */
   matrix_print(solution);

   /* 代入並顯示利用 QS 所估算 X=5.0 時的數值 */
   printf("s(%f) = %f\n", 5.0, quadratic_splines(solution, xftable, 5.0));
}

嚴格說起來『矩陣最佳化』函式，並不能稱之為最佳化，因為它並不是將空白『行』或『列』的部份刪除以減輕運算，而是補上空白『行』或『列』使矩陣可順利做消去法。雖然刪除比加入好，但因做法上要平移矩陣內數值比較複雜了點，所以就不傷這腦筋只求可運作，不過，有興趣的人可以自行改進。

話說回來，這矩陣函式庫只是用來解解聯立方程式的話，還真是一點用處也沒有呀，應該還要再加上一些實用的運算函式才是。

留言

張貼留言

這個網誌中的熱門文章

有趣的邏輯問題：是誰在說謊

開心的是期中考接近尾聲，而且今天也是我的生日，更開心的是老天爺待我很好，放我一天假沒給我排考。今天，在網路上逛啊逛，看到了一些判斷說謊的邏輯問題，雖然之前就曾看過，但不管再看多少遍都還是覺得很有趣。也想到好幾天沒寫東西了，就貼上來讓大伙一同分享吧！準備好接受考驗了嗎？ Ready? Go! 第一題： A 說B 在說謊，B 說 C 在說謊，C 說 A、B 都在說謊。請問到底誰在說謊？ 第二題： 現在審問四名竊賊嫌疑犯。已知當中有一名是竊賊，還知道這四個人不是誠實就是說謊，請根據他們回答問題的結果中，判斷誰是竊賊。 甲說：『乙沒有偷，是丁偷的。』 乙說：『我沒有偷，是丙偷的。』 丙說：『甲沒有偷，是乙偷的。』 丁說：『我沒有偷。』 解答：一：只有 B 說真話。二：乙是竊賊。
繼續閱讀

Web 技術中的 Session 是什麼？

圖片
很久沒更新了，自從十二月底退伍之後，一直忙著將書完稿給出版社，沒空坐下來好好在 Blog 上寫點東西。沒辦法，在當兵時閒不下來，就答應了出版社的邀約開始寫 Node.js 的書，因此欠下了書的債。還好，如預期的寫完了，預計這一兩個月應該就會上市，還請各界多多捧場！ 會寫這篇文章的原故，是因為在 Facebook 的社群上，有人問起了關於 Session 和 Cookie 相關性的問題。其實，社群的人都很好，有許多熱心的網友願意幫助別人。不過，為了幫忙回答問題，大家七嘴八舌的東說一句西說一句，甚至引起莫明的論戰，實在不是一個好的現象。況且，使用不同語言、有不同技術背景的人，都各自用自己的角度去描述和解說，讓原本就已經不懂的人更有如瞎子摸象，見不到全貌，搞更糊塗了。對發問的人來說，也不知該聽誰的好。 說來湊巧，筆者剛好在剛截稿的 Node.js 書籍中，對 Session 概念也有許多著墨。所以之前稍微在網路上查過，發現鮮少有人提出對 Session 技術的完整說明，而多半是討論 ASP、.Net、Java、PHP 等語言或框架中如何使用 Session，以及一些延伸的特性，以致很多人對 Session 的概念相當片面。因此這次的事件後，一些網友也建議我把在 Facebook 上回應的完整內容，在 Blog 發表記錄下來，讓更多人對 Session 有個比較完整且正確的認識。 本文目標 在 Web 的世界裡，Session 已經是被廣泛使用的一種技術，大多數的 Web 開發者，肯定都能運用自如。在現今的各類 Web 應用程式的開發框架和工具中，Session 也已經被包裝成容易使用的模式，所以本文也就不再多討論 Session 怎麼去使用，而是著重在 Session 的原理和實現上。要知道，無論是使用哪一種語言，這 Session 概念上都是通用的，沒有什麼不一樣。 Session 是什麼？ Session 之所以會存在，是因為 HTTP 為 stateless 的設計，Server 和 Client 不會一直保持連線狀態，也不會有雙方狀態的即時更新，所以，Server 並不知道 Client 的狀態（像是是否已經登入）。因此，後來的網站開發者，採用 Session 這樣的設計來解決這問題。有趣的是，Se
繼續閱讀

淺談 USB 通訊架構之定義（一）

圖片
USB(Universal Serial Bus) 的應用範疇以及價值已經不需要再多做說明，其身影從 HID(Human Interface Device)、隨身碟（USB Stick、Flash drive、Pen drive）到各種類比式或非連續訊號接收器（如： Webcam、Microphone、DVB Receiver），更進一步可以講到各種通訊裝置或外接設備，甚至是取代和模擬舊式通訊傳輸線（Serial Port）。USB 所有的 Spec 和相關資訊 ，可參考 http://www.usb.org/ 。 USB 整體運作模式 一般來說 USB 的通訊結構有如 Server/Client，以 PC 上的情形為例，位於主機上的 USB 裝置稱為『USB Host』，我們可以在上面外接上數個裝置（與 USB Host 相連的裝置通常被稱為 USB Device 或 USB Client）。 底層上，『Host 負責主導整個 USB 結構的通訊』，它會輪詢所有的 USB Deivce 以檢查是否有裝置需要傳送資料，所有 USB Device 都必需要等待 Host 的命令，唯有 Host 同意時，USB Device 才可以開始傳送資料。 Host 的行為通常是由 USB Host Controller 在做控制，其控制器是一組控制晶片，而對於現代 PC ，USB Host Controller 的晶片通常都被 Layout 在主機板上，以 PCI Bus 的形式存在供作業系統控制，我們可以使用 lspci 看到它的存在： 00:1a.0 USB Controller: Intel Corporation 82801I (ICH9 Family) USB UHCI Controller #4 (rev 03) 00:1a.1 USB Controller: Intel Corporation 82801I (ICH9 Family) USB UHCI Controller #5 (rev 03) 00:1a.2 USB Controller: Intel Corporation 82801I (ICH9 Family) USB UHCI Controller #6 (rev 03) 00:1a.7 USB Controller: Intel Corporati
繼續閱讀

淺談 USB 通訊架構之定義（二）

圖片
之前『 淺談 USB 通訊架構之定義（一） 』已經提過 USB 大致上的概念和一些特性，以及如何從 Linux 上去嘗試驗證 USB Device 上的 Information，但還尚未談及 USB Device 是如何與 USB Host 做溝通。而且，初始化的過程也很重要，該過程會讓一個 USB Device 描述自己的『模式(Configuration)』和『介面(Interface)』以及更多配置上的關鍵訊息。雖然這些初始化過程和定義繁雜，但是拜 OS 優良的包裝所賜，一般驅動程式都可以假設已知或是透過簡單的 API 去取得這些資訊以進行開發，因為大多數初始化動作會自動由 OS Kernel 和 Firmware 所完成。 簡單的 USB Device 初始化流程 插上 USB 裝置 Host 請 裝置 報告 『 裝置 描述資訊 (Device Descriptors )』 （此過程又稱為 Enumeration ） 由 Device Descriptors 得知 『 模式的數量 (bNumConfigurations) 』 Host 請 裝置 回報 『 模式的描述資訊( Configuration Descriptors) 』 由 Configuration Descriptors 得知 『 介面的數量( bNumInterfaces ) 』 Host 請 裝置 回報 『 介面 的描述資訊( Interface Descriptors) 』 由 Interface Descriptors 得知 『 端點的數量( bNumEndpoints ) 』 Host 請 裝置 回報 『 端點 的描述資訊( Endpoint Descriptors) 』 由 Endpoint Descriptors 得知 『 該端點的 資料傳輸模式( bmAttributes ) 』等訊息 等待驅動程式進行後續處理 整個 USB Device 的描述結構可以圖表示： 當然，若不是 USB Device 設計者或是 USB Host Controller 晶片的 Firmware/Driver 開發人員，可以不用知道這麼多過程細節，因為無論是 Linux 還是其他作業系統，它們都會將這些最低階的過程包裝起來，等著驅動程式去使用，而 Linux Kernel 上負責此任務的就被稱為
繼續閱讀

Reverse SSH Tunnel 反向打洞實錄

近來經手了幾個案子，其目標都是設計一台做特定用途的系統，但往往這台機器都會被鎖在防火牆後面，甚至不會連上網際網路，是完全封閉的環境。每當接近結案，將系統交到客戶手中以做最後的壓力測試時，便是痛苦的開始。怎麼說呢？因為當面對最後的大量測試，各種大大小小的臭蟲(Bugs)就會原形畢露，一一浮上臺面。而每次問題發生，便要千里迢迢趕到客戶那了解情況、解決問題。如果問題不大倒還好，一旦問題嚴重，被客戶關到三更半夜也是常有的事，更誇張的還必須每天都去報到。所以，為避免車程來回的時間和人力浪費，是否有可以遠端連線做處理呢？便思考起這問題。 雖然這些機器通常都不會有真實 IP，客戶也不會為了我們特別請 MIS 去開啟 NAT 或防火牆，但還好 SSH 提供了反向(Reverse)的機制讓我們可以連進去。其做法就是讓在客戶那的系統，透過 SSH 先連回我們自家有真實 IP 的主機，然後建立反向的通道即可。如此，我們便可從自家的主機，從 SSH 連回放在客戶那的系統。 實際的操作步驟： # 首先，在客戶那理的機器下指令連回我們自己的 Server，並設定自己 Server 上的 12345 port 會對應到幾器上的 SSH port ssh -NfR 12345:localhost:22 fred@myhost.com # 然後在 myhost 的機器上連自己的 12345 port，就可以連回在客戶那的機器 ssh localhost -p 12345 後記 可以將建立 Reverse SSH tunnel 的命令設定在開機後執行，並且透過設定 ssh config 讓斷線時自動重連，如此就可以確保 SSH 連線不會中斷。此外，若是客戶端那的機器是處於完全無法對外連線的環境，則可以考慮使用 3G 網卡，並將網卡暫時借放於客戶那，然後等結案後再取回，雖然 3G 通訊月租不便宜，但和經常來回的通勤費和時間花費相比，可滑算多了。
繼續閱讀