C语言矩阵运算

这是一个用C语言编写的矩阵运算头文件，话不多说，上代码

matrix.h

/**
 * File name: matrix.h
 * Author: Tiancai Xu
 * Version: 1.0
 * Data: 2021-10-29
 * Copyright: 2021.10, Tiancai Xu
 * This file defines the matrix and its basic operations
 * This file is for study or reference only
*/

#ifndef _STDARG_H
#include <stdarg.h>
#endif // _STDARG_H

#ifndef _STDIO_H
#include <stdio.h>
#endif // _STDIO_H

#ifndef _STDLIB_H
#include <stdlib.h>
#endif // _STDLIB_H

#ifndef _TIME_H
#include <time.h>
#endif // _TIME_H

#ifndef _MATRIX_H
#define _MATRIX_H

#define size_m unsigned int

/**< Matrix variable */
typedef struct
{
    double** M; /**< Data of Matrix */
    size_m row; /**< The number of row */
    size_m col; /**< The number of column */
} Matrix;
/**< Tips: This format is recommended for declarations: Matrix M = {NULL, 0, 0}; */

Matrix from(double* arr, size_m row, size_m col);
Matrix zeros(size_m row, size_m col);
Matrix ones(size_m row, size_m col);
Matrix diag(int order, ...);
Matrix diag_M(Matrix A);
Matrix iden(int order);
Matrix randM(size_m row, size_m col);
Matrix add(Matrix A, Matrix B);
Matrix sub(Matrix A, Matrix B);
Matrix mult_MM(Matrix A, Matrix B);
Matrix mult_kM(double k, Matrix B);
Matrix T(Matrix A);
Matrix inv(Matrix A);
Matrix cof(Matrix A, int r, int c);
double mode(Matrix A);
void scanM(Matrix* A);
void printM(Matrix A);
void clearM(Matrix* A);

#endif // _MATRIX_H

matrix.c

/**
 * File name: matrix.c
 * Author: Tiancai Xu
 * Version: 1.0
 * Data: 2021-10-29
 * Copyright: 2021.10, Tiancai Xu
 * This file includes the implementation of the functions declarated in matrix.h
 * This file is for study or reference only
*/

#include "matrix.h"

/** \brief Create a block of 2-D array memory
 *
 * \param row The number of row
 * \param col The number of column
 * \return M 2-D array pointer
 *
 */
static double** getM(size_m row, size_m col)
{
    int i;
    double** M;

    M = (double**)malloc(sizeof(double*) * row);

    for (i = 0; i < row; i++)
        *(M + i) = (double*)malloc(sizeof(double) * col);

    return M;
}

/** \brief Transform a 2-D array to Matrix
 *
 * \param arr The first adress of the a-D array
 * \param row The number of row
 * \param col The number of column
 * \return M Matrix
 *
 * For a 2-D array arr[r][c], use from((double*)arr,r,c) to avoid the gcc compiler warning
 */
Matrix from(double* arr, size_m row, size_m col)
{
    int i, j;
    Matrix M;
    M.row = row;
    M.col = col;
    M.M = getM(row, col);

    for (i = 0; i < row; i++)
        for (j = 0; j < col; j++)
            M.M[i][j] = *(arr + i * col + j);

    return M;
}

/** \brief Create a zeros Matrix
 *
 * \param row The number of row
 * \param col The number of column
 * \return M Matrix
 *
 */
Matrix zeros(size_m row, size_m col)
{
    int i, j;
    Matrix M;
    M.row = row;
    M.col = col;
    M.M = getM(row, col);

    for (i = 0; i < row; i++)
        for (j = 0; j < col; j++)
            M.M[i][j] = 0;

    return M;
}

/** \brief Create a ones Matrix
 *
 * \param row The number of row
 * \param col The number of column
 * \return M Matrix
 *
 */
Matrix ones(size_m row, size_m col)
{
    int i, j;
    Matrix M;
    M.row = row;
    M.col = col;
    M.M = getM(row, col);

    for (i = 0; i < row; i++)
        for (j = 0; j < col; j++)
            M.M[i][j] = 1;

    return M;
}

/** \brief Create a diagonal Matrix
 *
 * \param order The order of the Matrix
 * \return M Matrix
 *
 */
Matrix diag(int order, ...)
{
    int i;
    va_list valist;
    Matrix M;
    M = zeros(order, order);

    va_start(valist, order);

    for (i = 0; i < order; i++)
        M.M[i][i] = va_arg(valist, double);

    return M;
}

/** \brief Create a diagonal Matrix
 *
 * \param A It's a row or column vector
 * \return M Matrix
 *
 */
Matrix diag_M(Matrix A)
{
    int i, order;
    Matrix M;
    order = A.row * A.col;
    M = zeros(order, order);

    if (A.row == 1)
        for (i = 0; i < order; i++)
            M.M[i][i] = A.M[0][i];
    else if (A.col == 1)
        for (i = 0; i < order; i++)
            M.M[i][i] = A.M[i][0];
    else
    {
        fprintf(stderr, "Error: Arguments must be a row or column vector.\n");
        exit(-1);
    }

    return M;
}

