用 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));
}

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

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