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向量初始化
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向量销毁
int destoryVector(Vector vec) {
free(vec);
return 1;
}
2.3向量相加
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向量相减
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)
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<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向量乘数值
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向量除以数值
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向量的均方根
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打印向量
void vectorPrint(const Vector vec, int len) {
for (int i = 0; i<len; i++) {
printf("%lf\t", vec[i]);
}
printf("\n");
}
3.矩阵操作
3.1创建矩阵
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销毁矩阵
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填充矩阵
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矩阵获取行
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矩阵获取列
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设置行
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设置列
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矩阵加法
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矩阵减法
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矩阵乘以数值
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;
}
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矩阵转置
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矩阵复制
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<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矩阵取逆
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;
}
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求每一列的均值
Matrix matrixMeanCol(const Matrix mat, int row, int col) {
int i, j;
Matrix result = initMatrix(1, col);
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<row; i++) {
for (int j = 0; j<col; j++) {
printf("%f\t", mat[i][j]);
}
printf("\n");
}
printf("\n");
}
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;
}