C语言使用lapack求解线性方程组

安装lapacke库

本文假设使用Ubuntu操作系统。为了简便，使用lapacke库代替lapack。
安装lapacke库：

sudo apt install liblapacke-dev -y

线性方程组

本文以求解以下线性方程组为例，记录如何使用lapack库求解形如 $\mathbf{A}\mathbf{x}=\mathbf{b}$ 的线性方程组。示例方程如下：

$$
\begin{cases}
2x+3z&=1 \\
x-y+2z&=1\\
x-3y+4z&=2\\
\end{cases}
$$

该方程的解为：

$$
\begin{cases}
x&=0.5\\
y&=-0.5\\
z&=0\\
\end{cases}
$$

包含头文件

使用lapacke库的头文件：

#include<lapacke.h>

初始化参数与变量

int n = 3, nrhs = 1, lda = 3, ldb = 1, info;  
int ipiv[3];
double M[9] = {
    2.0,  0.0, 3.0, 
    1.0, -1.0, 2.0, 
    1.0, -3.0, 4.0
};

double V[3] = {
    1.0, 
    1.0, 
    2.0
};

注意其中参数n,nrhs,lda,ldb的意义。

调用函数求解

lapacke使用函数LAPACKE_dgesv求解线性方程组，若成功则返回0：

info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, nrhs, M, lda, ipiv, V, ldb);

其中参数LAPACK_ROW_MAJOR表示按行进行排列，也可以使用LAPACK_COL_MAJOR按列排列，则需要对系数矩阵进行转置。

编译

编写完成后使用gcc编译，如下：

gcc test.c -llapacke

编译完成得到a.out可执行文件，使用命令./a.out执行得到如下结果：

x[0] = 0.500000
x[1] = -0.500000
x[2] = 0.000000

符合正确结果。

完整代码

#include <stdio.h>
#include <stdlib.h>
#include <lapacke.h>
#include <time.h>

#define N 3
#define NRHS 1
#define LDA N
#define LDB NRHS

int main (){
    int n = N, nrhs = NRHS, lda = LDA, ldb = LDB, info;
    int ipiv[N];

    double M[LDA * N] = {
        2.0,  0.0, 3.0, 
        1.0, -1.0, 2.0, 
        1.0, -3.0, 4.0
    };

    double V[LDB * N] = {
        1.0, 
        1.0, 
        2.0
    };

    // Record the start time
    clock_t start = clock();

    // Call LAPACKE_dgesv
    info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, nrhs, M, lda, ipiv, V, ldb);
    if(info){
        exit(-1);
    }
    // Record the end time
    clock_t end = clock();

    // Calculate the elapsed time in seconds
    double elapsed_time = (end - start) / (double)CLOCKS_PER_SEC;

    printf("Solution to the system of linear equations:");
    for (int i = 0; i < N; i++) {
        printf("\nx[%d] = %f", i, V[i]);
    }
    printf("\nElapsed time: %e s\n", elapsed_time);
    printf("\n");
    return 0;
}