在C中实现矩阵运算

               

本文博客链接:http://blog.youkuaiyun.com/jdh99,作者:jdh,转载请注明.

环境：

主机:XP

开发环境:mingw

功能:

在C++中实现矩阵运算

说明：

将这篇文章(http://blog.youkuaiyun.com/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  

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; }

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;}