2008年6月27日 星期五

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

Standard
業界上常用的 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));
}

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

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