[C語言數值分析] 架構 matrix 基本函式－Edison.X. Blog

一些 matrix 之範例碼在 main function 那裡覺得老是放不是重點的 code ，很麻煩，

建議自己先包過一些常用的 function 出來。

Basic unity

double** mat_new(const size_t row, const size_t col)
{
    size_t i;
    double *ptr;
    double** mat = (double**)malloc(sizeof(double) * row * col + 
                                    sizeof(double*) * row);
    if(mat==NULL) return NULL;

    ptr = (double*)( (char*)mat + sizeof(double*)*row );
    for(i=0; i<row; ++i)
        mat[i] = ptr, ptr+=col;
    return mat;
}

double** mat_new_arg(const size_t row, const size_t col, ...)
{
    size_t i, j;
    va_list p;
    double** mat = mat_new(row, col);
    
    if(!mat) return NULL;
    va_start(p, col);
    for(i=0; i<row; ++i) for(j=0; j<col; ++j)
        mat[i][j] = va_arg(p, double);
    va_end(p);
    return mat;
}

void mat_print(const double** mat, const size_t row, const size_t col)
{
    size_t i, j;
    for(i=0; i<row; ++i){
        for(j=0; j<col; ++j){
            printf("% 6.2lf ", mat[i][j]);
        }
        puts("");
    }
    puts("");
}

以前習慣會再多下面這幾個

Rix : 將 mat 第 i 列乘上 x。
Rij : 交換 mat 之第 i 列與第 j 列。
Rijx : 將 mat 第 i 列乘上 x，加到第 j 列。

一樣會做 Cix, Cij, Cijx 版本。

後來沒再要求一定要包這六個 function 原因是，在真正的分析時，這些函式可能會對於不同 matrix 呼叫二次，但這些行為是可以在同一個 for-loop 裡面完成的，所以沒再強求。

至於 matrix 要用一維還是二維表示可以探討一段時間。考慮到 heap 有限時，一維可以表示的範圍稍稍大一些 (只有大一點而已)；二維有個好處，在做 Rij ( swap row ) 時，只需要交換 pointer 即可，不需逐一 assigned (當然也是必須要用 allocate 方式為前提才可行)。

至於 Cij 有沒有辦法交換 pointer 即可達成需求，答案是不可能。同一份程式語言只會有 management in row 或 management in column，不可能二個同時存在。而 C / C++ 是 management in row ，故這想法應是不可行。

matrix 問題大多實數係與複數係之演算法是一樣的，只有一點細節需要注意而已，但若已包成 class 型態，建議不妨以 template 撰之，有些演算法可以一次解掉 for real / for complex 問題。

另以前在寫 matrix 的時候，發現 M$ 似乎有封裝了一份簡易的 Matrix class ( 開發時從 tips 上看到 )，有興趣可到 msdn 上看這篇。

Error Handle

例外處理它永遠是個討論不完的議題，特別是 team work 時一定要確定這部份要從哪些做起，銜接的函式在前端是必做的，那後端要不要再做 double check ? 藝術。

一些較基本的錯誤例外處理，較建議先從 caller 那裡自己做判定，如：維度是否大於 0？rowi 是否小於 row ？pointer 是否為 NULL 等問題。

筆者手邊以前開發時之例外處理範例供參考，可再自定義自己的例外處理部份。

xmat.h 片斷

//////////////////////////////////////////////////////////////////////////
// error handle
// last_error = 0 : success.
// last_error = 1 : dim error, row count / col count error.
// last_error = 2 : allocate memory fail.
// last_error = 3 : input param pointer is NULL
// last_error = 4 : has no solution
// last_error = 5 : have infinite solution
#define MAT_SUCCESS 0
#define MAT_DIM_ERR 1
#define MAT_ALLOCATE_FAIL 2
#define MAT_NULL_POINTER 3
#define MAT_HAS_NO_SOl 4
#define MAT_HAVE_INF_SOL 5
#define MAT_SPEC_MATRIX_ERR 6

// static int last_error ; // in math.c
int   mat_last_error();
const char * mat_last_error_msg();

xmat.c 片斷

//////////////////////////////////////////////////////////////////////////
// error message
static int last_error ; 
int   mat_last_error()
{
    return last_error;
}
const char * mat_last_error_msg()
{
    const char* err_msg[] = {
        "success", 
        "dim error ( row/col = 0 , or ri >= row / ci >= col )", 
        "allocate memory fail",
        "input param is a NULL pointer", 
        "has no solutions", 
        "have infinite solutions"
    };
    return err_msg[last_error];
}

此後所有 xmat.h 裡之 function，都必須對 static int last_error; 變數做嚴格之定義與控管，這部份我倒認為沒什麼大問題。強調的是這種方法在 multi-thread 時會死很慘。

Bad Architecture

第一次開發 xmat.h / xmat.c 時，架構其實弄得不算很好 (當下認為是應是較佳的格局)，

認為可改善地方於此記下 note。

1. xmat.c 太過於包山包海

事實上所有關於 matrix 運算我都寫了進來，每份 matrix 運算還分 for real / for complex，連 basic unity operator 也都加上進來，一份 xmat.h 含蓋了超過 50 個 function declare，而這 50 個 function define，全都在 xmat.c 實作。

如果對於 matrix 議題有要細做的話，建議所有 algorithm 都獨立包成一個 .c 便可。

2. 太過強調例外處理

這問題見個見智。我拿 complex 做除法問題為例，有二份 prototype 做為選擇

int comp_div(struct comp * dst, struct comp a, struct comp b); // a / b 存到 dst, 失敗傳回 0, 成功傳回 1

struct comp comp_div( struct comp a, struct comp b) ; // 直接將 a/b 結果傳回去，不做防呆措施。

一開始在規劃時認為直接 return value 沒做防呆很危險，於是所有函式均採用第一種方法實作下去。

開發到後面，發現沒辦法用結合性真的很麻煩[Lemma]！再選一次架構的話會挑第二種方式實作，防呆措施一律由 caller 保證。

最後建議，matrix 議題真的非常多，光是乘法就可以研究一小陣子，不建議在一個小地方卡著，先把所有的東西搞完一遍後，還有時間再針對有興趣的地方做深入即可。

[Lemma]

又想到一個面試題：為什麼 strcpy 要設計成 return pointer ?

我不知道這題有沒有正確答案，我的想法是：結合性。

len = strlen ( strcpy (buf, "EdisonX") );