C语言实现常用的矩阵运算

#include <stdio.h> #include <stdlib.h> #include <math.h> #define MAT_LEGAL_CHECKING #define min(a, b) ((a) > (b) ? (b) : (a)) #define equal(a, b) ((a-b)<1e-20 && (a-b)>-(1e-20)) //定义一个矩阵结构体，包括矩阵的行，列以及元素 typedef struct { int row, col; float **element; }Mat; /************************************************************************/ /* Private Function */ /************************************************************************/ void swap(int *a, int *b) { int m; m = *a; *a = *b; *b = m; } void perm(int list[], int k, int m, int* p, Mat* mat, float* det) { int i; if (k > m){ float res = mat->element[0][list[0]]; for (i = 1; i < mat->row; i++){ res *= mat->element[i][list[i]]; } if (*p % 2){ //odd is negative *det -= res; } else{ //even is positive *det += res; } } else{ // if the element is 0, we don't need to calculate the value for this permutation if (!equal(mat->element[k][list[k]], 0.0f)) perm(list, k + 1, m, p, mat, det); for (i = k + 1; i <= m; i++) { if (equal(mat->element[k][list[i]], 0.0f)) continue; swap(&list[k], &list[i]); *p += 1; perm(list, k + 1, m, p, mat, det); swap(&list[k], &list[i]); *p -= 1; } } } /************************************************************************/ /* Public Function */ /************************************************************************/ //创建一个矩阵 Mat* MatCreate(Mat* mat, int row, int col) { int i; mat->element = (float**)malloc(row * sizeof(float*)); if (mat->element == NULL){ printf("mat create fail!\n"); return NULL; } for (i = 0; i < row; i++){ mat->element[i] = (float*)malloc(col * sizeof(float)); if (mat->element[i] == NULL){ int j; printf("mat create fail!\n"); for (j = 0; j < i; j++) free(mat->element[j]); free(mat->element); return NULL; } } mat->row = row; mat->col = col; return mat; } //删除一个矩阵 void MatDelete(Mat* mat) { int i; for (i = 0; i<mat->row; i++) free(mat->element[i]); free(mat->element); } //把向量赋值给矩阵赋值 Mat* MatSetVal(Mat* mat, float* val) { int row, col; for (row = 0; row < mat->row; row++){ for (col = 0; col < mat->col; col++){ mat->element[row][col] = val[col + row * mat->col]; } } return mat; } //矩阵的输出函数 void MatDump(Mat* mat) { int row, col; #ifdef MAT_LEGAL_CHECKING if (mat == NULL){ return; } #endif printf("Mat %dx%d:\n", mat->row, mat->col); for (row = 0; row < mat->row; row++){ for (col = 0; col < mat->col; col++){ printf("%.4f\t", mat->element[row][col]); } printf("\n"); } } //给矩阵全部赋0 Mat* MatZeros(Mat* mat) { int row, col; for (row = 0; row < mat->row; row++){ for (col = 0; col < mat->col; col++){ mat->element[row][col] = 0.0f; } } return mat; } //给矩阵全部赋1 Mat* MatEye(Mat* mat) { int i; MatZeros(mat); for (i = 0; i < min(mat->row, mat->col); i++){ mat->element[i][i] = 1.0f; } return mat; } //两个矩阵的和运算 Mat* MatAdd(Mat* src1, Mat* src2, Mat* dst) { int row, col; #ifdef MAT_LEGAL_CHECKING if (!(src1->row == src2->row && src2->row == dst->row && src1->col == src2->col && src2->col == dst->col)){ printf("err check, unmatch Matix for MatAdd\n"); MatDump(src1); MatDump(src2); MatDump(dst); return NULL; } #endif for (row = 0; row < src1->row; row++){ for (col = 0; col < src1->col; col++){ dst->element[row][col] = src1->element[row][col] + src2->element[row][col]; } } return dst; } //两个矩阵的减运算 Mat* MatSub(Mat* src1, Mat* src2, Mat* dst) { int row, col; #ifdef MAT_LEGAL_CHECKING if (!(src1->row == src2->row && src2->row == dst->row && src1->col == src2->col && src2->col == dst->col)){ printf("err check, unmatch Matix for MatSub\n"); MatDump(src1); MatDump(src2); MatDump(dst); return NULL; } #endif for (row = 0; row < src1->row; row++){ for (col = 0; col < src1->col; col++){ dst->element[row][col] = src1->element[row][col] - src2->element[row][col]; } } return dst; } //矩阵的乘法运算 Mat* MatMul(Mat* src1, Mat* src2, Mat* dst) { int row, col; int i; float temp; #ifdef MAT_LEGAL_CHECKING if (src1->col != src2->row || src1->row != dst->row || src2->col != dst->col){ printf("err check, unmatch Matix for MatMul\n"); MatDump(src1); MatDump(src2); MatDump(dst); return NULL; } #endif for (row = 0; row < dst->row; row++){ for (col = 0; col < dst->col; col++){ temp = 0.0f; for (i = 0; i < src1->col; i++){ temp += src1->element[row][i] * src2->element[i][col]; } dst->element[row][col] = temp; } } return dst; } //矩阵的转置运算 Mat* MatTrans(Mat* src, Mat* dst) { int row, col; #ifdef MAT_LEGAL_CHECKING if (src->row != dst->col || src->col != dst->row){ printf("err check, unmatch Matix for MatTranspose\n"); MatDump(src); MatDump(dst); return NULL; } #endif for (row = 0; row < dst->row; row++){ for (col = 0; col < dst->col; col++){ dst->element[row][col] = src->element[col][row]; } } return dst; } //返回矩阵的行列式|A| float MatDet(Mat* mat) { float det = 0.0f; int plarity = 0; int i; #ifdef MAT_LEGAL_CHECKING if (mat->row != mat->col){ printf("err check, not a square Matix for MatDetermine\n"); MatDump(mat); return 0.0f; } #endif int *list; list = (int*)malloc(sizeof(int)*mat->col); if (list == NULL){ printf("malloc list fail\n"); return 0.0f; } for (i = 0; i < mat->col; i++) list[i] = i; perm(list, 0, mat->row - 1, &plarity, mat, &det); free(list); return det; } // dst = adj(src)伴随矩阵 Mat* MatAdj(Mat* src, Mat* dst) { Mat smat; int row, col; int i, j, r, c; float det; #ifdef MAT_LEGAL_CHECKING if (src->row != src->col || src->row != dst->row || src->col != dst->col){ printf("err check, not a square Matix for MatAdj\n"); MatDump(src); MatDump(dst); return NULL; } #endif MatCreate(&smat, src->row - 1, src->col - 1); for (row = 0; row < src->row; row++){ for (col = 0; col < src->col; col++){ r = 0; for (i = 0; i < src->row; i++){ if (i == row) continue; c = 0; for (j = 0; j < src->col; j++){ if (j == col) continue; smat.element[r][c] = src->element[i][j]; c++; } r++; } det = MatDet(&smat); if ((row + col) % 2) det = -det; dst->element[col][row] = det; } } MatDelete(&smat); return dst; } //矩阵的逆运算 /* Mat* MatInv(Mat* src, Mat* dst) { Mat adj_mat; float det; int row, col; #ifdef MAT_LEGAL_CHECKING if (src->row != src->col || src->row != dst->row || src->col != dst->col){ printf("err check, not a square Matix for MatInv\n"); MatDump(src); MatDump(dst); return NULL; } #endif MatCreate(&adj_mat, src->row, src->col); MatAdj(src, &adj_mat); det = MatDet(src); if (equal(det, 0.0f)){ printf("err, determinate is 0 for MatInv\n"); return NULL; } for (row = 0; row < src->row; row++){ for (col = 0; col < src->col; col++) dst->element[row][col] = adj_mat.element[row][col] / det; } MatDelete(&adj_mat); return dst; } */ Mat* MatInv(Mat* a, int n){ int *is = new int[n]; int *js = new int[n]; int i, j, k; float d, p; for (k = 0; k < n; k++) { d = 0.0; for (i = k; i <= n - 1; i++) for (j = k; j <= n - 1; j++) { p = abs(a->element[i][j]); if (p>d) { d = p; is[k] = i; js[k] = j; } } if (0.0 == d) { free(is); free(js); printf("err**not inv\n"); return(0); } if (is[k] != k) for (j = 0; j <= n - 1; j++) { p = a->element[k][j]; a->element[k][j] = a->element[is[k]][j];