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

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

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

Vector initVector(int len) {
    Vector res = (ElemType*)malloc(sizeof(ElemType)*len);
    if (res == NULL) {
        printf("VECTOR CREATE ERRORn");
        return NULL;
    }
    return res;
}

int destoryVector(Vector vec) {
    free(vec);
    return 1;
}

Vector vectorAdd(const Vector vecA, const Vector vecB, int len) {
    Vector C = initVector(len);
    for (int i = 0; i 
2.4向量相减 
Vector vectorSubtract(const Vector vecA, const Vector vecB, int len) {
    Vector C = initVector(len);
    for (int i = 0; i 
2.5向量相乘(1,N)*(N,1) 
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) 
Matrix vecDotVecToMatrix(const Vector vecA, int aLen, const Vector vecB, int bLen) {
    Matrix res=(ElemType**)malloc(sizeof(ElemType*)*aLen);
    for(int i=0;i 
2.7向量乘数值 
Vector vectorMulNum(const Vector vec, int len, ElemType num) {
    Vector res = initVector(len);
    for (int i = 0; i 
2.8向量除以数值 
Vector vectorDivNum(const Vector vec, int len, ElemType num) {
    Vector res = initVector(len);
    for (int i = 0; i 
2.9向量的均方根 
ElemType vectorRootMeanSquare(const Vector vec, int len) {
    ElemType nor;
    for (int i = 0; i 
2.10打印向量 
void vectorPrint(const Vector vec, int len) {
    for (int i = 0; i 
3.矩阵操作 
3.1创建矩阵 
Matrix initMatrix(int row, int col) {
    Matrix mat = (ElemType**)malloc(sizeof(ElemType*)*row);
    if (mat == NULL) {
        printf("MATRIX CREATE ERRORn");
        return NULL;
    }
    for (int i = 0; i 
3.2销毁矩阵 
int destroyMatrix(Matrix mat, int row) {
    if (!mat) return 1;
    for (int i = 0; i 
3.3填充矩阵 
Matrix fillMatrix(Matrix mat, int row, int col, int type) {
    int i, j, k;
    switch (type) {
    case ZERO_MATRIX:
        for (i = 0; i 
3.4矩阵获取行 
Vector getRow(const Matrix mat, int row, int col, int loc) {
    if (loc >= row) {
        printf("GET ROW ERRORn");
        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矩阵获取列 
Vector getCol(const Matrix mat, int row, int col, int loc) {
    if (loc >= col) {
        printf("GET COLUMN ERRORn");
        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设置行 
int setRow(Matrix A, int row, int col, Vector V, int loc){
    if(loc<0 || loc>row){
        return 0;
    }
    for(int i=0;i 
3.7设置列 
int setCol(Matrix A, int row, int col, Vector V, int loc){
    if(loc<0 || loc>col){
        return 0;
    }
    for(int i=0;i 
3.8矩阵加法 
Matrix matrixAdd(const Matrix x, const Matrix y, int xrow, int xcol, int yrow, int ycol) {
    if (xrow != yrow || xcol != ycol) {
        printf("MATRIX ADD ERRORn");
        return NULL;
    }
    Matrix res = (ElemType**)malloc(sizeof(ElemType*)*xrow);
    for (int i = 0; i 
3.9矩阵减法 
Matrix matrixSubtract(const Matrix x, const Matrix y, int xrow, int xcol, int yrow, int ycol) {
    if (xrow != yrow || xcol != ycol) {
        printf("MATRIX SUBTRACT ERRORn");
        return NULL;
    }
    Matrix res = (ElemType**)malloc(sizeof(ElemType*)*xrow);
    for (int i = 0; i 
3.10矩阵乘以数值 
Matrix matrixMulNum(const Matrix mat, int row, int col, ElemType num) {
    Matrix res = (ElemType**)malloc(sizeof(ElemType*)*row);
    for (int i = 0; i 
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= 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;
}


Matrix matrixMultiply(const Matrix x, int xrow, int xcol, const Matrix y, int yrow, int ycol) {
    if (xcol != yrow) {
        printf("MATRIX MULTIPLY ERRORn");
        return NULL;
    }
    if (xcol < 10) {
        return _marixDotMatrixSmall(x, y, xrow, xcol, ycol);
    }
    else {
        return _marixDotMatrixBig(x, y, xrow, xcol, ycol);
    }
}
 
3.12矩阵转置 
Matrix matrixTranspose(const Matrix mat, int row, int col) {
    Matrix res = (ElemType**)malloc(sizeof(ElemType*)*col);
    for (int i = 0; i 
3.13矩阵复制 
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求行列式 
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= 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* ERRORn");
    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= 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 ERRORn");
        return NULL;
    }
    else {
        _getMatrixStart(mat, row, col, temp);
        for (int i = 0; i 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求每一列的均值 
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矩阵标准化？ 
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打印矩阵 
void matrixPrint(const Matrix mat, int row, int col) {
    for (int i = 0; i 
3.19对角化 
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;
}