matrixf.h

简单的最小二乘曲线拟合方法

#include<string.h>
#include <math.h>
#include <vector>
using namespace std;
typedef vector<vector<double> > MatrixTwo;
typedef vector<double> MatrixOne;

void matrixadd(double *p1,double *p2,double *p3,int n,int m)
{
	double *pend1;
	pend1=p1+n*m-1;
	for(;p1<=pend1;p1++,p2++,p3++)
	{
		(*p3)=(*p1)+(*p2);
		
    }

}


/*void matrixmulti(double *a,double *b,double *c,int n,int m,int p)
{
	for(int i=0;i<n;i++)
		for(int j=0;j<p;j++)
		{
			c[i*m+j]=0;
			for(int k=0;k<m;k++)
				c[i*m+j]=c[i*m+j]+a[i*m+k]*b[k*p+j];
		}

}
double matrixmulti(double *A,double *B,int N)
{
	double SumOne = 0.0;
    for(int i=0; i<N; i++)
    {
		SumOne+=(*(A+i))*(*(B+i));
    }
    return(SumOne);
}*/

void MatrixTrans(MatrixTwo &A,MatrixTwo &B,int M,int N)
{
	for (int i =0; i<M; i++)
    for (int j =0; j<N; j++)
    B[j][i]=A[i][j];
}


void MatrixMulti(MatrixTwo &A, MatrixTwo &B, MatrixTwo &C, int M, int N, int P)
{
	int i,j, k;
	for (i=0;i<M;i++)
		for (j=0;j<P;j++)
		{
			C[i][j]=0;
			for (k=0;k<N;k++)
				C[i][j]+=A[i][k]*B[k][j];
		}
}





void AntiMatrix(MatrixTwo &A,MatrixOne T,int M,int N)
{
	int i,j,s;
	double a,w,u;




	for(int k=0;k<=M-1;k++)//选主元
	{
		w=A[k][k];
		s=k;
		T[k]=k;
		for(i=k+1;i<=M-1;i++)
		{
			if(abs(A[i][k])>abs(w))
			{
				w=A[i][k];
				s=i;
			}
		}
		if(w==0)
		{
			printf("singular");
			exit(0);
		}
		if(s!=k)//换行
		{
			T[k]=s;
			for(j=0;j<=M-1;j++)
			{
				u=A[k][j];
				A[k][j]=A[s][j];
				A[s][j]=u;

			}
		}

		a=1/A[k][k];
		A[k][k]=a;
		for(i=0;i<=M-1;i++)//变换第K列的元素
		{
			if(i==k)
			{
			}
			else
			{
				A[i][k]=-1*a*A[i][k];
				//s1=A[i][k];//printf("%lf\n",A[i][k]);

			}

		}


		for(i=0;i<=M-1;i++)//变换除第k行k列的所有元素
		{
			if(i==k)
			{
			}
			else
			{
				for(j=0;j<=N-1;j++)
				{
					if(j==k)
					{
					}
					else
					{
						A[i][j]=A[i][j]+A[i][k]*A[k][j];//s2=A[i][j];
					}

				}
			}
		}

		for(j=0;j<=M-1;j++)//变换第K行的元素
		{
			if(j==k)
			{
			}
			else
			{
				A[k][j]=a*A[k][j];//s3=A[k][j];
			}
		}

		/*for(i = 0; i < M; i++)
	{
		for (j = 0; j < N; j++)
		{
			printf("%6.2lf ", A[i][j]);
		}
		printf("\n");
	}*/

	}

	for(k=M-1;k>=0;k--)//列交换
	{
		if(T[k]!=k)
		{
			s=T[k];
			for(i=0;i<=M-1;i++)
			{
				u=A[i][s];
				A[i][s]=A[i][k];
				A[i][k]=u;
			}
		}
	}

}