基于C语言的矩阵运算库

最近本着锻炼自己编程能力的目的，捣鼓了一下矩阵运算，看到网上很多的矩阵运算库都是基于C++的，很少有基于C语言的，于是自己想要写一个纯C语言的矩阵运算库。大神看到了不要喷我，我只是个小白。
个人感觉矩阵运算最重要的是数据结构的安排，我看到有些使用C语言编写的矩阵运算库里面的矩阵元素是使用以为数组存储的，这样的好处是确定了矩阵的行和列，直接将矩阵结构中的指针指向对应的数组即可；其形式如下：

typedef struct 
{
  int row;
  int column;
  float* data;
}Matrix_t;

这种一维数组有个缺点就是，运算的时候下标不好操作，需要程序员编程功底较好，对于我这种半吊子，还是使用简单点儿的二维数组吧，下标和各个元素一一对应，方便操作。这里使用动态内存分配的方式来管理矩阵的数据存储的分配和释放。

typedef unsigned int uint16_t;
typedef struct
{
    uint16_t row;
    uint16_t column;
    float **data;
}Matrix_t;

首先创建矩阵，创建矩阵时只需要矩阵的行数和列数，矩阵创建函数会自动根据矩阵的行和列分配一片动态内存给矩阵。为了方便起见，将所有的矩阵元素全部置为0；

/*
create a matrix with 0
*/
Matrix_t create_mat(uint16_t row, uint16_t column)
{
    Matrix_t mat;
    if (row <= 0||column<=0)
    {
        printf("error, in create_mat: row <= 0||column<=0\n");
        exit(1);
    }
    if (row > 0 && column > 0)
    {
        mat.row = row;
        mat.column = column;
        mat.data = (float **)malloc(row*sizeof(float *));//先指针的指针
        if (mat.data == NULL)
        {
            printf("error, in create_mat: mat.data==NULL");
            exit(1);
        }
        uint16_t i;
        for (i = 0; i < row; i++)
        {
            *(mat.data + i) = (float *)malloc(column*sizeof(float));//再分配每行的指针
            if (mat.data[i] == NULL)
            {
              printf("error, in create_mat: mat.data==NULL");
              exit(1);
            }
        }
        clear_mat(&mat);
    }
    return mat;
}

矩阵元素置0函数很简单，输入指向对应矩阵的指针，将每个元素设置为0；

/*
set all datas in matrix to zero
*/
void clear_mat(Matrix_t* mat)
{
    uint16_t i, j;
    for (i = 0; i < mat->row; i++)
    {
        for (j = 0; j < mat->column; j++)
        {
            mat->data[i][j] = 0;
        }
    }
}

由于创建矩阵过程中使用malloc函数分配了内存，在矩阵完成运算之后，一定要将矩阵的内存用free释放掉，以防止出现大量内存碎片，使程序崩溃。

/*
free a matrix
*/
void free_mat(Matrix_t *mat)
{
    uint16_t i;
    for (i = 0; i < mat->row; i++)
        free(mat->data[i]);/*释放行*/
    free(mat->data);/*释放头指针*/
}

需要注意的是释放的顺序必须是后申请的内存先释放，否则释放不彻底。
对于矩阵的元素设定需要设计一个函数，该函数将待设置矩阵和一维数组作为输入，将一维数组中的元素一次拷贝到矩阵的二维数组中。之前试过使用指向二维数组的指针，觉得不是很好，还是用一维数组吧。

/*
set datas to the matrix
*/
void set_mat_data(Matrix_t* mat,const float *data)
{
    uint16_t i, j;
    for (i = 0; i < mat->row; i++)
    {
        for (j = 0; j < mat->column; j++)
        {
            mat->data[i][j] = data[i*mat->column+j];
        }
    }
}

矩阵的加减实现方法区别不大，只需对应下标的元素相互加减即可，以下是矩阵相加的情形：

/*
mat=mat1+mat2
*/
Matrix_t add_mat(const Matrix_t* mat1, const Matrix_t* mat2)
{

    if (mat1->row != mat2->row)
    {
        printf("error, in add_mat: mat1->row != mat2->row\n");
        exit(1);
    }
    if (mat1->column != mat2->column)
    {
        printf("error, in add_mat: mat1->column != mat2->column\n");
        exit(1);
    }
    Matrix_t mat;
    uint16_t i, j;
    mat = create_mat(mat1->row, mat1->column);
    for (i = 0; i < mat1->row; i++)
    {
        for (j = 0; j < mat1->column; j++)
            mat.data[i][j] = mat1->data[i][j] + mat2->data[i][j];
    }
    return mat;
}

矩阵转置时将矩阵的行和列对换；

/*
transpose the matrix, mat=mat'
*/
Matrix_t transpose_mat(const Matrix_t* mat)//mat'
{
    Matrix_t mat_T;
    mat_T = create_mat(mat->row, mat->column);
    uint16_t i, j;
    for (i = 0; i < mat->row; i++)
    {
        for (j = 0; j < mat->column; j++)
        {
            mat_T.data[i][j] = mat->data[j][i];
        }
    }
    return mat_T;
}

两个矩阵相乘，根据线性代数中的方法，第一个矩阵的第一行乘以第二个矩阵的每一列，得到结果矩阵的第一行的值。。。。。实现起来并不困难。

