计算数学
矩阵向量基本运算
C语言基本线性代数运算
矩阵运算是数值运算中的基本工具。这里以C语言为例，开发常见的运算函数，并可以在此基础上进一步扩展。并给出一个线性最小二乘范例。拟合的Makefile文件如下。
CC = gcc
cflag = -o -g
test.exe : test.o matlib.o
$(CC) -o test.exe test.o matlib.o
test.o :test.c
$(CC) -g -c test.c
matlib.o: matlib.c matlib.h
$(CC) -g -c matlib.c matlib.h
clean:
rm -f *.o *.gch test.exe
/*------------------------------------------------------
Versions and Changes :
v1.0----2024.02.24
C 语言线性代数
主要包括
1. 简单的矩阵向量加减乘法运算，坐标旋转矩阵
2. 常用的矩阵分解 LDL, QR(MGS正交化，householder变换)
3. 最小二乘问题求解
4. 一般线性方程求解(正交化方法)
5. 矩阵求逆
6. 其他数值分析方法
--------------------------------------------------------
Notes:
1. 头文件为 matlib.h
2. 调用该函数库需要加入头文件
#include "matlib.h"
3. 编译方式为 gcc main.cpp matlib.cpp
测试环境 HP unix 编译器 HP aCC
移植到其他操作系统及编译器（如g++或Vc++）可能会有少量语法需要修改。
4. 如有bug，请联系 song.yz@foxmail.com
--------------------------------------------------------
Author : 宋叶志
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2024
-------------------------------------------------------*/
#ifndef matlib_h
#define matlib_h
#include
#include
#include
// 定义向量结构体
typedef struct
{
int size; // 向量的大小
double *data; // 指向向量数据的指针
} Vector;
// 定义矩阵结构体
typedef struct
{
int rows; // 矩阵的行数
int cols; // 矩阵的列数
double **data; // 指向矩阵数据的二维指针
} Matrix;
// 初始化向量
Vector createVector(int size);
// 释放向量内存
void freeVector(Vector *vector);
// 初始化矩阵
Matrix createMatrix(int rows, int cols);
// 释放矩阵内存
void freeMatrix(Matrix *matrix);
//-----------------------------------------------------------
//设置向量为一个数值
Vector vecSet(int size, double number);
//复制向量
Vector vecCopy(Vector v1);
//向量所有元素加、减、乘、除一个数值
Vector vecAddNum(Vector v1, double number);
Vector vecMinusNum(Vector v1, double number);
Vector vecMulNum(Vector v1, double number);
Vector vecDivNum(Vector v1, double number);
double vectorNorm(Vector v1);
double vectorDot(Vector v1, Vector V2);
// Vector vecAddVec();
// Vector vecMinusVec();
// 矩阵转置
void matrixTranspose(Matrix A, Matrix *AT);
//矩阵乘以向量
void matMulVec(Matrix matrix, Vector vector, Vector *result);
//矩阵乘以矩阵
void matMulMat(Matrix A, Matrix B, Matrix *C);
//累加法方程
void AddNormalEqation(Matrix *ATA, Vector *ATb, Vector H, double oc);
// 法方程求解最小二乘
void LeastSquareEquation(Matrix ATA, Vector ATb, Vector *X, Matrix *P);
//LDL分解
//根据LDL分解计算对称正定矩阵逆
void covariance(Matrix L, Vector D, Matrix *P);
void LDL(Matrix A, Matrix *L, Vector *D);
//上三角矩阵求解
void uptriEqation(Matrix A, Vector b, Vector *x);
//下三角矩阵求解
void downtriEqation(Matrix A, Vector b, Vector *x);
//下三角矩阵的逆
void inversLowTriMatrix(Matrix L, Matrix *R);
//上三角矩阵的逆
void inversUperTriMatrix(Matrix U, Matrix *R);
//对称正定矩阵逆
void inversSymmetricPositiveMatrix(Matrix A, Matrix *invA);
//旋转矩阵
void rotMatrix(Matrix *Rmat, double angle, char axID);
//打印矩阵
void printMatrix(Matrix matrix, char format);
void printVector(Vector vector, char format);
// lagrange插值
void interpLagrange(Vector x, Vector y, Vector xx, Vector *yy);
#endif
/*------------------------------------------------------
Versions and Changes :
v1.0----2024.02.24
c linear algebra library
--------------------------------------------------------
Notes:
1. 头文件为 matlib.h
2. 调用该函数库需要加入头文件
#include "matlib.h"
3. 编译方式为 $(cc) main.cpp matlib.cpp
--------------------------------------------------------
Author : 宋叶志
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2024
-------------------------------------------------------*/
#include "matlib.h"
#include
#include
Vector createVector(int size)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 20:02
初始化向量，并设置所有元素为0
--------------------------------------------------------
Input Parameters :
size -- 向量长度
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
Vector vector;
vector.size = size;
vector.data = (double*) malloc(size * sizeof(double));
int i;
for (i = 0; i < size; i++)
{
vector.data[i] = 0.0;
}
return vector;
}
void freeVector(Vector *vector)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 20:02
释放向量内存
在C99中允许变长数组，为防止部分编译器不支持，因此采用动态方式分配数组
用法 freeVector(&vec)
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
free(vector->data);
vector->data = NULL;
vector->size = 0;
}
Matrix createMatrix(int rows, int cols)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 20:05
初始化矩阵，并设置所有元素为0.
Changes :
--------------------------------------------------------
Input Parameters :
rows
cols
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
Matrix matrix;
matrix.rows = rows;
matrix.cols = cols;
// 分配指向行指针的内存
matrix.data = (double**) malloc(rows * sizeof(double*));
if (matrix.data == NULL)
{
perror("Error allocating memory for matrix rows");
exit(EXIT_FAILURE);
}
// 分配每行元素的内存
for (int i = 0; i < rows; i++)
{
matrix.data[i] = (double*) malloc(cols * sizeof(double));
if (matrix.data[i] == NULL)
{
perror("Error allocating memory for matrix elements");
// 释放之前已分配的内存
for (int j = 0; j < i; j++)
{
free(matrix.data[j]);
}
free(matrix.data);
exit(EXIT_FAILURE);
}
}
int i, j;
for (i = 0; i < rows; i++)
for (j = 0; j < cols; j++)
{
matrix.data[i][j] = 0.0;
}
return matrix;
}
void freeMatrix(Matrix *matrix)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 20:06
释放矩阵内存
用法 freeMatrix(&matrix)
Changes :
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
for (int i = 0; i < matrix->rows; i++)
{
free(matrix->data[i]);
}
free(matrix->data);
matrix->data = NULL;
matrix->rows = 0;
matrix->cols = 0;
}
Vector vecSet(int size, double number)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 20:07
设置向量
--------------------------------------------------------
Input Parameters :
size ----向量长度
number ----要设置的double型浮点数
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
Vector vector = createVector(size);
int i;
for (i = 0; i < size; i++)
vector.data[i] = number;
return vector;
}
Vector vecCopy(Vector v1)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 20:08
向量复制
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
Vector v2 = createVector(v1.size);
int i;
for (i = 0; i < v1.size; i++)
{
v2.data[i] = v1.data[i];
}
return v2;
}
Vector vecAddNum(Vector v1, double number)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 20:08
向量每个元素加一个数值
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
Vector v2 = createVector(v1.size);
int i;
for (i = 0; i < v1.size; i++)
{
v2.data[i] = v1.data[i] + number;
}
return v2;
}
Vector vecMinusNum(Vector v1, double number)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 20:09
向量每个元素减一个数值
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
Vector v2 = createVector(v1.size);
int i;
for (i = 0; i < v1.size; i++)
{
v2.data[i] = v1.data[i] - number;
}
return v2;
}
Vector vecMulNum(Vector v1, double number)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 20:09
向量每个元素乘以一个数值
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
Vector v2 = createVector(v1.size);
int i;
for (i = 0; i < v1.size; i++)
{
v2.data[i] = v1.data[i] * number;
}
return v2;
}
Vector vecDivNum(Vector v1, double number)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 20:09
向量每个元素除以一个数值
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
Vector v2 = createVector(v1.size);
int i;
for (i = 0; i < v1.size; i++)
{
v2.data[i] = v1.data[i] / number;
}
return v2;
}
double vectorDot(Vector v1, Vector v2)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 20:09
两个向量内积
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
double sum = 0;
int i;
for (i = 0; i < v1.size; i++)
{
sum = sum + v1.data[i] * v2.data[i];
}
return sum;
}
double vectorNorm(Vector v1)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 20:10
向量长度
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
double sum = vectorDot(v1, v1);
return sqrt(sum);
}
// Vector vecAddVec();
// Vector vecMinusVec();
// Vector vecDot();
// vector vecCross();
void matMulVec(Matrix matrix, Vector vector, Vector *result)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 20:10
矩阵乘以向量
--------------------------------------------------------
Input Parameters :
matrix
vector
Output Parameters :
result
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
//Vector result = createVector(matrix.rows);
for (int i = 0; i < matrix.rows; i++)
{
result->data[i] = 0;
for (int j = 0; j < matrix.cols; j++)
{
result->data[i] += matrix.data[i][j] * vector.data[j];
}
}
}
void matMulMat(Matrix A, Matrix B, Matrix *C)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 20:10
矩阵乘以矩阵
用法 matMulMat(A,B,&C)
--------------------------------------------------------
Input Parameters :
A--------
B
Output Parameters :
C
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
if (A.cols != B.rows)
{
printf(
"Error: Number of columns in A must be equal to number of rows in B for multiplication.\n");
exit(1);
}
//Matrix C = createMatrix(A.rows, B.cols);
for (int i = 0; i < A.rows; i++)
{
for (int j = 0; j < B.cols; j++)
{
C->data[i][j] = 0;
for (int k = 0; k < A.cols; k++)
{
C->data[i][j] += A.data[i][k] * B.data[k][j];
}
}
}
}
void matrixTranspose(Matrix A, Matrix *AT)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 16:09
矩阵转置
--------------------------------------------------------
Input Parameters :
A
Output Parameters :
AT = A'
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
if (A.rows != AT->cols || A.cols != AT->rows)
{
printf("Incorrect matrix dimensions");
}
int i, j;
for (i = 0; i < A.cols; i++)
for (j = 0; j < A.rows; j++)
{
AT->data[i][j] = A.data[j][i];
}
}
void AddNormalEqation(Matrix *ATA, Vector *ATb, Vector H, double oc)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 13:51
add a new observation to normal equations
(ATA+ H'H) x = ATb + H'oc
ATA = ATA + H'H
ATb = ATb + H'oc
--------------------------------------------------------
Input Parameters :
ATA
ATb
H
oc
Output Parameters :
ATA -- updated
ATb -- updated
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
int N = ATA->rows;
int i, j;
for (i = 0; i < N; i++)
{
for (j = 0; j < N; j++)
{
ATA->data[i][j] = ATA->data[i][j] + H.data[i] * H.data[j];
}
ATb->data[i] = ATb->data[i] + H.data[i] * oc;
}
}
void LeastSquareEquation(Matrix ATA, Vector ATb, Vector *X, Matrix *P)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 14:04
法方程求解最小二乘
--------------------------------------------------------
Input Parameters :
ATA
ATb
Output Parameters :
X --- 方程解
P ---- 协因素 inv(ATA)
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
int N = ATA.rows;
Matrix L = createMatrix(N, N);
Vector D = createVector(N);
LDL(ATA, &L, &D);
double tmp1;
int i, k;
Vector y = createVector(N);
y.data[0] = ATb.data[0];
for (i = 1; i < N; i++)
{
tmp1 = 0.0;
for (k = 0; k <= i - 1; k++)
{
tmp1 = tmp1 + L.data[i][k] * y.data[k];
}
y.data[i] = ATb.data[i] - tmp1;
}
X->data[N - 1] = y.data[N - 1] / D.data[N - 1];
//----------------------
for (i = N - 1; i >= 0; i--)
{
tmp1 = 0.0;
for (k = i + 1; k < N; k++)
{
tmp1 = tmp1 + L.data[k][i] * X->data[k];
}
X->data[i] = y.data[i] / D.data[i] - tmp1;
}
covariance(L, D, P);
freeMatrix(&L);
freeVector(&D);
freeVector(&y);
}
void LDL(Matrix A, Matrix *L, Vector *D)
/*------------------------------------------ comment
Author : 宋叶志
Date : 2024.02.24
---------------------------------------------------
Desciption : 不开平方的CHOLESKY分解
* LDL'分解
Post Script :
*
input :
* A----------对称正定矩阵
output :
* L----------L矩阵
* D----------对角矩阵-（这里实际用N维向量存储）
*
-------------------------------------------------*/
{
int N = A.rows;
Matrix G = createMatrix(N, N);
int i, j, k;
//设置初值
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
{
L->data[i][j] = 0.0;
G.data[i][j] = 0.0;
}
//D(0)=A(0,0);
D->data[0] = A.data[0][0];
double tmp1;
for (i = 1; i < N; i++)
{
for (j = 0; j <= i - 1; j++)
//此层循环计算g(i,j)
{
tmp1 = 0.0;
for (k = 0; k <= j - 1; k++)
//tmp1=tmp1+G(i,k)*L(j,k);
tmp1 = tmp1 + G.data[i][k] * L->data[j][k];
//G(i,j)=A(i,j)-tmp1;
G.data[i][j] = A.data[i][j] - tmp1;
}
//此循环计算L(i,j)
for (j = 0; j <= i - 1; j++)
//L(i,j)=G(i,j)/D(j);
L->data[i][j] = G.data[i][j] / D->data[j];
//以下计算D(i)
tmp1 = 0.0;
for (k = 0; k <= i - 1; k++)
//tmp1=tmp1+G(i,k)*L(i,k);
tmp1 = tmp1 + G.data[i][k] * L->data[i][k];
//D(i)=A(i,i)-tmp1;
D->data[i] = A.data[i][i] - tmp1;
}
//设置对角元素
for (i = 0; i < N; i++)
//L(i, i) = 1.0;
L->data[i][i] = 1.0;
freeMatrix(&G);
}
void covariance(Matrix L, Vector D, Matrix *P)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 17:03
根据LDL分解计算矩阵协方差
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
int N = L.rows;
Matrix R = createMatrix(N, N);
inversLowTriMatrix(L, &R);
Matrix iU = createMatrix(N, N);
matrixTranspose(R, &iU);
int i, j;
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
{
iU.data[i][j] = iU.data[i][j] / D.data[j];
}
matMulMat(iU, R, P);
freeMatrix(&R);
freeMatrix(&iU);
}
void inversLowTriMatrix(Matrix L, Matrix *R)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 15:59
Inverse of Lower triangular matrix
下三角矩阵的逆矩阵
Changes :
*
a11 b11
a21 a22 b21 b22
a31 a32 a33 b31 b32 b33
a41 a42 a43 a44 b41 b42 b43 b44
a51 a52 a53 a54 a55 b51 b52 b53 b54 b55
c11
c21 c22
c31 c32 c33
c41 c42 c43 c44
c51 c52 c53 c54 c55
0+a(4,2)*b(2,2)+a(4,3)*b(3,2)+a(4,4)*b(4,2)+0=0
a(4,4)*b(4,2)= a(4,k)*b(k,2) k=2,4
--------------------------------------------------------
Input Para :
L--- lower triangula matrix
Output Para :
R---inverse of L(also lower triangula matrix)
Ref:
1.
Chapter5,p228
2.
P58.
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
int M = L.rows;
int i, j, k;
double sum;
for (i = 0; i < M; i++)
{
R->data[i][i] = 1.0 / L.data[i][i];
}
for (i = 0; i < M; i++)
{
for (j = 0; j < M; j++)
{
sum = 0.0;
for (k = i; k < j; k++)
{
sum = sum + L.data[j][k] * R->data[k][i];
R->data[j][i] = -sum * R->data[j][j];
}
}
}
}
void inversUperTriMatrix(Matrix U, Matrix *R)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 16:31
Inverse of Uper triangular matrix
上三角矩阵的逆矩阵
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
int N = U.rows;
Matrix UT = createMatrix(N, N);
matrixTranspose(U, &UT);
Matrix R1 = createMatrix(N, N);
inversLowTriMatrix(UT, &R1);
//======================
matrixTranspose(R1, R);
//======================
freeMatrix(&R1);
freeMatrix(&UT);
}
void inversSymmetricPositiveMatrix(Matrix A, Matrix *invA)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 16:43
invers of Symmetric positive-definite matrix
对称正定矩阵求逆
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
int N = A.rows;
Matrix L = createMatrix(N, N);
Vector D = createVector(N);
LDL(A, &L, &D);
Matrix R = createMatrix(N, N);
inversLowTriMatrix(L, &R);
Matrix iU = createMatrix(N, N);
matrixTranspose(R, &iU);
int i, j;
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
{
iU.data[i][j] = iU.data[i][j] / D.data[j];
}
matMulMat(iU, R, invA);
freeMatrix(&L);
freeVector(&D);
freeMatrix(&R);
freeMatrix(&iU);
}
void uptriEqation(Matrix A, Vector b, Vector *x)
/*------------------------------------------ comment
Author : 宋叶志
Date : 2024.02.25
---------------------------------------------------
Desciption :
* 上三角矩阵方程解法
A * x = b
Post Script :
*
input :
* A-----系数矩阵
* b-----右向量
output :
* x-----方程的解
-------------------------------------------------*/
{
int N;
N = A.rows;
//获取方程维数
//x(N - 1) = b(N - 1) / A(N - 1, N - 1);
x->data[N - 1] = b.data[N - 1] / A.data[N - 1][N - 1];
for (int i = N - 2; i >= 0; i--)
{
//x(i) = b(i);
x->data[i] = b.data[i];
for (int j = i + 1; j <= N - 1; j++)
//x(i) = x(i) - A(i, j) * x(j);
x->data[i] = x->data[i] - A.data[i][j] * x->data[j];
//x(i) = x(i) / A(i, i);
x->data[i] = x->data[i] / A.data[i][i];
}
}
void downtriEqation(Matrix A, Vector b, Vector *x)
/*------------------------------------------ comment
Author : 宋叶志
Date : 2024.02.25
Ref: song
---------------------------------------------------
Desciption :
* 下三角矩阵方程
A * x = b
Post Script :
* 要用到 矩阵类MAT及向量类 VEC
input parameters :
* A----方程系数
* b----右向量
output parameters :
* x----计算结果
-------------------------------------------------*/
{
int N;
N = A.rows;
//x(0) = b(0) / A(0, 0);
x->data[0] = b.data[0] / A.data[0][0];
for (int i = 1; i < N; i++)
{
//x(i) = b(i);
x->data[i] = b.data[i];
//逐元素减
for (int j = 0; j <= i - 1; j++)
x->data[i] = x->data[i] - A.data[i][j] * x->data[j];
//x(i) = x(i) - A(i, j) * x(j);
//x(i) = x(i) / A(i, i);
x->data[i] = x->data[i] / A.data[i][i];
}
}
void printMatrix(Matrix matrix, char format)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 1:25
简易矩阵输出屏幕函数
Changes :
--------------------------------------------------------
Input Parameters :
matrix
format 1. l或L 小数表示法， long format
2. s或S 小数表示法 short format
3. e 科学计数法 long format
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
for (int i = 0; i < matrix.rows; i++)
{
for (int j = 0; j < matrix.cols; j++)
{
switch (format)
{
case 's':
case 'S':
printf("%18.4f ", matrix.data[i][j]);
break;
case 'l':
case 'L':
printf("%18.10f ", matrix.data[i][j]);
break;
case 'e':
case 'E':
printf("%18e ", matrix.data[i][j]);
break;
default:
printf("wrong format");
}
}
printf("\n");
}
}
void printVector(Vector vector, char format)
/*------------------------------------------------------
Created : Song Yezhi 2024-2-25 1:48
简易向量输出函数
Changes :
--------------------------------------------------------
Input Parameters :
vector
format 1. l或L 小数表示法， long format
2. s或S 小数表示法 short format
3. e 科学计数法 long format
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2024
-------------------------------------------------------*/
{
if (vector.size <= 6)
{
for (int i = 0; i < vector.size; i++)
{
switch (format)
{
case 's':
case 'S':
printf("%18.4f ", vector.data[i]);
break;
case 'l':
case 'L':
printf("%18.10f ", vector.data[i]);
break;
case 'e':
case 'E':
printf("%18e ", vector.data[i]);
break;
default:
printf("wrong format");
}
printf("\n");
}
}
if (vector.size > 6)
{
for (int i = 0; i < vector.size; i++)
{
switch (format)
{
case 's':
case 'S':
printf("%18.4f \n ", vector.data[i]);
break;
case 'l':
case 'L':
printf("%18.10f \n ", vector.data[i]);
break;
case 'e':
case 'E':
printf("%18e \n ", vector.data[i]);
break;
default:
printf("wrong format");
}
}
}
}
void rotMatrix(Matrix *Rmat, double angle, char axID)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2024.02.25
rotations matrix
--------------------------------------------------------
Input Parameters :
angle----input angle (rad)
char axID----the ax to of rotation
Output Parameters :
Rmat----rotation matrix
Notes :
** it is very important for the break in each switch case
in C and C++ language
--------------------------------------------------------
Author : 宋叶志
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
double C = cos(angle);
double S = sin(angle);
switch (axID)
{
case 'x':
Rmat->data[0][0] = 1.0;
Rmat->data[0][1] = 0.0;
Rmat->data[0][2] = 0.0;
Rmat->data[1][0] = 0.0;
Rmat->data[1][1] = +C;
Rmat->data[1][2] = +S;
Rmat->data[2][0] = 0.0;
Rmat->data[2][1] = -S;
Rmat->data[2][2] = +C;
break; //it is very important for the "break" at every case
case 'y':
Rmat->data[0][0] = +C;
Rmat->data[0][1] = 0.0;
Rmat->data[0][2] = -S;
Rmat->data[1][0] = 0.0;
Rmat->data[1][1] = 1.0;
Rmat->data[1][2] = 0.0;
Rmat->data[2][0] = +S;
Rmat->data[2][1] = 0.0;
Rmat->data[2][2] = +C;
break;
case 'z':
Rmat->data[0][0] = +C;
Rmat->data[0][1] = +S;
Rmat->data[0][2] = 0.0;
Rmat->data[1][0] = -S;
Rmat->data[1][1] = +C;
Rmat->data[1][2] = 0.0;
Rmat->data[2][0] = 0.0;
Rmat->data[2][1] = 0.0;
Rmat->data[2][2] = 1.0;
break;
default:
printf("the wrong rotation ax \n");
break;
}
}
void interpLagrange(Vector x, Vector y, Vector xx, Vector *yy)
/*------------------------------------------------------
Versions and Changes :
v1.0----2024.02.25
lagrange 插值方法， 对多点一次性完成插值
C语言中，允许 xx为 Vector(1) yy为 Vector(1)
--------------------------------------------------------
paramtegers :
* x-----------样本数据
* y-----------样本函数值
* xx----------要计算的点
* 返回值为-----插值结果，维数同xx
Ref . :
scientific computation in C# Tsinghua Univ. Press
------------------------------------------------------
Author : 宋叶志
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
//获得样本维数
int N = x.size;
//获得被插值点的数据
int M = xx.size;
int p, i;
for (p = 0; p < M; p++)
//对所有的要插值的点做循环
{
//至各分式之和为0
double sum = 0;
for (i = 0; i < N; i++)
{
double tmp = y.data[i];
for (int k = 0; k < N; k++)
{
if (k == i)
continue;
//如果k与i相等，则退出当次循环
//但并不退出循环体
//tmp = tmp * (xx(p) - x(k)) / (x(i) - x(k));
tmp = tmp * (xx.data[p] - x.data[k]) / (x.data[i] - x.data[k]);
}
sum = sum + tmp;
}
//结果赋值给要计算的点
yy->data[p] = sum;
}
}
#include
#include
#include "matlib.h"
#include
void test_LeastSquareEquation(FILE *inputFile, FILE *outputFile)
{
char line[200];
char tmp[20];
// Matrix A;
// Vector b;
// Vector X;
int rows, cols, i, j;
while (1)
{
if (feof(inputFile))
break;
fgets(line, sizeof(line), inputFile);
if (strncmp(line, "LeastSqaureEquation", 19) == 0)
{
sscanf(line, "%20s%d%d", tmp, &rows, &cols);
break;
}
}
Matrix ATA = createMatrix(cols, cols);
Vector ATb = createVector(cols);
Vector H = createVector(cols);
double oc;
fgets(line, sizeof(line), inputFile);
//跳过一行注释
for (i = 0; i < rows; i++)
{
for (j = 0; j < cols; j++)
{
fscanf(inputFile, "%lf", &H.data[j]);
}
fscanf(inputFile, "%lf", &oc);
// printf("H = \n");
// printVector(H,'l');
// printf("oc = %lf \n",oc);
//累加法方程
AddNormalEqation(&ATA, &ATb, H, oc);
// printf("ATA = \n ");
// printMatrix(ATA,'l');
// printf("ATb = \n ");
// printVector(ATb,'l');
}
printf("ATA = \n");
printMatrix(ATA, 'l');
printf("ATb = \n");
printVector(ATb, 'l');
Vector X = createVector(cols);
Matrix P = createMatrix(cols, cols);
LeastSquareEquation(ATA, ATb, &X, &P);
printf("X = \n");
printVector(X, 'l');
printf("Cov = \n");
printMatrix(P, 'l');
inversSymmetricPositiveMatrix(ATA, &P);
printf("invATA = \n");
printMatrix(P, 'l');
freeMatrix(&ATA);
freeVector(&ATb);
freeVector(&H);
freeVector(&X);
freeMatrix(&P);
}
int main()
{
FILE *inputFile;
FILE *outputFile;
inputFile = fopen("testdata.txt", "r");
outputFile = fopen("outputdata.txt", "w");
//test_vecSet();
//test_vecCopy();
//test_vecAddnum( );
//test_vecMinusNum( );
//test_vecMulNum( );
//test_vecDivNum( );
//test_matMulVec();
//test_matMulMat();
//test_vectorDot();
//test_LeastSquareEquation( inputFile, outputFile);
//test_LDL( inputFile, outputFile);
//test_interp();
fclose(inputFile);
fclose(outputFile);
return 0;
}
linear algebra test data sets
##################################################################
LDL 6
10.743651420771521 6.665748929147245 7.324667285683409 6.737351580871993 5.753170110245249 6.478842535949054
6.665748929147245 8.102999126026232 6.559760228773046 6.146728226995120 5.434970979629590 5.461718378432780
7.324667285683409 6.559760228773046 7.737843465942078 6.181422895127592 5.273400381539975 6.044952695243294
6.737351580871993 6.146728226995120 6.181422895127592 7.115394200317208 5.494739962412434 5.164683353703877
5.753170110245249 5.434970979629590 5.273400381539975 5.494739962412434 5.866541080893321 5.374679283386545
6.478842535949054 5.461718378432780 6.044952695243294 5.164683353703877 5.374679283386545 8.135734690511807
MATLAB ANS LDL
l =
1.000000000000000 0 0 0 0 0
0.620436075975049 1.000000000000000 0 0 0 0
0.681767026759829 0.507967173772834 1.000000000000000 0 0 0
0.627100723674463 0.495706995053581 0.342433308901797 1.000000000000000 0 0
0.535494859701256 0.470214886196089 0.234514940582084 0.480820502509731 1.000000000000000 0
0.603039160729193 0.363471518460910 0.520448357440642 0.046894962098067 0.690262318715572 1.000000000000000
d =
10.743651420771521 0 0 0 0 0
0 3.967328016991228 0 0 0 0
0 0 1.720434603054151 0 0 0
0 0 0 1.713783641431062 0 0
0 0 0 0 1.417737738533325 0
0 0 0 0 0 2.559333903885087
###################################################################
LeastSqaureEquation 18 5
[A(18,5) ,b(18,1)]
0.181847028302852 0.417267069084370 0.575208595078466 0.368484596490336 0.350727103576883 5.969580855552750
0.263802916521990 0.049654430325742 0.059779542947156 0.625618560729690 0.939001561999887 7.739932458933137
0.145538980384717 0.902716109915281 0.234779913372406 0.780227435151377 0.875942811492984 10.155934738402925
0.136068558708664 0.944787189721646 0.353158571222071 0.081125768865785 0.550156342898422 6.160403441773421
0.869292207640089 0.490864092468080 0.821194040197959 0.929385970968730 0.622475086001227 11.144521827051184
0.579704587365570 0.489252638400019 0.015403437651555 0.775712678608402 0.587044704531417 7.642494414210966
0.549860201836332 0.337719409821377 0.043023801657808 0.486791632403172 0.207742292733028 4.340248419730341
0.144954798223727 0.900053846417662 0.168990029462704 0.435858588580919 0.301246330279491 5.701698585168293
0.853031117721894 0.369246781120215 0.649115474956452 0.446783749429806 0.470923348517591 7.680622845138857
0.622055131485066 0.111202755293787 0.731722385658670 0.306349472016557 0.230488160211558 5.417466488172673
0.350952380892271 0.780252068321138 0.647745963136307 0.508508655381127 0.844308792695389 10.110272991944919
0.513249539867053 0.389738836961253 0.450923706430945 0.510771564172110 0.194764289567049 5.662406037606080
0.401808033751942 0.241691285913833 0.547008892286345 0.817627708322262 0.225921780972399 6.926337020589685
0.075966691690842 0.403912145588115 0.296320805607773 0.794831416883453 0.170708047147859 5.805619302963495
0.239916153553658 0.096454525168389 0.744692807074156 0.644318130193692 0.227664297816554 6.382497634970438
0.123318934835166 0.131973292606335 0.188955015032545 0.378609382660268 0.435698684103899 4.647061516306039
0.183907788282417 0.942050590775485 0.686775433365315 0.811580458282477 0.311102286650413 8.930168536311305
0.239952525664903 0.956134540229802 0.183511155737270 0.532825588799455 0.923379642103244 9.450955639050354
ATA =
3.390035635975875 2.803314241005550 3.055791110742465 3.824672558330794 2.985394845803748
2.803314241005550 6.247895194369039 3.539144538853189 5.057853433992208 4.750703394200390
3.055791110742465 3.539144538853189 4.281595235000124 4.294108273380072 3.156080646448181
3.824672558330794 5.057853433992208 4.294108273380072 6.649206245407066 4.897452377637229
2.985394845803748 4.750703394200390 3.156080646448181 4.897452377637229 5.237629093375518
ATb =
48.389701912556291
69.901468953273991
55.935702219210405
77.906791116269844
67.732998550975566
inv(ATA) =
1.130272617301131 0.211862852877081 -0.478366199297830 -0.316818608528521 -0.251915751507541
0.211862852877081 0.626277384496548 -0.197373964767773 -0.163982049998924 -0.416549147594521
-0.478366199297830 -0.197373964767773 0.885191220739288 -0.276863312760076 0.177173133639693
-0.316818608528521 -0.163982049998924 -0.276863312760076 0.869693512585635 -0.317055840693516
-0.251915751507541 -0.416549147594521 0.177173133639693 -0.317055840693516 0.902042228253418
The ans of LeastSqaureEquation is [1 2 3 4 5]