向量与矩阵计算的C语言代码

最新推荐文章于 2024-08-04 17:17:39 发布

小滕摸摸头

最新推荐文章于 2024-08-04 17:17:39 发布

阅读量3.3k

点赞数 4

分类专栏： c 文章标签：矩阵 c语言线性代数

本文链接：https://blog.csdn.net/weixin_42248256/article/details/120886121

版权

c 专栏收录该内容

2 篇文章 1 订阅

订阅专栏

1.初始化定义

typedef double ElemType;
typedef ElemType* Vector;
typedef ElemType** Matrix;
typedef ElemType*** MatrixStack;

#define ZERO_MATRIX 0
#define ONES_MATRIX 1
#define UNIT_MATRIX 2

2.向量操作

2.1向量初始化

/* brief：创建向量
** param：长度
** output：Vector
*/
Vector initVector(int len) {
    Vector res = (ElemType*)malloc(sizeof(ElemType)*len);
    if (res == NULL) {
        printf("VECTOR CREATE ERROR\n");
        return NULL;
    }
    return res;
}

2.2向量销毁

/* brief：销毁向量
** param：向量
** output：整数
*/
int destoryVector(Vector vec) {
    free(vec);
    return 1;
}

2.3向量相加

/* 说明：两个向量的加法
** param[0]：向量A
** param[1]：向量B
** param[2]：长度
** output：向量
*/
Vector vectorAdd(const Vector vecA, const Vector vecB, int len) {
    Vector C = initVector(len);
    for (int i = 0; i<len; i++) {
        C[i] = vecA[i] + vecB[i];
    }
    return C;
}

2.4向量相减

/* 说明：两个向量的减法
** param[0]：向量A
** param[1]：向量B
** param[2]：长度
** output：向量
*/
Vector vectorSubtract(const Vector vecA, const Vector vecB, int len) {
    Vector C = initVector(len);
    for (int i = 0; i<len; i++) {
        C[i] = vecA[i] - vecB[i];
    }
    return C;
}

2.5向量相乘(1,N)*(N,1)

/* brief：向量相乘(1,N)*(N,1)
** param[0]：向量A
** param[1]：向量B
** param[2]：向量长度
** output：浮点数
*/
ElemType vecDotVecToElem(const Vector vecA, const Vector vecB, int n) {
    int i, i1;
    ElemType ret = vecA[n - 1] * vecB[n - 1];
    for (i = n - 2; i >= 1; i -= 2) {
        i1 = i - 1;
        ret += vecA[i] * vecB[i] + vecA[i1] * vecB[i1];
    }
    if (i == 0) ret += vecA[0] * vecB[0]; 
    return ret;
}

2.6向量相乘(N,1)*(1,N)

/* brief：向量相乘(N,1)*(1,N)
** param[0]：向量A
** param[1]：向量A长度
** param[2]：向量B
** param[3]：向量B长度
** output：矩阵
*/
Matrix vecDotVecToMatrix(const Vector vecA, int aLen, const Vector vecB, int bLen) {
    Matrix res=(ElemType**)malloc(sizeof(ElemType*)*aLen);
    for(int i=0;i<aLen;i++){
        res[i]=(ElemType*)malloc(sizeof(ElemType)*bLen);
        for(int j=0;j<bLen;j++){
            res[i][j]=vecA[i]*vecB[j];
        }
    }
    return res;
}

2.7向量乘数值

/* brief：Vector乘数值
** param[0]：向量
** param[1]：长度
** param[2]：数值
** output：向量
*/
Vector vectorMulNum(const Vector vec, int len, ElemType num) {
    Vector res = initVector(len);
    for (int i = 0; i<len; i++) {
        res[i] = vec[i] * num;
    }
    return res;
}

2.8向量除以数值

/* brief：Vector除以数值
** param[0]：向量
** param[1]：长度
** param[2]：数值
** output：向量
*/
Vector vectorDivNum(const Vector vec, int len, ElemType num) {
    Vector res = initVector(len);
    for (int i = 0; i<len; i++) {
        res[i] = vec[i] / num;
    }
    return res;
}

2.9向量的均方根

/* brief：向量的均方根
** param[0]：向量
** param[1]：长度
** output：均方根
*/
ElemType vectorRootMeanSquare(const Vector vec, int len) {
    ElemType nor;
    for (int i = 0; i<len; i++) {
        nor += vec[i] * vec[i];
    }
    return sqrt(nor);
}

2.10打印向量