/** \brief Create a identity Matrix
 *
 * \param order The order of the Matrix
 * \return M Matrix
 *
 */
Matrix iden(int order)
{
    int i;
    Matrix M;
    M = zeros(order, order);

    for (i = 0; i < order; i++)
        M.M[i][i] = 1;

    return M;
}

/** \brief Create a random Matrix between 0 and 1.0
 *
 * \param row The number of row
 * \param col The number of column
 * \return M Matrix
 *
 */
Matrix randM(size_m row, size_m col)
{
    time_t t;
    int i, j;
    Matrix M = zeros(row, col);
    srand((unsigned) time(&t));

    for (i = 0; i < row; i++)
        for (j = 0; j < col; j++)
            M.M[i][j] = ((double) rand()) / ((double)RAND_MAX);
    return M;
}

/** \brief Matrix addition
 *
 * \param A Matrix
 * \param B Matrix
 * \return M Matrix, M = A + B
 *
 */
Matrix add(Matrix A, Matrix B)
{
    if (!(A.row == B.row && A.col == B.col))
    {
        fprintf(stderr, "Error: The two matrices must have the same size.\n");
        exit(-1);
    }

    int i, j;
    Matrix M;
    M.row = A.row;
    M.col = A.col;
    M.M = getM(M.row, M.col);

    for (i = 0; i < M.row; i++)
        for (j = 0; j < M.col; j++)
            M.M[i][j] = A.M[i][j] + B.M[i][j];

    return M;
}

/** \brief Matrix subtraction
 *
 * \param A Matrix
 * \param B Matrix
 * \return M Matrix, M = A - B
 *
 */
Matrix sub(Matrix A, Matrix B)
{
    if (!(A.row == B.row && A.col == B.col))
    {
        fprintf(stderr, "Error: The two matrices must have the same size.\n");
        exit(-1);
    }

    int i, j;
    Matrix M;
    M.row = A.row;
    M.col = A.col;
    M.M = getM(M.row, M.col);

    for (i = 0; i < M.row; i++)
        for (j = 0; j < M.col; j++)
            M.M[i][j] = A.M[i][j] - B.M[i][j];

    return M;
}

/** \brief Matrix multiplication
 *
 * \param A Matrix
 * \param B Matrix
 * \return M Matrix, M = A * B
 *
 */
Matrix mult_MM(Matrix A, Matrix B)
{
    if (A.col != B.row)
    {
        fprintf(stderr, "Error: The number of columns in matrix A must be equal to the number of rows in matrix B.\n");
        exit(-1);
    }

    int i, j, k;
    Matrix M;
    M.row = A.row;
    M.col = B.col;
    M.M = getM(M.row, M.col);

    for (i = 0; i < M.row; i++)
        for (j = 0; j < M.col; j++)
        {
            M.M[i][j] = 0;

            for (k = 0; k < A.col; k++)
                M.M[i][j] += A.M[i][k] * B.M[k][j];
        }

    return M;
}

/** \brief Matrix multiplication
 *
 * \param k Real number
 * \param A Matrix
 * \return M Matrix, M = k * A
 *
 */
Matrix mult_kM(double k, Matrix B)
{
    int i, j;
    Matrix M;
    M.row = B.row;
    M.col = B.col;
    M.M = getM(M.row, M.col);

    for (i = 0; i < M.row; i++)
        for (j = 0; j < M.col; j++)
            M.M[i][j] = k * B.M[i][j];

    return M;
}

/** \brief Matrix transpose
 *
 * \param A Matrix
 * \return M Matrix, M = A^T
 *
 */
Matrix T(Matrix A)
{
    int i, j;
    Matrix M;
    M.row = A.col;
    M.col = A.row;
    M.M = getM(M.row, M.col);

    for (i = 0; i < M.row; i++)
        for (j = 0; j < M.col; j++)
            M.M[i][j] = A.M[j][i];

    return M;
}

/** \brief Matrix inversion
 *
 * \param A Matrix
 * \return M Inverse Matrix, M = A^-1
 *
 */
Matrix inv(Matrix A)
{
    if (A.row != A.col)
    {
        fprintf(stderr, "Error: The matrix must be square.\n");
        exit(-1);
    }

    int m = mode(A);

    if (m == 0)
    {
        fprintf(stderr, "Error: The inverse of the target matrix does not exist.\n");
        exit(-1);
    }

    int i, j, s;
    Matrix M, temp;
    M.row = A.row;
    M.col = A.col;
    M.M = getM(M.row, M.col);

    for (i = 0; i < M.row; i++)
        for (j = 0; j < M.col; j++)
        {
            temp = cof(A, i, j);
            s = (i + j) % 2 == 0 ? 1 : -1;
            M.M[j][i] = s * mode(temp) / m;
            clearM(&temp);
        }

    return M;
}

/** \brief Complementary minor, cofactor
 *
 * \param A Matrix
 * \param r The index of row
 * \param c The index of column
 * \return M The cofactor Matrix M[r][c]
 *
 */

