C语言实现复数矩阵基本操作+LU分解求逆

C语言实现（复数）矩阵运算

代码并不是很严谨，没有考虑内存释放之类的问题，初次尝试还有待完善，提供矩阵运算思路、希望能够帮助到有需要的同学。

1.首先创建复数类型结构体

typedef struct
{
    float real;
    float imag;
}Complex_F;

结构体Complex_F中定义了复数的实部与虚部，这里我使用Complex_F *A或者Complex_F A[row*col]来表示一个矩阵，由于是行排列，所以没有使用二维数组，只需知道A[i][j]=A[i*col+j]即可。
下面是所使用的代码，重点是LU分解求逆矩阵，
参考了https://blog.csdn.net/weixin_41645983/

代码实现功能

代码实现了复数的加减乘除，复数矩阵的初始化、拷贝、打印，复数矩阵的乘法、求逆以及转置（复数矩阵一般用共轭转置情形较多，比如要求H矩阵每列范数，可用H共轭转置乘以H，取对角线元素即可）。

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
typedef struct
{
    float real;
    float imag;
}Complex_F;//我们定义一个复数类型，构建Complex_F类型数组

void initMatrix(Complex_F* A, int row, int col);
Complex_F MulComplex(const Complex_F a, const Complex_F b);
Complex_F addComplex(const Complex_F a, const Complex_F b);
void addMatrix(const Complex_F* A, const Complex_F* B, Complex_F* C, int row, int col);
void PrintMatrix(const Complex_F* A, int row, int col);
void MulMatrix(const Complex_F* A, const Complex_F* B, Complex_F* C, int row, int AB_rc, int col);
void TransMatrix(const Complex_F* A, Complex_F* B, int row, int col);
void MatrixCopy(const Complex_F* A, Complex_F* B, int row, int col);
void Inverse_LU(Complex_F* A, int row, int col);


void initMatrix(Complex_F* A, int row, int col)
{
    int i, j;
    for (i = 0; i < row; i++)
    {
        for (j = 0; j < col; j++)
        {
            A[i * col + j].real = 0.0;
            A[i * col + j].imag = 0.0;
        }
    }
    //以上写法是将数组看成矩阵，行排列读取第i行第j列用i*col+j；
}
//当传进去的矩阵元素在函数中不发生改变时，用const Complex_F* A
//只要不想改变实参时，都可以加上const避免出错
Complex_F MulComplex(const Complex_F a, const Complex_F b)
{
    Complex_F c;
    c.real = a.real * b.real - a.imag * b.imag;
    c.imag = a.real * b.imag + a.imag * b.real;
    return c;
}
Complex_F addComplex(const Complex_F a, const Complex_F b)
{
    Complex_F c;
    c.real = a.real + b.real;
    c.imag = a.imag + b.imag;
    return c;
}
Complex_F subComplex(const Complex_F a, const Complex_F b)
{
    Complex_F c;
    c.real = a.real - b.real;
    c.imag = a.imag - b.imag;
    return c;
}
Complex_F divComplex(const Complex_F a, const Complex_F b)
{
    Complex_F c;
    float norm = b.real * b.real + b.imag * b.imag;
    c.real = (a.real * b.real + a.imag * b.imag) / norm;
    c.imag = (a.imag * b.real - a.real * b.imag) / norm;
    return c;
}
void addMatrix(const Complex_F* A, const Complex_F* B, Complex_F* C, int row, int col)
//传递指针可以不用带返回值，我们想得到C的结果，在传入进去前给C分配空间即可，A和B加上const
{
    int i, j;
    for (i = 0; i < row; i++) {
        for (j = 0; j < col; j++)
        {
            C[i * col + j].real = A[i * col + j].real + B[i * col + j].real;
            C[i * col + j].imag = A[i * col + j].imag + B[i * col + j].imag;
        }
    }
}
void PrintMatrix(const Complex_F* A, int row, int col)
{
    int i, j;
    for (i = 0; i < row; i++)
    {
        for (j = 0; j < col; j++)
        {
            printf("%.3f+(%.3f)i\t", A[i * col + j].real, A[i * col + j].imag);
        }
        printf("\n");
    }
}

void MulMatrix(const Complex_F* A, const Complex_F* B, Complex_F* C, int row, int AB_rc, int col)
{
    int i, j, k;
    Complex_F sum, tmp;
    for (i = 0; i < row; i++)
    {
        for (j = 0; j < col; j++)
        {

            sum.real = 0.0;
            sum.imag = 0.0;
            for (k = 0; k < AB_rc; k++)
            {
                sum = addComplex(sum, MulComplex(A[i * AB_rc + k], B[k * col + j]));
            }
            C[i * col + j].real = sum.real;
            C[i * col + j].imag = sum.imag;
        }
    }
}