/* brief：打印向量
** param[0]：向量
** param[1]：长度
** output：无
*/
void vectorPrint(const Vector vec, int len) {
    for (int i = 0; i<len; i++) {
        printf("%lf\t", vec[i]);
    }
    printf("\n");
}

3.矩阵操作

3.1创建矩阵

/* brief：创建矩阵
** param[0]：行
** param[1]：列
** output：矩阵
*/
Matrix initMatrix(int row, int col) {
    Matrix mat = (ElemType**)malloc(sizeof(ElemType*)*row);
    if (mat == NULL) {
        printf("MATRIX CREATE ERROR\n");
        return NULL;
    }
    for (int i = 0; i<row; i++) {
        mat[i] = (ElemType*)malloc(sizeof(ElemType)*col);
        if (mat[i] == NULL) {
            printf("MATRIX CREATE ERROR\n");
            return NULL;
        }
    }
    return mat;
}

3.2销毁矩阵

/* brief：销毁矩阵
** param[0]：数组
** param[1]：行
** output：整数
*/
int destroyMatrix(Matrix mat, int row) {
    if (!mat) return 1;
    for (int i = 0; i<row; i++) {
        free(mat[i]);
    }
    free(mat);
    mat = NULL;
    return 1;
}

3.3填充矩阵

/* brief：填充矩阵
** param[0]：矩阵
** param[1]：行
** param[2]：列
** param[3]：类型[ZERO_MATRIX, ONES_MATRIX, UNIT_MATRIX]
** 输出：矩阵
*/
Matrix fillMatrix(Matrix mat, int row, int col, int type) {
    int i, j, k;
    switch (type) {
    case ZERO_MATRIX:
        for (i = 0; i<row; ++i) {
            for (j = 0; j<col; ++j) {
                mat[i][j] = 0.0;
            }
        }
        break;
    case UNIT_MATRIX:
        for (i = 0; i<row; ++i) {
            for (j = 0; j<col; ++j) {
                mat[i][j] = 0.0;
            }
        }
        k = (col<row) ? col : row;
        for (i = 0; i<k; ++i) {
            mat[i][i] = 1.0;
        }
        break;
    case ONES_MATRIX:
        for (i = 0; i<row; ++i) {
            for (j = 0; j<col; ++j)
                mat[i][j] = 1.0;
        }
        break;
    }
    return mat;
}

3.4矩阵获取行

/* brief：浮点矩阵获取行
** param[0]：矩阵
** param[1]：行
** param[2]：列
** param[3]：行定位
** output：Vector
*/
Vector getRow(const Matrix mat, int row, int col, int loc) {
    if (loc >= row) {
        printf("GET ROW ERROR\n");
        return NULL;
    }
    int i;
    Vector x = initVector(col);
    for (i = col - 1; i>0; --i) {
        x[i] = mat[loc][i];
        --i;
        x[i] = mat[loc][i];
    }
    if (i == 0) x[0] = mat[loc][0];
    return x;
}

3.5矩阵获取列

/* brief：浮点矩阵获取列
** param[0]：矩阵
** param[1]：行
** param[2]：列
** param[3]：列定位
** output：Vector
*/
Vector getCol(const Matrix mat, int row, int col, int loc) {
    if (loc >= col) {
        printf("GET COLUMN ERROR\n");
        return NULL;
    }
    int i;
    Vector x = initVector(row);
    for (i = row - 1; i>0; --i) {
        x[i] = mat[i][loc];
        --i;
        x[i] = mat[i][loc];
    }
    if (i == 0) x[0] = mat[0][loc];
    return x;
}

3.6设置行

/* brief：设置行
** param[0]：矩阵
** param[1]：行
** param[2]：列
** param[3]：插入的行
** param[4]：插入的位置
** output：1
*/
int setRow(Matrix A, int row, int col, Vector V, int loc){
    if(loc<0 || loc>row){
        return 0;
    }
    for(int i=0;i<col;i++){
        A[loc][i]=V[i];
    }
    return 1;
}

3.7设置列

/* brief：设置列
** param[0]：矩阵
** param[1]：行
** param[2]：列
** param[3]：插入的列
** param[4]：插入的位置
** output：1
*/
int setCol(Matrix A, int row, int col, Vector V, int loc){
    if(loc<0 || loc>col){
        return 0;
    }
    for(int i=0;i<row;i++){
        A[i][loc]=V[i];
    }
    return 1;
}

3.8矩阵加法

