矩matrix八卦
下面这个文件是进行矩阵运算的一些操作，你使用最后一个函数EquationResolution就可以进行最小二乘求解了。当然，在计算前，你需要自己构造上述的系数矩阵和常数矩阵，其实就是两个数组。
// Matrix.hpp
#if !defined(_MATRIXCALCULATE_HPP__INCLUDED_)
#define _MATRIXCALCULATE_HPP__INCLUDED_
#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000
#include <math.h>
class CMatrixCalculate
{
public:
////求矩阵matrix的行列式值, n为维数
double CalculateLiner(double *matrix,int n)
{
double liner=0;
int i=0,j=0,k=0;
int p=0,q=0;
if(n==1)
return matrix[0];
else
{
double *tr=(double *)calloc((n-1)*(n-1),sizeof(double));
for(k=0;k<n;k++)
{
p=0;
for(i=0;i<n;i++)
{
if(i!=0)
{
q=0;
for(j=0;j<n;j++)
{
if(j!=k)
{
tr[p*(n-1)+q]=matrix[i*n+j];
q++;
}
}
p++;
}
}
liner+=matrix[k]*pow(-1,k)*CalculateLiner(tr,n-1);
}
free(tr);
return liner;
}
}
////求矩阵matrix的第i行,j列的代数余子式, n为维数
double CalculateCofacter(double *matrix,int i,int j,int n)
{
int x=0,y=0;
int p=0,q=0;
if(n==1)
return matrix[0];
else
{
double *tr=(double *)calloc((n-1)*(n-1),sizeof(double));
p=0;
for(x=0;x<n;x++)
{
if(x!=i)
{
q=0;
for(y=0;y<n;y++)
{
if(y!=j)
{
tr[p*(n-1)+q]=matrix[x*n+y];
q++;
}
}
p++;
}
}
double cc=pow(-1,i+j)*CalculateLiner(tr,n-1);
free(tr);
return cc;
}
}
////矩阵转置, matrixAT=(matrixA)T, m,n为matrixA的行,列数
void CalculateAT(double *matrixA,double *matrixAT,int m,int n)
{
for (int i=0;i<m;i++)
{
for (int j=0;j<n;j++)
{
matrixAT[j*m+i]=matrixA[i*n+j];
}
}
}
////矩阵相乘, matrixAB=matrixA*matrixB, i,j为matrixA的行,列数, j,k为为matrixB的行,列数
void CalculateAB(double *matrixA,double *matrixB,double *matrixAB,int i,int j,int k)
{
for (int m=0;m<i;m++)
{
for (int n=0;n<k;n++)
{
matrixAB[m*k+n]=0;
for (int s=0;s<j;s++)
{
matrixAB[m*k+n]+=matrixA[m*j+s]*matrixB[s*k+n];
}
}
}
}
////matrixATA=(matrixA)T*matrixA, m,n分别为matrixA的行,列数
void CalculateATA(double *matrixA,double *matrixATA,int m,int n)
{
double *at=(double *)calloc((m)*(n),sizeof(double));
CalculateAT(matrixA,at,m,n);
CalculateAB(at,matrixA,matrixATA,n,m,n);
free(at);
}
////matrixA_为matrixA的逆, m为阶数
void CalculateA_(double *matrixA,double *matrixA_,int m)
{
double liner=CalculateLiner(matrixA,m);
for(int i=0;i<m;i++)
{
for(int j=0;j<m;j++)
matrixA_[j*m+i]=CalculateCofacter(matrixA,i,j,m)/liner;
}
}
////////////////////////////////////////////////////////////////////
////求正定矩阵a的逆矩，n为阶数
int MatrixInversion(double *a, int n)
{
int i, j, k, m;
double w, g, *b;
b = new double [n];
for(k = 0; k <= n - 1; k++)
{
w = a[0];
w=a[0]+1.0e-15;
/*
if(fabs(w)+1.0 == 1.0)
{
delete b;
printf("fail\n");
return(-2);
}
*/
m = n - k - 1;
for(i = 1; i <= n - 1; i++)
{
g = a[i * n];
b[i] = g / w;
if(i <= m)
b[i] = -b[i];
for(j = 1; j <= i; j++)
a[(i - 1) * n + j - 1] = a[i * n + j] + g * b[j];
}
a[n * n - 1] = 1.0 / w;
for(i = 1; i <= n - 1; i++)
a[(n - 1) * n + i - 1] = b[i];
}
for(i = 0; i <= n - 2; i++)
for(j = i + 1; j <= n - 1; j++)
a[i * n + j] = a[j * n + i];
delete b;
return(2);
}
////求正定矩阵a的逆矩，n为阶数
void MatInversion(double *a,int n)
{
int i,j,k;
for(k=0;k<n;k++)
{
for(i=0;i<n;i++)
{
if(i!=k) *(a+i*n+k) = -*(a+i*n+k)/(*(a+k*n+k));
}
*(a+k*n+k)=1/(*(a+k*n+k));
for(i=0;i<n;i++)
{
if(i!=k)
{
for(j=0;j<n;j++)
{
if(j!=k) *(a+i*n+j) += *(a+k*n+j)* *(a+i*n+k);
}
}
}
for(j=0;j<n;j++)
{
if(j!=k) *(a+k*n+j)*=*(a+k*n+k);
}
}
}
/*矩阵求逆子程序(Good)*/
int Invers_matrix(double *m1,int n)
{
int *is,*js;
int i,j,k,l,u,v;
double temp,max_v;
is=(int *)malloc(n*sizeof(int));
js=(int *)malloc(n*sizeof(int));
if(is==NULL||js==NULL)
{
printf("out of memory!\n");
return(0);
}
for(k=0;k<n;k++)
{
max_v=0.0;
for(i=k;i<n;i++)
{
for(j=k;j<n;j++)
{
temp=fabs(m1[i*n+j]);
if(temp>max_v)
{
max_v=temp; is[k]=i; js[k]=j;
}
}
}
if(max_v==0.0)
{
free(is); free(js);
printf("invers is not availble!\n");
return(0);
}
if(is[k]!=k)
{
for(j=0;j<n;j++)
{
u=k*n+j; v=is[k]*n+j;
temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
}
}
if(js[k]!=k)
for(i=0;i<n;i++)
{
u=i*n+k; v=i*n+js[k];
temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
}
l=k*n+k;
m1[l]=1.0/m1[l];
for(j=0;j<n;j++)
{
if(j!=k)
{
u=k*n+j;
m1[u]*=m1[l];
}
}
for(i=0;i<n;i++)
{
if(i!=k)
{
for(j=0;j<n;j++)
{
if(j!=k)
{
u=i*n+j;
m1[u]-=m1[i*n+k]*m1[k*n+j];
}
}
}
}
for(i=0;i<n;i++)
{
if(i!=k)
{
u=i*n+k;
m1[u]*=-m1[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;
temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
}
if(is[k]!=k)
for(i=0;i<n;i++)
{
u=i*n+k; v=i*n+is[k];
temp=m1[u]; m1[u]=m1[v]; m1[v]=temp;
}
}
free(is); free(js);
return(1);
}
////求矩阵方程matrixA*matrixX=matrixL的最小二乘解, m,n分别为matrixA矩阵的行,列数
void EquationResolution(double *matrixA,double *matrixL,double *matrixX,int m,int n)
{
if (m<n) return;
double *at=(double *)malloc((m)*(n)*sizeof(double));
double *ata=(double *)malloc((n)*(n)*sizeof(double));
double *atl=(double *)malloc((n)*sizeof(double));
CalculateATA(matrixA,ata,m,n);
MatrixInversion(ata,n);
CalculateAT(matrixA,at,m,n);
CalculateAB(at,matrixL,atl,n,m,1);
CalculateAB(ata,atl,matrixX,n,n,1);
free(at);
free(ata);
free(atl);
}
}
#endif // !defined(_MATRIXCALCULATE_HPP__INCLUDED_)