如何用GSL數學庫求解矩陣方程組？

無雨森發表於2019-05-11
矩陣
gsl_matrix *CalculateAxb(const gsl_matrix *A,const gsl_matrix *b)
{
    int i;
    int size1=A->size1;
    gsl_matrix *tempA=gsl_matrix_alloc(size1,size1);
    gsl_matrix_memcpy(tempA,A);
    gsl_matrix *tempb=gsl_matrix_alloc(1,size1);
    gsl_matrix_memcpy(tempb,b);

    gsl_vector *vectorbx=gsl_vector_alloc(size1);
    for(i=0;i<size1;i++)
        gsl_vector_set(vectorbx,i,gsl_matrix_get(tempb,0,i));
    gsl_vector *x=gsl_vector_alloc(size1);
    int s;
    gsl_permutation *px=gsl_permutation_alloc(size1);
    gsl_linalg_LU_decomp(tempA,px,&s);
    gsl_linalg_LU_solve(tempA,px,vectorbx,x);
    gsl_matrix *temp=gsl_matrix_alloc(1,size1);
    for(i=0;i<size1;i++)
        gsl_matrix_set(temp,0,i,gsl_vector_get(x,i));
    gsl_matrix_free(tempA);
    gsl_matrix_free(tempb);
    gsl_vector_free(vectorbx);
    gsl_vector_free(x);
    gsl_permutation_free(px);
    return temp;
}

講解

gsl_linalg_LU_decomp

函式原型：

int gsl_linalg_LU_decomp(gsl_matrix *A,gsl_permutation *p,int *signum)

1.這個函式將矩陣A進行LU分解。

線上性代數中，LU分解是矩陣分解的一種，可以將一個矩陣分解為一個下三角矩陣和一個上三角矩陣的乘積（有時是它們和一個置換矩陣的乘積）。

這個函式的LU分解則是將一個矩陣和置換矩陣的乘積分解為一個下三角矩陣和一個上三角矩陣即PA=LU。

P:置換矩陣;
L:the lower triangular matrices;
U:the upper triangular matrices.

2.LU分解的演算法

部分選主元的高斯消元法 Gaussian Elimination with partial pivoting
(Golub & Van Loan, Matrix Computations, Algorithm 3.4.1)

相關文章