/* brief：矩阵加法
** param[0]：矩阵A
** param[1]：矩阵B
** param[2]：矩阵A行
** param[3]：矩阵A列
** param[4]：矩阵B行
** param[5]：矩阵B列
** output：矩阵
*/
Matrix matrixAdd(const Matrix x, const Matrix y, int xrow, int xcol, int yrow, int ycol) {
    if (xrow != yrow || xcol != ycol) {
        printf("MATRIX ADD ERROR\n");
        return NULL;
    }
    Matrix res = (ElemType**)malloc(sizeof(ElemType*)*xrow);
    for (int i = 0; i<xrow; i++) {
        res[i] = (ElemType*)malloc(sizeof(ElemType)*xcol);
        for (int j = 0; j<xcol; j++) {
            res[i][j] = x[i][j] + y[i][j];
        }
    }
    return res;
}

3.9矩阵减法

/* brief：矩阵减法
** param[0]：矩阵A
** param[1]：矩阵B
** param[2]：矩阵A行
** param[3]：矩阵A列
** param[4]：矩阵B行
** param[5]：矩阵B列
** output：矩阵
*/
Matrix matrixSubtract(const Matrix x, const Matrix y, int xrow, int xcol, int yrow, int ycol) {
    if (xrow != yrow || xcol != ycol) {
        printf("MATRIX SUBTRACT ERROR\n");
        return NULL;
    }
    Matrix res = (ElemType**)malloc(sizeof(ElemType*)*xrow);
    for (int i = 0; i<xrow; i++) {
        res[i] = (ElemType*)malloc(sizeof(ElemType)*xcol);
        for (int j = 0; j<xcol; j++) {
            res[i][j] = x[i][j] - y[i][j];
        }
    }
    return res;
}

3.10矩阵乘以数值

/* brief：矩阵乘以数值
** param[0]：矩阵
** param[1]：行
** param[2]：列
** param[3]：数值
** output：数组
*/
Matrix matrixMulNum(const Matrix mat, int row, int col, ElemType num) {
    Matrix res = (ElemType**)malloc(sizeof(ElemType*)*row);
    for (int i = 0; i<row; i++) {
        res[i] = (ElemType*)malloc(sizeof(ElemType)*col);
        for (int j = 0; j<col; j++) {
            res[i][j] = mat[i][j] * num;
        }
    }
    return res;
}

3.11矩阵乘法

Matrix _marixDotMatrixSmall(const Matrix x, const Matrix y, int xrow, int xcol, int ycol) {
    Matrix ret = (ElemType**)malloc(sizeof(ElemType*)*xrow);
    for (int i = 0; i<xrow; i++) {
        ret[i] = (ElemType*)malloc(sizeof(ElemType)*ycol);
    }
    for (int i = xrow - 1; i >= 0; i--) {
        for (int k = ycol - 1; k >= 0; k--) {
            int j;
            ElemType woo = x[i][xcol - 1] * y[xcol - 1][k];
            for (j = xcol - 2; j >= 1; j -= 2) {
                int i0 = j - 1;
                woo += x[i][j] * y[j][k] + x[i][i0] * y[i0][k];
            }
            if (j == 0) woo += x[i][0] * y[0][k];
            ret[i][k] = woo;
        }
    }
    return ret;
}

Matrix _marixDotMatrixBig(const Matrix x, const Matrix y, int xrow, int xcol, int ycol) {
    int m = xrow, n = ycol;
    Matrix A = initMatrix(m, n);
    for (int i = n - 1; i != -1; --i) {
        Vector v = getCol(y, xrow, xcol, i);
        for (int j = m - 1; j != -1; --j) {
            Vector xj = getRow(x, xrow, xcol, j);
            A[j][i] = vecDotVecToElem(xj, v, n);
            destoryVector(xj);
        }
        destoryVector(v);
    }
    return A;
}

/* brief：矩阵乘法
** param[0]：矩阵A
** param[1]：矩阵B
** param[2]：矩阵A行
** param[3]：矩阵A列
** param[4]：矩阵B行
** param[5]：矩阵B列
** 输出：矩阵相乘结果
*/
Matrix matrixMultiply(const Matrix x, int xrow, int xcol, const Matrix y, int yrow, int ycol) {
    if (xcol != yrow) {
        printf("MATRIX MULTIPLY ERROR\n");
        return NULL;
    }
    if (xcol < 10) {
        return _marixDotMatrixSmall(x, y, xrow, xcol, ycol);
    }
    else {
        return _marixDotMatrixBig(x, y, xrow, xcol, ycol);
    }
}

