本文博客链接:http://blog.csdn.net/jdh99,作者:jdh
环境:
主机:XP
开发环境:mingw
功能:
在C++中实现矩阵运算
说明:
将这篇文章(http://blog.csdn.net/jdh99/article/details/7360091)中在C++下实现的矩阵运算移植到C下,以便在单片机中运算.
源代码:
matrix.h:
- #ifndef _MATRIX_H
- #define _MATRIX_H
- //头文件
- #include <stdio.h>
- #include <stdlib.h>
- //矩阵数据结构
- //二维矩阵
- struct _Matrix
- {
- int m;
- int n;
- float *arr;
- };
- //矩阵方法
- //设置m
- void matrix_set_m(struct _Matrix *m,int mm);
- //设置n
- void matrix_set_n(struct _Matrix *m,int nn);
- //初始化
- void matrix_init(struct _Matrix *m);
- //释放
- void matrix_free(struct _Matrix *m);
- //读取i,j坐标的数据
- //失败返回-31415,成功返回值
- float matrix_read(struct _Matrix *m,int i,int j);
- //写入i,j坐标的数据
- //失败返回-1,成功返回1
- int matrix_write(struct _Matrix *m,int i,int j,float val);
- //矩阵运算
- //成功返回1,失败返回-1
- int matrix_add(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C);
- //C = A - B
- //成功返回1,失败返回-1
- int matrix_subtract(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C);
- //C = A * B
- //成功返回1,失败返回-1
- int matrix_multiply(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C);
- //行列式的值,只能计算2 * 2,3 * 3
- //失败返回-31415,成功返回值
- float matrix_det(struct _Matrix *A);
- //求转置矩阵,B = AT
- //成功返回1,失败返回-1
- int matrix_transpos(struct _Matrix *A,struct _Matrix *B);
- //求逆矩阵,B = A^(-1)
- //成功返回1,失败返回-1
- int matrix_inverse(struct _Matrix *A,struct _Matrix *B);
- #endif
#ifndef _MATRIX_H
#define _MATRIX_H
//头文件
#include <stdio.h>
#include <stdlib.h>
//矩阵数据结构
//二维矩阵
struct _Matrix
{
int m;
int n;
float *arr;
};
//矩阵方法
//设置m
void matrix_set_m(struct _Matrix *m,int mm);
//设置n
void matrix_set_n(struct _Matrix *m,int nn);
//初始化
void matrix_init(struct _Matrix *m);
//释放
void matrix_free(struct _Matrix *m);
//读取i,j坐标的数据
//失败返回-31415,成功返回值
float matrix_read(struct _Matrix *m,int i,int j);
//写入i,j坐标的数据
//失败返回-1,成功返回1
int matrix_write(struct _Matrix *m,int i,int j,float val);
//矩阵运算
//成功返回1,失败返回-1
int matrix_add(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C);
//C = A - B
//成功返回1,失败返回-1
int matrix_subtract(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C);
//C = A * B
//成功返回1,失败返回-1
int matrix_multiply(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C);
//行列式的值,只能计算2 * 2,3 * 3
//失败返回-31415,成功返回值
float matrix_det(struct _Matrix *A);
//求转置矩阵,B = AT
//成功返回1,失败返回-1
int matrix_transpos(struct _Matrix *A,struct _Matrix *B);
//求逆矩阵,B = A^(-1)
//成功返回1,失败返回-1
int matrix_inverse(struct _Matrix *A,struct _Matrix *B);
#endif
matrix.c:
- #include "matrix.h"
- //矩阵方法
- //设置m
- void matrix_set_m(struct _Matrix *m,int mm)
- {
- m->m = mm;
- }
- //设置n
- void matrix_set_n(struct _Matrix *m,int nn)
- {
- m->n = nn;
- }
- //初始化
- void matrix_init(struct _Matrix *m)
- {
- m->arr = (float *)malloc(m->m * m->n * sizeof(float));
- }
- //释放
- void matrix_free(struct _Matrix *m)
- {
- free(m->arr);
- }
- //读取i,j坐标的数据
- //失败返回-31415,成功返回值
- float matrix_read(struct _Matrix *m,int i,int j)
- {
- if (i >= m->m || j >= m->n)
- {
- return -31415;
- }
- return *(m->arr + i * m->n + j);
- }
- //写入i,j坐标的数据
- //失败返回-1,成功返回1
- int matrix_write(struct _Matrix *m,int i,int j,float val)
- {
- if (i >= m->m || j >= m->n)
- {
- return -1;
- }
- *(m->arr + i * m->n + j) = val;
- return 1;
- }
- //矩阵运算
- //成功返回1,失败返回-1
- int matrix_add(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C)
- {
- int i = 0;
- int j = 0;
- //判断是否可以运算
- if (A->m != B->m || A->n != B->n || \
- A->m != C->m || A->n != C->n)
- {
- return -1;
- }
- //运算
- for (i = 0;i < C->m;i++)
- {
- for (j = 0;j < C->n;j++)
- {
- matrix_write(C,i,j,matrix_read(A,i,j) + matrix_read(B,i,j));
- }
- }
- return 1;
- }
- //C = A - B
- //成功返回1,失败返回-1
- int matrix_subtract(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C)
- {
- int i = 0;
- int j = 0;
- //判断是否可以运算
- if (A->m != B->m || A->n != B->n || \
- A->m != C->m || A->n != C->n)
- {
- return -1;
- }
- //运算
- for (i = 0;i < C->m;i++)
- {
- for (j = 0;j < C->n;j++)
- {
- matrix_write(C,i,j,matrix_read(A,i,j) - matrix_read(B,i,j));
- }
- }
- return 1;
- }
- //C = A * B
- //成功返回1,失败返回-1
- int matrix_multiply(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C)
- {
- int i = 0;
- int j = 0;
- int k = 0;
- float temp = 0;
- //判断是否可以运算
- if (A->m != C->m || B->n != C->n || \
- A->n != B->m)
- {
- return -1;
- }
- //运算
- for (i = 0;i < C->m;i++)
- {
- for (j = 0;j < C->n;j++)
- {
- temp = 0;
- for (k = 0;k < A->n;k++)
- {
- temp += matrix_read(A,i,k) * matrix_read(B,k,j);
- }
- matrix_write(C,i,j,temp);
- }
- }
- return 1;
- }
- //行列式的值,只能计算2 * 2,3 * 3
- //失败返回-31415,成功返回值
- float matrix_det(struct _Matrix *A)
- {
- float value = 0;
- //判断是否可以运算
- if (A->m != A->n || (A->m != 2 && A->m != 3))
- {
- return -31415;
- }
- //运算
- if (A->m == 2)
- {
- value = matrix_read(A,0,0) * matrix_read(A,1,1) - matrix_read(A,0,1) * matrix_read(A,1,0);
- }
- else
- {
- value = matrix_read(A,0,0) * matrix_read(A,1,1) * matrix_read(A,2,2) + \
- matrix_read(A,0,1) * matrix_read(A,1,2) * matrix_read(A,2,0) + \
- matrix_read(A,0,2) * matrix_read(A,1,0) * matrix_read(A,2,1) - \
- matrix_read(A,0,0) * matrix_read(A,1,2) * matrix_read(A,2,1) - \
- matrix_read(A,0,1) * matrix_read(A,1,0) * matrix_read(A,2,2) - \
- matrix_read(A,0,2) * matrix_read(A,1,1) * matrix_read(A,2,0);
- }
- return value;
- }
- //求转置矩阵,B = AT
- //成功返回1,失败返回-1
- int matrix_transpos(struct _Matrix *A,struct _Matrix *B)
- {
- int i = 0;
- int j = 0;
- //判断是否可以运算
- if (A->m != B->n || A->n != B->m)
- {
- return -1;
- }
- //运算
- for (i = 0;i < B->m;i++)
- {
- for (j = 0;j < B->n;j++)
- {
- matrix_write(B,i,j,matrix_read(A,j,i));
- }
- }
- return 1;
- }
- //求逆矩阵,B = A^(-1)
- //成功返回1,失败返回-1
- int matrix_inverse(struct _Matrix *A,struct _Matrix *B)
- {
- int i = 0;
- int j = 0;
- int k = 0;
- struct _Matrix m;
- float temp = 0;
- float b = 0;
- //判断是否可以运算
- if (A->m != A->n || B->m != B->n || A->m != B->m)
- {
- return -1;
- }
- /*
- //如果是2维或者3维求行列式判断是否可逆
- if (A->m == 2 || A->m == 3)
- {
- if (det(A) == 0)
- {
- return -1;
- }
- }
- */
- //增广矩阵m = A | B初始化
- matrix_set_m(&m,A->m);
- matrix_set_n(&m,2 * A->m);
- matrix_init(&m);
- for (i = 0;i < m.m;i++)
- {
- for (j = 0;j < m.n;j++)
- {
- if (j <= A->n - 1)
- {
- matrix_write(&m,i,j,matrix_read(A,i,j));
- }
- else
- {
- if (i == j - A->n)
- {
- matrix_write(&m,i,j,1);
- }
- else
- {
- matrix_write(&m,i,j,0);
- }
- }
- }
- }
- //高斯消元
- //变换下三角
- for (k = 0;k < m.m - 1;k++)
- {
- //如果坐标为k,k的数为0,则行变换
- if (matrix_read(&m,k,k) == 0)
- {
- for (i = k + 1;i < m.m;i++)
- {
- if (matrix_read(&m,i,k) != 0)
- {
- break;
- }
- }
- if (i >= m.m)
- {
- return -1;
- }
- else
- {
- //交换行
- for (j = 0;j < m.n;j++)
- {
- temp = matrix_read(&m,k,j);
- matrix_write(&m,k,j,matrix_read(&m,k + 1,j));
- matrix_write(&m,k + 1,j,temp);
- }
- }
- }
- //消元
- for (i = k + 1;i < m.m;i++)
- {
- //获得倍数
- b = matrix_read(&m,i,k) / matrix_read(&m,k,k);
- //行变换
- for (j = 0;j < m.n;j++)
- {
- temp = matrix_read(&m,i,j) - b * matrix_read(&m,k,j);
- matrix_write(&m,i,j,temp);
- }
- }
- }
- //变换上三角
- for (k = m.m - 1;k > 0;k--)
- {
- //如果坐标为k,k的数为0,则行变换
- if (matrix_read(&m,k,k) == 0)
- {
- for (i = k + 1;i < m.m;i++)
- {
- if (matrix_read(&m,i,k) != 0)
- {
- break;
- }
- }
- if (i >= m.m)
- {
- return -1;
- }
- else
- {
- //交换行
- for (j = 0;j < m.n;j++)
- {
- temp = matrix_read(&m,k,j);
- matrix_write(&m,k,j,matrix_read(&m,k + 1,j));
- matrix_write(&m,k + 1,j,temp);
- }
- }
- }
- //消元
- for (i = k - 1;i >= 0;i--)
- {
- //获得倍数
- b = matrix_read(&m,i,k) / matrix_read(&m,k,k);
- //行变换
- for (j = 0;j < m.n;j++)
- {
- temp = matrix_read(&m,i,j) - b * matrix_read(&m,k,j);
- matrix_write(&m,i,j,temp);
- }
- }
- }
- //将左边方阵化为单位矩阵
- for (i = 0;i < m.m;i++)
- {
- if (matrix_read(&m,i,i) != 1)
- {
- //获得倍数
- b = 1 / matrix_read(&m,i,i);
- //行变换
- for (j = 0;j < m.n;j++)
- {
- temp = matrix_read(&m,i,j) * b;
- matrix_write(&m,i,j,temp);
- }
- }
- }
- //求得逆矩阵
- for (i = 0;i < B->m;i++)
- {
- for (j = 0;j < B->m;j++)
- {
- matrix_write(B,i,j,matrix_read(&m,i,j + m.m));
- }
- }
- //释放增广矩阵
- matrix_free(&m);
- return 1;
- }
#include "matrix.h"
//矩阵方法
//设置m
void matrix_set_m(struct _Matrix *m,int mm)
{
m->m = mm;
}
//设置n
void matrix_set_n(struct _Matrix *m,int nn)
{
m->n = nn;
}
//初始化
void matrix_init(struct _Matrix *m)
{
m->arr = (float *)malloc(m->m * m->n * sizeof(float));
}
//释放
void matrix_free(struct _Matrix *m)
{
free(m->arr);
}
//读取i,j坐标的数据
//失败返回-31415,成功返回值
float matrix_read(struct _Matrix *m,int i,int j)
{
if (i >= m->m || j >= m->n)
{
return -31415;
}
return *(m->arr + i * m->n + j);
}
//写入i,j坐标的数据
//失败返回-1,成功返回1
int matrix_write(struct _Matrix *m,int i,int j,float val)
{
if (i >= m->m || j >= m->n)
{
return -1;
}
*(m->arr + i * m->n + j) = val;
return 1;
}
//矩阵运算
//成功返回1,失败返回-1
int matrix_add(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C)
{
int i = 0;
int j = 0;
//判断是否可以运算
if (A->m != B->m || A->n != B->n || \
A->m != C->m || A->n != C->n)
{
return -1;
}
//运算
for (i = 0;i < C->m;i++)
{
for (j = 0;j < C->n;j++)
{
matrix_write(C,i,j,matrix_read(A,i,j) + matrix_read(B,i,j));
}
}
return 1;
}
//C = A - B
//成功返回1,失败返回-1
int matrix_subtract(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C)
{
int i = 0;
int j = 0;
//判断是否可以运算
if (A->m != B->m || A->n != B->n || \
A->m != C->m || A->n != C->n)
{
return -1;
}
//运算
for (i = 0;i < C->m;i++)
{
for (j = 0;j < C->n;j++)
{
matrix_write(C,i,j,matrix_read(A,i,j) - matrix_read(B,i,j));
}
}
return 1;
}
//C = A * B
//成功返回1,失败返回-1
int matrix_multiply(struct _Matrix *A,struct _Matrix *B,struct _Matrix *C)
{
int i = 0;
int j = 0;
int k = 0;
float temp = 0;
//判断是否可以运算
if (A->m != C->m || B->n != C->n || \
A->n != B->m)
{
return -1;
}
//运算
for (i = 0;i < C->m;i++)
{
for (j = 0;j < C->n;j++)
{
temp = 0;
for (k = 0;k < A->n;k++)
{
temp += matrix_read(A,i,k) * matrix_read(B,k,j);
}
matrix_write(C,i,j,temp);
}
}
return 1;
}
//行列式的值,只能计算2 * 2,3 * 3
//失败返回-31415,成功返回值
float matrix_det(struct _Matrix *A)
{
float value = 0;
//判断是否可以运算
if (A->m != A->n || (A->m != 2 && A->m != 3))
{
return -31415;
}
//运算
if (A->m == 2)
{
value = matrix_read(A,0,0) * matrix_read(A,1,1) - matrix_read(A,0,1) * matrix_read(A,1,0);
}
else
{
value = matrix_read(A,0,0) * matrix_read(A,1,1) * matrix_read(A,2,2) + \
matrix_read(A,0,1) * matrix_read(A,1,2) * matrix_read(A,2,0) + \
matrix_read(A,0,2) * matrix_read(A,1,0) * matrix_read(A,2,1) - \
matrix_read(A,0,0) * matrix_read(A,1,2) * matrix_read(A,2,1) - \
matrix_read(A,0,1) * matrix_read(A,1,0) * matrix_read(A,2,2) - \
matrix_read(A,0,2) * matrix_read(A,1,1) * matrix_read(A,2,0);
}
return value;
}
//求转置矩阵,B = AT
//成功返回1,失败返回-1
int matrix_transpos(struct _Matrix *A,struct _Matrix *B)
{
int i = 0;
int j = 0;
//判断是否可以运算
if (A->m != B->n || A->n != B->m)
{
return -1;
}
//运算
for (i = 0;i < B->m;i++)
{
for (j = 0;j < B->n;j++)
{
matrix_write(B,i,j,matrix_read(A,j,i));
}
}
return 1;
}
//求逆矩阵,B = A^(-1)
//成功返回1,失败返回-1
int matrix_inverse(struct _Matrix *A,struct _Matrix *B)
{
int i = 0;
int j = 0;
int k = 0;
struct _Matrix m;
float temp = 0;
float b = 0;
//判断是否可以运算
if (A->m != A->n || B->m != B->n || A->m != B->m)
{
return -1;
}
/*
//如果是2维或者3维求行列式判断是否可逆
if (A->m == 2 || A->m == 3)
{
if (det(A) == 0)
{
return -1;
}
}
*/
//增广矩阵m = A | B初始化
matrix_set_m(&m,A->m);
matrix_set_n(&m,2 * A->m);
matrix_init(&m);
for (i = 0;i < m.m;i++)
{
for (j = 0;j < m.n;j++)
{
if (j <= A->n - 1)
{
matrix_write(&m,i,j,matrix_read(A,i,j));
}
else
{
if (i == j - A->n)
{
matrix_write(&m,i,j,1);
}
else
{
matrix_write(&m,i,j,0);
}
}
}
}
//高斯消元
//变换下三角
for (k = 0;k < m.m - 1;k++)
{
//如果坐标为k,k的数为0,则行变换
if (matrix_read(&m,k,k) == 0)
{
for (i = k + 1;i < m.m;i++)
{
if (matrix_read(&m,i,k) != 0)
{
break;
}
}
if (i >= m.m)
{
return -1;
}
else
{
//交换行
for (j = 0;j < m.n;j++)
{
temp = matrix_read(&m,k,j);
matrix_write(&m,k,j,matrix_read(&m,k + 1,j));
matrix_write(&m,k + 1,j,temp);
}
}
}
//消元
for (i = k + 1;i < m.m;i++)
{
//获得倍数
b = matrix_read(&m,i,k) / matrix_read(&m,k,k);
//行变换
for (j = 0;j < m.n;j++)
{
temp = matrix_read(&m,i,j) - b * matrix_read(&m,k,j);
matrix_write(&m,i,j,temp);
}
}
}
//变换上三角
for (k = m.m - 1;k > 0;k--)
{
//如果坐标为k,k的数为0,则行变换
if (matrix_read(&m,k,k) == 0)
{
for (i = k + 1;i < m.m;i++)
{
if (matrix_read(&m,i,k) != 0)
{
break;
}
}
if (i >= m.m)
{
return -1;
}
else
{
//交换行
for (j = 0;j < m.n;j++)
{
temp = matrix_read(&m,k,j);
matrix_write(&m,k,j,matrix_read(&m,k + 1,j));
matrix_write(&m,k + 1,j,temp);
}
}
}
//消元
for (i = k - 1;i >= 0;i--)
{
//获得倍数
b = matrix_read(&m,i,k) / matrix_read(&m,k,k);
//行变换
for (j = 0;j < m.n;j++)
{
temp = matrix_read(&m,i,j) - b * matrix_read(&m,k,j);
matrix_write(&m,i,j,temp);
}
}
}
//将左边方阵化为单位矩阵
for (i = 0;i < m.m;i++)
{
if (matrix_read(&m,i,i) != 1)
{
//获得倍数
b = 1 / matrix_read(&m,i,i);
//行变换
for (j = 0;j < m.n;j++)
{
temp = matrix_read(&m,i,j) * b;
matrix_write(&m,i,j,temp);
}
}
}
//求得逆矩阵
for (i = 0;i < B->m;i++)
{
for (j = 0;j < B->m;j++)
{
matrix_write(B,i,j,matrix_read(&m,i,j + m.m));
}
}
//释放增广矩阵
matrix_free(&m);
return 1;
}
main.c:(测试代码)
- #include <stdio.h>
- #include "matrix.h"
- //打印2维矩阵
- void printf_matrix(struct _Matrix *A)
- {
- int i = 0;
- int j = 0;
- int m = 0;
- int n = 0;
- m = A->m;
- n = A->n;
- for (i = 0;i < m;i++)
- {
- for (j = 0;j < n;j++)
- {
- printf("%f\t",matrix_read(A,i,j));
- }
- printf("\n");
- }
- }
- int main()
- {
- int i = 0;
- int j = 0;
- int k = 0;
- struct _Matrix m1;
- struct _Matrix m2;
- struct _Matrix m3;
- //初始化内存
- matrix_set_m(&m1,3);
- matrix_set_n(&m1,3);
- matrix_init(&m1);
- matrix_set_m(&m2,3);
- matrix_set_n(&m2,3);
- matrix_init(&m2);
- matrix_set_m(&m3,3);
- matrix_set_n(&m3,3);
- matrix_init(&m3);
- //初始化数据
- k = 1;
- for (i = 0;i < m1.m;i++)
- {
- for (j = 0;j < m1.n;j++)
- {
- matrix_write(&m1,i,j,k++);
- }
- }
- for (i = 0;i < m2.m;i++)
- {
- for (j = 0;j < m2.n;j++)
- {
- matrix_write(&m2,i,j,k++);
- }
- }
- //原数据
- printf("A:\n");
- printf_matrix(&m1);
- printf("B:\n");
- printf_matrix(&m2);
- printf("A:行列式的值%f\n",matrix_det(&m1));
- //C = A + B
- if (matrix_add(&m1,&m2,&m3) > 0)
- {
- printf("C = A + B:\n");
- printf_matrix(&m3);
- }
- //C = A - B
- if (matrix_subtract(&m1,&m2,&m3) > 0)
- {
- printf("C = A - B:\n");
- printf_matrix(&m3);
- }
- //C = A * B
- if (matrix_multiply(&m1,&m2,&m3) > 0)
- {
- printf("C = A * B:\n");
- printf_matrix(&m3);
- }
- //C = AT
- if (matrix_transpos(&m1,&m3) > 0)
- {
- printf("C = AT:\n");
- printf_matrix(&m3);
- }
- if (matrix_inverse(&m1,&m3) > 0)
- {
- printf("C = A^(-1):\n");
- printf_matrix(&m3);
- }
- getchar();
- return 0;
- }