csr矩阵A利用选主元lu分解求解ax=b的c代码
时间: 2023-08-04 22:09:19 浏览: 60
以下是使用选主元 LU 分解求解稀疏矩阵 Ax=b 的 C 代码:
```
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 100 // 矩阵的最大维数
#define eps 1e-10 // 用于判断是否为零的阈值
typedef struct {
int i, j;
double val;
} node;
int n; // 矩阵的维数
node a[N * N]; // 稀疏矩阵存储结构
int ia[N + 1], ja[N * N + 1]; // 行和列指针
// 选主元 LU 分解函数
void lu_decomp(double *A, int *ipiv) {
int i, j, k, p;
double maxa, tmp;
for (i = 0; i < n; i++) {
ipiv[i] = i;
}
for (k = 0; k < n - 1; k++) {
maxa = fabs(A[k * n + k]);
p = k;
for (i = k + 1; i < n; i++) {
if (fabs(A[i * n + k]) > maxa) {
maxa = fabs(A[i * n + k]);
p = i;
}
}
if (p != k) {
for (j = 0; j < n; j++) {
tmp = A[k * n + j];
A[k * n + j] = A[p * n + j];
A[p * n + j] = tmp;
}
i = ipiv[k];
ipiv[k] = ipiv[p];
ipiv[p] = i;
}
for (i = k + 1; i < n; i++) {
A[i * n + k] /= A[k * n + k];
for (j = k + 1; j < n; j++) {
A[i * n + j] -= A[i * n + k] * A[k * n + j];
}
}
}
}
// LU 分解求解函数
void lu_solve(double *A, int *ipiv, double *b, double *x) {
int i, j;
double sum;
for (i = 0; i < n; i++) {
sum = b[ipiv[i]];
for (j = 0; j < i; j++) {
sum -= A[i * n + j] * x[j];
}
x[i] = sum;
}
for (i = n - 1; i >= 0; i--) {
sum = x[i];
for (j = i + 1; j < n; j++) {
sum -= A[i * n + j] * x[j];
}
x[i] = sum / A[i * n + i];
}
}
int main() {
int i, j, k, cnt = 0;
double A[N][N], b[N], x[N];
int ipiv[N];
// 读入矩阵和向量
scanf("%d", &n);
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
scanf("%lf", &A[i][j]);
if (fabs(A[i][j]) > eps) {
a[cnt].i = i;
a[cnt].j = j;
a[cnt].val = A[i][j];
cnt++;
}
}
}
for (i = 0; i < n; i++) {
scanf("%lf", &b[i]);
}
// 构造行和列指针
k = 0;
for (i = 0; i <= n; i++) {
ia[i] = k;
for (j = 0; j < n; j++) {
if (a[k].i == i) {
ja[k] = a[k].j;
k++;
}
}
}
// LU 分解求解
lu_decomp(a, ipiv);
lu_solve(a, ipiv, b, x);
// 输出结果
for (i = 0; i < n; i++) {
printf("%lf ", x[i]);
}
return 0;
}
```
其中,稀疏矩阵存储结构使用了三元组表示法,同时构造了行和列指针。选主元 LU 分解的代码实现和普通的 LU 分解类似,只是在选主元时需要额外的代码。LU 分解求解的代码也和普通的 LU 分解类似,只是需要传入行和列指针。
相关推荐
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)