3.12矩阵转置

/* brief：矩阵转置
** param[0]：矩阵
** param[1]：行
** param[2]：列
** output：转置结果
*/
Matrix matrixTranspose(const Matrix mat, int row, int col) {
    Matrix res = (ElemType**)malloc(sizeof(ElemType*)*col);
    for (int i = 0; i<col; i++) {
        res[i] = (ElemType*)malloc(sizeof(ElemType)*row);
        for (int j = 0; j<row; j++) {
            res[i][j] = mat[j][i];
        }
    }
    return res;
}

3.13矩阵复制

/* brief：矩阵复制
** param[0]：待复制矩阵
** param[1]：行
** param[2]：列
** output：复制结果
*/
Matrix matrixCopy(const Matrix mat, int row, int col) {
    Matrix result = initMatrix(row, col);
    for (int i = 0; i < row; i++) {
        for (int j = 0; j < col; j++) {
            result[i][j] = mat[i][j];
        }
    }
    return result;
}

3.14求行列式

/* brief：求行列式
** param[1]：输入方阵
** param[2]：方阵的行数
** 输出：行列式结果
*/
ElemType matrixDet(const Matrix mat, int n) {
    if (n == 1) {
        return mat[0][0];
    }
    ElemType ans = 0.0;
    Matrix temp = initMatrix(n, n);
    int i, j, k;
    for (i = 0; i<n; i++) {
        for (j = 0; j<n - 1; j++) {
            for (k = 0; k<n - 1; k++) {
                temp[j][k] = mat[j + 1][(k >= i) ? k + 1 : k];

            }
        }
        ElemType t = matrixDet(temp, n - 1);
        if (i % 2 == 0) {
            ans += mat[0][i] * t;
        }
        else {
            ans -= mat[0][i] * t;
        }
    }
    return ans;
}

3.15矩阵取逆

//计算每一行每一列的每个元素所对应的余子式，组成A*
void _getMatrixStart(Matrix arcs, int row, int col, Matrix ans) {
    if (row != col) printf("GET A* ERROR\n");
    int n = row;
    if (n == 1) {
        ans[0][0] = 1;
        return;
    }
    int i, j, k, t;
    Matrix temp = initMatrix(row, col);
    for (i = 0; i<n; i++) {
        for (j = 0; j<n; j++) {
            for (k = 0; k<n - 1; k++) {
                for (t = 0; t<n - 1; t++) {
                    temp[k][t] = arcs[k >= i ? k + 1 : k][t >= j ? t + 1 : t];
                }
            }
            ans[j][i] = matrixDet(temp, n - 1);
            if ((i + j) % 2 == 1) {
                ans[j][i] = -ans[j][i];
            }
        }
    }
}

Matrix _matrixInverse(Matrix mat, int row, int col) {
    ElemType det = matrixDet(mat, row);
    Matrix result = initMatrix(row, col);
    Matrix temp = initMatrix(row, col);
    if (det == 0) {
        printf("GET INVERSE ERROR\n");
        return NULL;
    }
    else {
        _getMatrixStart(mat, row, col, temp);
        for (int i = 0; i<row; i++) {
            for (int j = 0; j<col; j++) {
                result[i][j] = temp[i][j] / det;
            }
        }
    }
    return result;
}

