new.cpp

利用平方根法（Cholesky法）求超定方程组的最小二乘解

// New.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include "math.h"


//求a矩阵的转置；
//其中输入矩阵a为m×n阶，结果保存在n×m阶矩阵b里；
void transpose(double a[], int m, int n, double b[])
{
	int i,j;
	for (i=0; i<n; i++)
		for(j=0; j<m; j++)
		{
//			s = a[j*n+i];
//			b[i*m+j] = s;
			b[i*m+j] = a[j*n+i];
		}
}



//矩阵相乘；
//m×n阶的矩阵a和n×k阶的矩阵b相乘，得到m×k阶的矩阵保存在c中；
void trmul(double a[], double b[], int m, int n, int k, double c[])
{	
	int i,j,l,u;
	for (i=0; i<=m-1; i++)
		for (j=0; j<=k-1; j++)
		{
			u=i*k+j; 
			c[u]=0;
			for (l=0; l<=n-1; l++)
				c[u]=c[u]+a[i*n+l]*b[l*k+j];
		}
		return;
}


//平方根法（Cholesky法）求实正定对称方程组的解；
//输入n×n阶正定对称系数矩阵a，n×m阶常数矩阵；
//返回分解后的U矩阵上三角部分存于a，方程的m组解向量存于d；
int chlk(double a[], int n, int m, double d[])
{
	int i,j,k,u,v;
	if ((a[0]+1.0==1.0)||(a[0]<0.0))
//	if (a<0)
	{ 
		printf("fail\n"); 
		return(false);
	}
	a[0]=sqrt(a[0]);
	for (j=1; j<=n-1; j++) a[j]=a[j]/a[0];
	for (i=1; i<=n-1; i++)
	{
		u=i*n+i;
		for (j=1; j<=i; j++)
		{ 
			v=(j-1)*n+i;
			a[u]=a[u]-a[v]*a[v];
		}
		if ((a[u]+1.0==1.0)||(a[u]<0.0))
		{ 
			printf("fail\n"); 
			return(false);
		}
		a[u]=sqrt(a[u]);
		if (i!=(n-1))
		{ 
			for (j=i+1; j<=n-1; j++)
			{ 
				v=i*n+j;
				for (k=1; k<=i; k++)
					a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
					a[v]=a[v]/a[u];
			}
		}
	}
	for (j=0; j<=m-1; j++)
	{ 
		d[j]=d[j]/a[0];
		for (i=1; i<=n-1; i++)
		{ 
			u=i*n+i; 
			v=i*m+j;
			for (k=1; k<=i; k++)
				d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
				d[v]=d[v]/a[u];
		}
	}
	for (j=0; j<=m-1; j++)
	{
		u=(n-1)*m+j;
		d[u]=d[u]/a[n*n-1];
		for (k=n-1; k>=1; k--)
		{ 
			u=(k-1)*m+j;
			for (i=k; i<=n-1; i++)
			{ 
				v=(k-1)*n+i;
				d[u]=d[u]-a[v]*d[i*m+j];
			}
			v=(k-1)*n+k-1;
			d[u]=d[u]/a[v];
		}
	}
	return(true);
}


//求超定方程的最小二乘解
//输入m×n阶系数矩阵r，m×1阶常数矩阵p；
//返回值为n×1阶解向量；
float *fun(float r[], float p[], int m, int n)
{
    int i,j;
	double *R = new double [m*n];
	double *P = new double [m*1];
	double *RT = new double [n*m];
	double *a = new double [n*n];
	double *b = new double [m*1];
	float *Q = new float [n*1];

	for (i=0; i<m*n; i++)
	{
		R[i] = r[i];
		P[i] = p[i];
	}
	//for (j=0; j<m; j++)	
	//	P[j] = p[j];
	
    transpose(R, m, n, RT);
	trmul(RT, R, n, m, n, a);
	trmul(RT, P, n, m, 1, b);
    int flag = chlk(a, n, 1, b);

	if (!flag)
		printf("Error!");
	else
		for (i=0; i<n; i++)
			Q[i] = (float)b[i];
	return Q;

	delete []R;
	delete []P;
	delete []RT;
	delete []a;
	delete []b;
	delete []Q;
}


int _tmain(int argc, _TCHAR* argv[])
{
	using namespace std;


	int i,j,m,n,k;
	//cin>>m>>n;
	m = 4;
	n = 3;
	int s=m*n;
	int t=m*1;

	float a[12] = {1,1,-1,2,1,0,1,-1,0,-1,2,1};
	float b[4] = {2,-3,1,4};
    float *p;
	p = fun(a, b, m, n);

	for (i=0; i<n; i++)
	{
		cout<<p[i]<<endl;
	}




	//double *p = new double [s];
	//double *q = new double [t];

	//int *p = new int [s];
	//int *q = new int [t];

	//for (i=0; i<s; i++)
	//	cin>>p[i];
	//for (i=0; i<t; i++)
	//	cin>>q[i];

	//float a[16] = {5.0,7.0,6.0,5.0,7.0,10.0,8.0,7.0,6.0,8.0,10.0,9.0,5.0,7.0,9.0,10.0};
	//float b[8] = {23.0,92.0,32.0,128.0,33.0,132.0,31.0,124.0};

	//for (j=0; j<s; j++)
	//{
	//	cout<<a[j]<<"	";
	//	p[j] = a[j];
	//}
	//cout<<endl;
	//for (j=0; j<t; j++)
	//{
	//	cout<<b[j]<<"	";
	//	q[j] = b[j];
	//}
	//cout<<endl;

	//double *p,*q;
	//p = (double *)a;
	//q = (double *)b;

	//transpose(p, m, n, q);
	//trmul(p, q, m, n, k, r);

	//for (j=0; j<s; j++)
	//	cout<<p[j]<<"	";
	//cout<<endl;
	//for (j=0; j<t; j++)
	//	cout<<q[j]<<"	";
	//cout<<endl;

	//for (j=0; j<m*k; j++)
	//	cout<<r[j]<<"	";
	//cout<<endl;
	

//	int flag = chlk(p, m, n, q);
	
//	if (flag)
//		for (i=0; i<m; i++)
//			printf("x(%d)=%f,   %f\n",i,q[i*n],q[i*n+1]);



//	delete []p;
//	delete []q;
//	delete []r;
	return 0;
}