宝藏就在这里,就看还在上这门课的学弟学妹能不能发现了
一、矩阵与矢量运算模块
1.1矢量运算
1.1.1矢量加法
1.原理
有两个三维空间矢量
我们可以得
2.代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 | bool VecAdd(const int n1, const double Vec1[], const int n2, const double Vec2[], double result[]) { int i = 0; if (n1!=n2) { printf("Vector's dimensions must agree\n"); return false; } for (i; i < n1;i++) { result[i] = Vec1[i] + Vec2[i]; } return true; } |
1.1.2矢量减法
1.原理
有两个三维空间矢量
我们可以得
2.代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 | bool VecSub(const int n1, const double Vec1[], const int n2, const double Vec2[], double result[]) { int i = 0; if (n1 != n2) { printf("Vector's dimensions must agree\n"); return false; } for (i; i < n1; i++) { result[i] = Vec1[i] - Vec2[i]; } return true; } |
1.1.3矢量叉积
1.原理
有两个三维空间矢量
我们可以得
2.代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | bool VecCrossPlus(const int n1, const double Vec1[], const int n2, const double Vec2[], double result[]) { int i = 0; if (n1 != n2) { printf("Vector's dimensions must agree\n"); return false; } result[0] = Vec1[1] * Vec2[2] - Vec2[1] * Vec1[2]; result[1] = Vec1[2] * Vec2[0] - Vec1[0] * Vec2[2]; result[2] = Vec1[0] * Vec2[1] - Vec2[0] * Vec1[1]; return true; } |
1.1.4矢量乘积
1.原理
有两个三维空间矢量
我们可以得
2.代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | bool VecPlus(const int n1, const double Vec1[], const int n2, const double Vec2[], double *result) { int i =0 ; *result = 0; if (n1!=n2) { printf("Vector's dimensions must agree\n"); return false; } for (i; i < n1;i++) { *result = Vec2[i] * Vec1[i] + *result; } return true; } |
1.2矩阵运算
1.2.2矩阵加法
1.原理
设矩阵
规定
2.代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | bool MatAdd(const int row1, const int col1, const int row2, const int col2, const double A[], const double B[], double output[]) { int i = 0; if (row1 != row2 || col1 != col2) { printf("MatAdd:Inner matrix dimensions must agree.\n"); return false; } else { for ( i = 0; i < row1*col1; i++) { output[i] = A[i] + B[i]; } return true; } } |
1.2.3矩阵减法
1.原理
设矩阵
规定
2.代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | bool MatSub(const int row1, const int col1, const int row2, const int col2, const double A[], const double B[], double output[]) { int i = 0; if (row1 != row2 || col1 != col2) { printf("MatAdd:Inner matrix dimensions must agree.\n"); return false; } else { for (int i = 0; i < row1*col1; i++) { output[i] = A[i] - B[i]; } return true; } } |
1.2.4矩阵数量乘法
1.原理
设矩阵 ,
是数域F中的任意一个数
规定
2.代码
1 2 3 4 5 6 7 8 9 10 | bool MatMul_Num(const double k, const int row, const int col,const double A[], double output[])//矩阵相乘 { int i = 0; for ( i = 0; i < row*col;i++) { output[i] = k*A[i]; } return true; } |
1.2.5矩阵乘法
1.原理
设矩阵
规定
2.代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 | bool MatMul_Mat(const int row1, const int col1, const int row2, const int col2, const double A[], const double B[], double output[]) { int i, j, l; if (col2!=row2) { printf("MatAdd:Inner matrix dimensions must agree.\n"); return false; } else {
int u = 0; for ( i = 0; i <= row1 - 1; i++) for ( j = 0; j <= col2 - 1; j++) { u = i*col2 + j; output[u] = 0.0; for ( l = 0; l <= col1 - 1; l++) { output[u] = output[u] + A[i*col1 + l] * B[l*col2 + j]; } } return true; } } |
1.2.6矩阵点乘
1.原理
设矩阵
规定
2.代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | bool MatMul_Poi(const int row1, const int col1, const int row2, const int col2, const double A[], const double B[], double output[]) { int i; if (row1!=row2||col1!=col2) { printf("MatAdd:Inner matrix dimensions must agree.\n"); return false; } else { for (i = 0; i < row1*col1; i++) { output[i] = A[i] * B[i]; } return ture; } } |
1.2.7矩阵求逆
1.原理
对于一个方阵 ,若有
,则称矩阵
为矩阵
的逆矩阵。
2.代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 | bool MatInver(int n, double a[], double b[]) { int i, j, k, l, u, v, is[10], js[10]; /* matrix dimension <= 10 */ double d, p;
if (n <= 0) { printf("Error dimension in MatrixInv!\n"); return false; }
/* 将输入矩阵赋值给输出矩阵b,下面对b矩阵求逆,a矩阵不变; */ for (i = 0; i < n; i++) { for (j = 0; j < n; j++) { b[i*n + j] = a[i*n + j]; } }
for (k = 0; k < n; k++) { d = 0.0; for (i = k; i < n; i++) /* 查找右下角方阵中主元素的位置 ;*/ { for (j = k; j<n; j++) { l = n*i + j; p = fabs(b[l]); if (p>d) { d = p; is[k] = i; js[k] = j; } } }
if (d < 2.2204460492503131e-16) /* 主元素接近于0,矩阵不可逆 */ { printf("Divided by 0 in MatrixInv!\n"); return false; }
if (is[k] != k) /* 对主元素所在的行与右下角方阵的首行进行调换; */ { for (j = 0; j < n; j++) { u = k*n + j; v = is[k] * n + j; p = b[u]; b[u] = b[v]; b[v] = p; } }
if (js[k] != k) /* 对主元素所在的列与右下角方阵的首列进行调换; */ { for (i = 0; i < n; i++) { u = i*n + k; v = i*n + js[k]; p = b[u]; b[u] = b[v]; b[v] = p; } }
l = k*n + k; b[l] = 1.0 / b[l]; /* 初等行变换; */ for (j = 0; j < n; j++) { if (j != k) { u = k*n + j; b[u] = b[u] * b[l]; } } for (i = 0; i < n; i++) { if (i != k) { for (j = 0; j < n; j++) { if (j != k) { u = i*n + j; b[u] = b[u] - b[i*n + k] * b[k*n + j]; } } } } for (i = 0; i < n; i++) { if (i != k) { u = i*n + k; b[u] = -b[u] * b[l]; } } }
for (k = n - 1; k >= 0; k--) /* 将上面的行列调换重新恢复; */ { if (js[k] != k) { for (j = 0; j < n; j++) { u = k*n + j; v = js[k] * n + j; p = b[u]; b[u] = b[v]; b[v] = p; } } if (is[k] != k) { for (i = 0; i < n; i++) { u = i*n + k; v = is[k] + i*n; p = b[u]; b[u] = b[v]; b[v] = p; } } }
return true; } |
1.2.8矩阵转置
1.原理
2.代码
1 2 3 4 5 6 7 8 9 10 11 12 | bool MatTrans(const int row, const int col, const double A[], double output[])//矩阵转置 { for (int i = 0; i < row; i++) { for (int j = 0; j < col;j++) { output[j*col + i] = A[i*col + j]; } } return true; } |