/* brief：矩阵取逆
** param[0]：待取逆矩阵
** param[1]：行
** param[2]：列
** output：取逆结果
*/
Matrix matrixInverse(const Matrix mat, int row, int col) {
    Matrix tmp1, tmp2, tmp3, tmp4;
    Matrix B, X, X_H, B_R;
    Matrix ret;

    if (row < col) {
        tmp1 = matrixTranspose(mat, row, col);
        B = matrixMultiply(mat, row, col, tmp1, row, col);

        tmp2 = matrixMultiply(B, row, col, B, row, col);
        tmp3 = _matrixInverse(tmp2, row, col);
        X_H = matrixMultiply(tmp3, row, col, B, row, col);

        X = matrixTranspose(X_H, row, col);
        tmp4 = matrixMultiply(X, row, col, B, row, col);
        B_R = matrixMultiply(tmp4, row, col, X_H, row, col);

        ret = matrixMultiply(tmp1, row, col, B_R, row, col);

        destroyMatrix(tmp1, row);
        destroyMatrix(tmp2, row);
        destroyMatrix(tmp3, row);
        destroyMatrix(tmp4, row);
        destroyMatrix(X_H, row);
        destroyMatrix(X, row);
        destroyMatrix(B, row);
        destroyMatrix(B_R, row);
    }
    else if (row > col) {
        tmp1 = matrixTranspose(mat, row, col);
        B = matrixMultiply(tmp1, row, col, mat, row, col);

        tmp2 = matrixMultiply(B, row, col, B, row, col);
        tmp3 = _matrixInverse(tmp2, row, col);
        X_H = matrixMultiply(tmp3, row, col, B, row, col);

        X = matrixTranspose(X_H, row, col);
        tmp4 = matrixMultiply(X, row, col, B, row, col);
        B_R = matrixMultiply(tmp4, row, col, X_H, row, col);

        ret = matrixMultiply(B_R, row, col, tmp1, row, col);

        destroyMatrix(tmp1, row);
        destroyMatrix(tmp2, row);
        destroyMatrix(tmp3, row);
        destroyMatrix(tmp4, row);
        destroyMatrix(X_H, row);
        destroyMatrix(X, row);
        destroyMatrix(B, row);
        destroyMatrix(B_R, row);
    }
    else {
        ret = _matrixInverse(mat, row, col);
    }
    return ret;
}

3.16求每一列的均值

/* brief：求每一列的均值
** param[0]：输入矩阵
** param[1]：行
** param[2]：列
** output：均值矩阵
*/
Matrix matrixMeanCol(const Matrix mat, int row, int col) {
    int i, j;
    Matrix result = initMatrix(1, col);
    // 第一行设置为0
    for (j = 0; j < col; j++)
        result[0][j] = 0.0;
    for (i = 0; i < row; i++) {
        for (j = 0; j < col; j++)
            result[0][j] += mat[i][j];
    }
    for (j = 0; j < col; j++)
        result[0][j] /= (ElemType)row;
    return result;
}

3.17矩阵标准化？

/* brief：矩阵标准化？
** param[0]: 输入矩阵
** param[1]: 均值矩阵
** param[2]：输入矩阵行
** param[3]：输入矩阵列
** param[4]：均值矩阵行
** param[5]：均值矩阵列
** output：标准化矩阵
*/
Matrix matrixStd(const Matrix mat, const Matrix mean, int rowMat, int colMat, int rowMean, int colMean) {
    int i, j, row, col;
    Matrix std = initMatrix(rowMean, colMean);
    for (j = 0; j < colMat; j++) {
        for (i = 0; i < rowMat; i++) {
            std[0][j] += (mat[i][j] - mean[0][j]) * (mat[i][j] - mean[0][j]);
        }
    }
    for (i = 0; i < colMean; i++) {
        std[0][i] = std[0][i] / (((ElemType)rowMean) - (ElemType)1.0);
        std[0][i] = sqrt(std[0][i]);
    }
    return std;
}

3.18打印矩阵

/* brief：打印矩阵
** param[0]: 输入矩阵
** param[1]：行
** param[2]：列
** output：无
*/
void matrixPrint(const Matrix mat, int row, int col) {
    for (int i = 0; i<row; i++) {
        for (int j = 0; j<col; j++) {
            printf("%f\t", mat[i][j]);
        }
        printf("\n");
    }
    printf("\n");
}

3.19对角化

/* brief：对角化
** param[0]: 输入方阵
** param[1]：行
** output：无
*/
Matrix diag(Vector d, int N){
    int i1,j;
    Matrix A = initMatrix(N,N);
    for(int i=N-1;i>=0;i--) {
        i1 = i+2;
        for(j=N-1;j>=i1;j-=2) {
            A[i][j] = 0;
            A[i][j-1] = 0;
        }
        if(j>i) { A[i][j] = 0; }
        A[i][i] = d[i];
        for(j=i-1;j>=1;j-=2) {
            A[i][j] = 0;
            A[i][j-1] = 0;
        }
        if(j==0) { A[i][0] = 0; }
    }
    return A;
}

小滕摸摸头

关注

4
点赞
踩
45

收藏

觉得还不错? 一键收藏
2
评论
向量与矩阵计算的C语言代码

1.初始化定义typedef double ElemType;typedef ElemType* Vector;typedef ElemType** Matrix;typedef ElemType*** MatrixStack;#define ZERO_MATRIX 0#define ONES_MATRIX 1#define UNIT_MATRIX 22.向量操作2.1向量初始化/* brief：创建向量** param：长度** output：Vector*/Vector
复制链接

扫一扫

专栏目录