[C語言數值分析] rank of matrix－Edison.X. Blog

這份副函式除了檢查 rank 之外，也會拿來被 applicate 成，

檢查是否為 Singular matrix ( 也就是 inverse matrix 是否存在)。

令 matrix[row][col]，對 matrix 進行高斯消去成上三角矩陣，只有幾點要注意。

(1) 最大進行次數 = max(row, col)。

(2) 考慮數值穩定性問題，pivot 不要只是 matrix[k][k]!=0 就好，而是要找整個 matrix 最大值 (絕對值的最大值)。

(3) 稍加修改後可拿來算 determine。

只宣告沒實作的程式碼 [C語言數值分析] 架構 matrix 基本函式。

演算法適用於 complex matrix。

範例碼

/************************************************************************/
/* */
/* filename : gauss_elimination_rank.c */
/* author : EdisonX / Edison.Shih */
/* compiler : Visual C++ 2010 */
/* date : 2012.03.07 */
/* */
/************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>

double** mat_new(const size_t row, const size_t col);
double** mat_new_arg(const size_t row, const size_t col, ...);
void mat_print(const double** mat, const size_t row, const size_t col);

size_t mat_rank(double** mat, size_t row, size_t col, const double eps)
{
    if(mat==NULL || !row || !col){
        return 0;
    } else {
        size_t i, j, k, p;
        size_t is, js, rank;
        double tmp, pivot, **mat_tmp;
        double *row_ptr;
        
        // avoid to destroy source data, copy to mat_tmp
        mat_tmp = mat_new(row, col);
        if(mat_tmp==NULL){
            puts("allocate fail!!");
            return 0;
        }
        memcpy((void*)*mat_tmp, (void*)*mat, sizeof(double)*row*col);
        
        p = (row < col) ? (row) : col;
        for(k=rank=0; k<p; ++k){
            // find pivot
            for(pivot=0.0, i=k; i<row; ++i) for(j=k; j<col; ++j){
                tmp = fabs(mat_tmp[i][j]);
                if(tmp > pivot)
                    is = i, js = j, pivot = tmp;
            }
            // check pivot
            if(pivot < eps){
                free(mat_tmp);
                return rank;
            }
            ++rank;
            if(is!=k){ // row swap
                row_ptr = mat_tmp[is];
                mat_tmp[is] = mat_tmp[k];
                mat_tmp[k]  = row_ptr;
            }
            if(js!=k){ // col swap
                for(i=0; i<row; ++i) {
                    tmp = mat_tmp[i][k];
                    mat_tmp[i][k] = mat_tmp[i][js];
                    mat_tmp[i][js] = tmp;
                }
            }
            // elimination
            for(i=k+1; i<row; ++i){
                tmp = mat_tmp[i][k] / mat_tmp[k][k];
                for(j=k+1; j<col; ++j)
                    mat_tmp[i][j] -= tmp*mat_tmp[k][j];
            }
            
        }
        free(mat_tmp);
        return rank;
    }
}

int main()
{
    const double eps = 1E-9;
    const size_t row = 3, col = 3;
    double **m1 = mat_new_arg(row,col, 
        1.0, 2.0, 3.0, 
        4.0, 5.0, 6.0, 
        7.0, 8.0, 9.0);
    puts("m1 ="), mat_print(m1, row, col);
    printf("rank(m1) = %u\n", mat_rank(m1, row, col, eps));
    free((void*)m1);
    getchar();
    return 0;
}