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

算 determine 不只一種方法，也有人用餘因子方式求，

但筆者認為效率較不佳，還是以高斯進行上三角化求。

演算法適用於 complex matrix。

範例碼

/************************************************************************/
/* */
/* filename : gauss_elimination_determine.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);

double mat_det(double** mat, size_t dim, const double eps)
{
    if(mat==NULL || !dim){
        return 0.0;
    } else {
        size_t i, j, k, is, js;
        double *row_ptr, **mat_tmp;
        double tmp, pivot, det, flag;

        // avoid to destroy source matrix
        mat_tmp = mat_new(dim, dim);
        if(!mat_tmp) {
            puts("allocate fail.");
            return 0.0;
        }
        memcpy((void*)*mat_tmp, (void*)*mat, sizeof(double)*dim*dim);

        flag = det = 1.0;
        for(k=0; k<dim-1; ++k){
            // find pivot
            for(pivot=0.0, i=k; i<dim; ++i) for(j=k; j<dim; ++j){
                tmp = fabs(mat_tmp[i][j]);
                if(tmp > pivot)
                    pivot = tmp, is=i, js=j;
            }
            // check pivot
            if(pivot<=eps){
                free((void*)mat_tmp);
                return 0.0;
            }
            // row swap
            if(is!=k){
                flag = -flag;
                row_ptr = mat_tmp[is];
                mat_tmp[is] = mat_tmp[k];
                mat_tmp[k] = row_ptr;
            }
            // col swap
            if(js!=k){
                flag = -flag;
                for(i=k; i<dim; ++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<dim; ++i){
                tmp = mat_tmp[i][k] / mat_tmp[k][k];
                for(j=k+1; j<dim; ++j)
                    mat_tmp[i][j] -= tmp * mat_tmp[k][j];
            }
            det*=mat_tmp[k][k];
        }
        det*=flag*mat_tmp[dim-1][dim-1];
        free((void*)mat_tmp);
        return det;
    }
}

int main()
{
    const double eps = 1E-9;
    const size_t dim = 4;
    double **m1 = mat_new_arg(dim,dim, 
         1., 2., 2., 1., 
         2., 7., 5., 2.,  
         1., 2., 4., 2.,
        -1., 4.,-6., 3.);
    puts("m1 ="), mat_print(m1, dim, dim);
    printf("det(m1) = %lf\n", mat_det(m1, dim, eps));
    free((void*)m1);
    getchar();
    return 0;
}