Matrix cof(Matrix A, int r, int c)
{
    if (A.row != A.col)
    {
        fprintf(stderr, "Error: The matrix must be square.\n");
        exit(-1);
    }

    if (r < 0 || r > A.row)
    {
        fprintf(stderr, "Error: Out of matrix index range.(r)\n");
        exit(-1);
    }

    if (c < 0 || c > A.col)
    {
        fprintf(stderr, "Error: Out of matrix index range.(c)\n");
        exit(-1);
    }

    int i, j, m, n;
    Matrix M;
    M.row = A.row - 1;
    M.col = A.col - 1;
    M.M = getM(M.row, M.col);

    for (i = 0, m = 0; i < A.row; i++)
    {
        if (i == r)
            continue;

        for (j = 0, n = 0; j < A.col; j++)
        {
            if (j == c)
                continue;

            M.M[m][n] = A.M[i][j];
            n++;
        }

        m++;
    }

    return M;
}

/** \brief Get the mode of Matrix
 *
 * \param A Matrix
 * \return result The mode, result = |A|
 *
 */
double mode(Matrix A)
{
    if (A.row != A.col)
    {
        fprintf(stderr, "Error: The matrix must be square.\n");
        exit(-1);
    }

    if (A.row == 1)
        return A.M[0][0];

    double result;
    int i, s;
    Matrix M;

    for (result = 0, i = 0; i < A.row; i++)
    {
        s = i % 2 == 0 ? 1 : -1;
        M = cof(A, i, 0);
        result += A.M[i][0] * s * mode(M);
        clearM(&M);
    }

    return result;
}

/** \brief Input a Matrix
 *
 * \param A Matrix pointer
 * \return void
 *
 */
void scanM(Matrix* A)
{
    int i, j;

    if (A->M == NULL)
    {
        printf("Enter the number of row&column, e.g. 5,5: ");
        scanf("%u,%u", &A->row, &A->col);
    }

    for (i = 0; i < A->row; i++)
        for (j = 0; j < A->col; j++)
        {
            printf("Enter M(%d,%d): ", i, j);
            scanf("%lf", &A->M[i][j]);
        }
}

/** \brief Print a Matrix
 *
 * \param A Matrix
 * \return void
 *
 */
void printM(Matrix A)
{
    int i, j;

    for (i = 0; i < A.row; i++)
    {
        for (j = 0; j < A.col; j++)
            printf("%-10.4g ", A.M[i][j]);

        printf("\n");
    }
}

/** \brief Clear a Matrix
 *
 * \param A Matrix pointer
 * \return void
 *
 * To release occupied memory
 */
void clearM(Matrix* A)
{
    int i;

    if (A->M != NULL)
    {
        for (i = 0; i < A->row; i++)
            free(A->M[i]);

        free(A->M);
    }

    A->M = NULL;
    A->row = 0;
    A->col = 0;
}

测试用test.c

/**< IDE: Code::Blocks 20.03 64bit */
#include <stdio.h>
#include "matrix.h"

int main()
{
    Matrix M = {NULL, 0, 0};
    Matrix M1 = {NULL, 0, 0};
    Matrix M2 = {NULL, 0, 0};

    printf("/*****************************/\n");

    printf("\nM = zeros(2, 3)\n");
    M = zeros(2, 3);
    printM(M);
    clearM(&M);

    printf("\n/*****************************/\n");

    printf("\nM = ones(3, 2)\n");
    M = ones(3, 2);
    printM(M);
    clearM(&M);

    printf("\n/*****************************/\n");

    printf("\nM = iden(3)\n");
    M = iden(3);
    printM(M);
    clearM(&M);

    printf("\n/*****************************/\n");

    printf("\nM = diag(3, 1.0, 2.0, 3.1)\n");
    M = diag(3, 1.0, 2.0, 3.1);
    printM(M);
    clearM(&M);

    printf("\n/*****************************/\n");

    printf("\nM = randM(2, 2)\n");
    M = randM(2, 2);
    printM(M);
    clearM(&M);

    printf("\n/*****************************/\n");

    printf("\nM = from((double*)arr, 2, 3)\n");
    double arr[3][3] = {{1.0, 2.0, 3.0}, {2.0, 2.0, 1.0}, {3.0, 4.0, 3.0}};
    M = from((double*)arr, 3, 3);
    printM(M);
    M1 = T(M);
    printf("\nM1 = T(M)\n");
    printM(M1);
    printf("\nM2 = mult_MM(M1, M)\n");
    M2 = mult_MM(M1, M);
    printM(M2);
    clearM(&M1);
    clearM(&M2);
    printf("\nM1 = inv(M)\n");
    M1 = inv(M);
    printM(M1);
    printf("\nM1 = mult_MM(M, M1)\n");
    M2 = mult_MM(M, M1);
    printM(M2);
    clearM(&M);
    clearM(&M1);
    clearM(&M2);

    printf("\n/*****************************/\n");

    return 0;
}