线性代数
本章描述了求解线性系统的函数。本库提供了直接作用于gsl_vector和gsl_matrix对象的线性代数操作。这些例程使用Golub & Van Loan的《矩阵计算》的标准算法,对有效率要求的可以使用一级和二级BLAS。
本章描述的函数声明在头文件gsl_linalg.h中。
14.1 LU分解
一个一般的m×n矩阵A有一个LU分解
PA = LU
其中P为M×M置换矩阵,L为M×min(M,N), U为min(M,N)×N。对于方阵,L是下单位三角形矩阵,U是上三角矩阵。当M > N时, L是一个单位更低的梯形矩阵大小为M×N。当M < N时, U是大小为M×N的上梯形。对于方阵,这种分解可以将线性方程组Ax = b转化为一对三角形方程组(Ly = Pb, Ux = y),通过前代和后代求解。注意,LU分解对于奇异矩阵是有效的。
int gsl_linalg_LU_decomp(gsl_matrix * A, gsl_permutation * p, int * signum)
int gsl_linalg_complex_LU_decomp(gsl_matrix_complex * A, gsl_permutation * p,
int * signum)
这两个函数将矩阵A进行LU分解,PA = LU。输入矩阵A的上三角(或梯形)和下三角(或梯形)部分包含矩阵U。输入矩阵的下三角(或梯形)部分(不包括对角线)包含L。L的对角线元素是统一的,不存储。
置换矩阵P在输出时被编码在置换矩阵变量p中。矩阵P的第j列由单位矩阵的第k列给出,其中k=pj是置换向量的第j个元素。置换的符号由signum给出。它的值是(-1)n,其中n是置换中的交换数。
分解中使用的算法是带有部分旋转的高斯消去(Golub & Van Loan,矩阵计算,算法3.4.1),结合基于3级BLAS的递归算法(Peise和Bientinesi,2016)。
int gsl_linalg_LU_solve(const gsl_matrix * LU, const gsl_permutation * p,
const gsl_vector * b, gsl_vector * x)
int gsl_linalg_complex_LU_solve(const gsl_matrix_complex * LU,
const gsl_permutation * p, const gsl_vector_complex * b,
gsl_vector_complex * x)
这两个函数使用gsl_linalg_LU_decomp()或gsl_linalg_complex_LU_decomp()给出的将A分解为(LU, p)的LU分解来求出方系统 Ax = b。
int gsl_linalg_LU_svx(const gsl_matrix * LU, const gsl_permutation * p,
gsl_vector * x)
int gsl_linalg_complex_LU_svx(const gsl_matrix_complex * LU,
const gsl_permutation * p, gsl_vector_complex * x)
这两个函数使用预先计算好的A的LU分解(LU, p)来解方程组Ax = b。输入时,x应该包含右边的b,输出时,它被解替代。