*
mat=mat1*mat2
*/
Matrix_t mult_mat(const Matrix_t *mat1, const Matrix_t* mat2)
{
    Matrix_t mat;
    if (mat1->column != mat2->row)
    {
        printf("error,In mult_mat: mat1->column != mat2->row\n");
        exit(1);
    }   
    else
    {
        mat = create_mat(mat1->row, mat2->column);
        clear_mat(&mat);
        uint16_t i, j;
        for (i = 0; i < mat1->row; i++)
        {
            for (j = 0; j < mat2->column; j++)
            {
                uint16_t m;
                for (m = 0; m < mat1->column; m++)
                {
                    mat.data[i][j] += mat1->data[i][m] * mat2->data[m][j];
                }
            }
        }
    }
    return mat;
}

举证行列式的值的计算是比较复杂的，在线性代数中，一般采用的方法是：将矩阵通过行变换（列主元消去法）变换为上三角矩阵，需要特别注意的是，矩阵两行互换时，行列式的值变号，编程时需要考虑。

/*
get matrix's derterminent value
*/
float det_mat(Matrix_t *m)
{
    uint16_t i, j, n, max_row;
    char swap_f;
    float max, k;
    float det=1;
    swap_f = 0;
    if (m->column != m->row)
    {
        printf("error:In det_mat (m->column != m->row)\n");
        exit(1);
    }
    Matrix_t mat = copy_mat(m);
    for (i = 0; i < mat.row-1; i++)
    {
        max = abs(mat.data[i][i]);
        max_row = i;
        for (j = i + 1; j < mat.row; j++)
        {
            if (max < abs(mat.data[j][i]))
            {
                max = abs(mat.data[j][i]);
                max_row = j;
            }
        }
        if (i != max_row)
        {
            swap_row_mat(&mat, i, max_row);
            swap_f++;
        }
        for (j = i + 1; j < mat.row; j++)
        {
            k = -mat.data[j][i]/mat.data[i][i];
            for (n= 0; n < mat.column; n++)
            {
                mat.data[j][n] = mat.data[i][n] * k + mat.data[j][n];
            }
        }
    }
    if (swap_f%2==1)swap_f = -1;
    else swap_f = 1;
    det = 1;
    for (i = 0; i < mat.column; i++)
        det *= mat.data[i][i];
    det *= swap_f;
    //show_mat("in det_mat: mat_upper",&mat);
    free_mat(&mat);
    return det;
}

矩阵的逆矩阵，这是难点所在，其实编程上容易实现的还是Gauss-Jordan列主元消去法，即：将 A|I 通过行变换为 I| $A^{-1}$

/*
get inverse matrix
use main column element of Gauss-Jordan algrithm: A|I  --->  I|A^(-1)
*/
Matrix_t inverse_mat(Matrix_t* m)
{
    if (det_mat(m) == 0)
    {
        printf("error: In inverse_mat(det_mat(mat) == 0)\n");
        exit(1);
    }

    Matrix_t mat = copy_mat(m);
    Matrix_t inv_mat = eye(m->row);

    int i, j, n,max_row;
    int swap_f = 0;
    float max,k;
    for (i = 0; i < mat.row - 1; i++)
    {
        max = abs(mat.data[i][i]);
        max_row = i;
        for (j = i + 1; j < mat.row; j++)
        {
            if (max < abs(mat.data[j][i]))
            {
                max = abs(mat.data[j][i]);
                max_row = j;
            }
        }
        if (i != max_row)
        {
            swap_row_mat(&mat, i, max_row);
            swap_row_mat(&inv_mat, i, max_row);
            swap_f++;
        }
        for (j = i + 1; j < mat.row; j++)
        {
            k = -mat.data[j][i] / mat.data[i][i];
            for (n = 0; n < mat.column; n++)
            {
                mat.data[j][n] = mat.data[i][n] * k + mat.data[j][n];
                inv_mat.data[j][n] = inv_mat.data[i][n] * k + inv_mat.data[j][n];
            }
        }
    }

    for (i = 0; i < mat.row; i++)
    {
        k = 1/mat.data[i][i];
        scale_row_mat(&mat,i, k);
        scale_row_mat(&inv_mat, i, k);
    }
    for (i = mat.row-1; i>0; i--)
    {
        for (j = i - 1; j >=0; j--)
        {
            k = -mat.data[j][i] / mat.data[i][i];
            for (n = 0; n < mat.column; n++)
            {
                mat.data[j][n] = mat.data[j][n] + k*mat.data[i][n];
                inv_mat.data[j][n] = inv_mat.data[j][n] + k*inv_mat.data[i][n];
            }
        }

    }
    //show_mat("in inverse_mat: mat", &mat);
    //show_mat("in inverse_mat: inv_mat", &inv_mat);
    free_mat(&mat);
    return inv_mat;
}

测试一下：

#include<iostream>
#include"myMatrix.h"

using namespace std;

int main()
{
    Matrix_t mat=create_mat(4,4);
    float data[16] = { 1, 2, 3, 4,
                        2, 1, 4, 5,
                        5, 4, 3, 5,
                        7, 6, 5, 4 };
    set_mat_data(&mat, data);
    show_mat("mat", &mat);
    Matrix_t inv_mat=inverse_mat(&mat);
    show_mat("inv_mat", &inv_mat);

    system("pause");
    return 0;
}

这里写图片描述
以上是实现该矩阵运算库的核心代码了，感兴趣的可以下载本代码：
https://download.csdn.net/download/shuoyueqishilove/10433912

基于C语言的矩阵运算库

猜你喜欢