void MatrixCopy(const Complex_F* A, Complex_F* B, int row, int col)
{
    int i, j;
    for (i = 0; i < row; i++)
    {
        for (j = 0; j < col; j++)
        {
            B[i * col + j] = A[i * col + j];
        }
    }
}
void TransMatrix(const Complex_F* A, Complex_F* B, int row, int col)
{
    int i, j;
    for (i = 0; i < row; i++)
    {
        for (j = 0; j < col; j++)
        {
            B[j * row + i] = A[i * col + j];
        }
    }
}
void TransConjMatrix(const Complex_F* A, Complex_F* B, int row, int col)
{
    int i, j;
    for (i = 0; i < row; i++)
    {
        for (j = 0; j < col; j++)
        {
            B[j * row + i].real = A[i * col + j].real;
            B[j * row + i].imag = -A[i * col + j].imag;
        }
    }
}
//LU分解求逆矩阵
//LU分解将一个矩阵分解成一个单位下三角矩阵和一个上三角矩阵的乘积
//我们希望L的对角线上的元素为1,使用Doolittle分解得到
//对于满秩矩阵A可以通过相乘一个消元矩阵得到一个上三角矩阵U
//A=LU div(A)=div(U)div(L)

void Inverse_LU(Complex_F* A, int row, int col)
{
    Complex_F* L;
    Complex_F* U;
    Complex_F* L_Inverse;
    Complex_F* U_Inverse;
    L = (Complex_F*)malloc(row * col * sizeof(Complex_F));
    U = (Complex_F*)malloc(row * col * sizeof(Complex_F));
    L_Inverse = (Complex_F*)malloc(row * col * sizeof(Complex_F));
    U_Inverse = (Complex_F*)malloc(row * col * sizeof(Complex_F));
    initMatrix(L, row, col);
    initMatrix(U, row, col);
    initMatrix(L_Inverse, row, col);
    initMatrix(U_Inverse, row, col);
    int i, j, k, t;
    Complex_F s;
    for (i = 0; i < row; i++)//对角元素
    {
        L[i * col + i].real = 1.0;
        L[i * col + i].imag = 0.0;
    }
    for (j = 0; j < col; j++)
    {
        U[j] = A[j];
    }
    for (i = 1; i < col; i++)
    {
        L[i * col] = divComplex(A[i * col], U[0]);

    }
    for (k = 1; k < row; k++)
    {
        for (j = k; j < col; j++)
        {
            s.imag = 0.0;
            s.real = 0.0;
            for (t = 0; t < k; t++) {
                s = addComplex(s, MulComplex(L[k * col + t], U[t * col + j]));
            }
            U[k * col + j] = subComplex(A[k * col + j], s);
        }
        for (i = k; i < col; i++)
        {
            s.imag = 0.0;
            s.real = 0.0;
            for (t = 0; t < k; t++)
            {
                s = addComplex(s, MulComplex(L[i * col + t], U[t * col + k]));
            }
            L[i * col + k] = divComplex(subComplex(A[i * col + k], s), U[k * col + k]);
        }
    }

    for (i = 0; i < row; i++)
    {
        L_Inverse[i * col + i].imag = 0.0;
        L_Inverse[i * col + i].real = 1.0;
    }
    for (j = 0; j < col; j++)
    {
        for (i = j + 1; i < row; i++)
        {
            s.imag = 0.0;
            s.real = 0.0;
            for (k = j; k < i; k++) {
                s = addComplex(s, MulComplex(L[i * col + k], L_Inverse[k * col + j]));
            }
            L_Inverse[i * col + j].real = (-1) * MulComplex(L_Inverse[j * col + j], s).real;
            L_Inverse[i * col + j].imag = (-1) * MulComplex(L_Inverse[j * col + j], s).imag;
        }
    }
    Complex_F di;
    di.real = 1.0;
    di.imag = 0.0;
    for (i = 0; i < col; i++)                    //按列序，列内按照从下到上，计算u的逆矩阵
    {
        U_Inverse[i * col + i] = divComplex(di, U[i * col + i]);

    }
    for (j = 0; j < col; j++)
    {
        for (i = j - 1; i >= 0; i--)
        {
            s.imag = 0.0;
            s.real = 0.0;
            for (k = i + 1; k <= j; k++)
            {
                s = addComplex(s, MulComplex(U[i * col + k], U_Inverse[k * col + j]));

            }
            s.imag = -s.imag;
            s.real = -s.real;
            U_Inverse[i * col + j] = divComplex(s, U[i * col + i]);
        }
    }
    MulMatrix(U_Inverse, L_Inverse, A, row, row, col);
}
/*Complex_F Det(const Complex_F* A,int row,int col)
{
    if(row!=col)
    {
        printf("error:please input a square matrix!\n");
        return;
    }
    else{


    }
}*/
int main()
{
    Complex_F A[4] = { {1,0},{10,0},{5,0},{2,0} };
    Complex_F B[4] = { {1.0,2.0},{1.0,2.0},{1.1,1.1},{1.0,2.0} };
    Complex_F* D;
    Complex_F* C;
    Complex_F* A_;
    C = (Complex_F*)malloc(4 * sizeof(Complex_F));//空指针一定要记得分配空间！！
    D = (Complex_F*)malloc(4 * sizeof(Complex_F));
    A_ = (Complex_F*)malloc(4 * sizeof(Complex_F));
    //MulMatrix(A,B,C,2,2,2);
    //PrintMatrix(C,2,2);
    //TransConjMatrix(A,D,2,2);//计算A矩阵每列的范数平方可以用A的转置乘以A，取对角元素得到，复数矩阵则是共轭转置
    //MulMatrix(D,A,C,2,2,2);
    //PrintMatrix(C,2,2);
    Inverse_LU(A, 2, 2);
    PrintMatrix(A, 2, 2);

    return 0;
}