leastsquare.cpp

本程序是我上课实验报告的源代码

#include"LeastSquare.h"
#include"stdio.h"


LeastSquare::LeastSquare()
{   ifstream infile("input.dat");
	    if(!infile)
		{
		cout << "Sorry! The input file failed!" << endl;
		exit(-1);
		}
	    else{
            infile>>m>>n>>N;//从文件读入行列数和非零元个数
	         int i,j;
             for(i=1; i<=m; i++) //将矩阵先初始化，全部置为零
			 { 
		         for(j=1;j<=n;j++)
	               {A[i][j]=0.;
	                b[i]=0.; 
	                }
			 }
		}
   
        for(int k=1; k<=N; k++)
		{	 
	          infile>>i>>j;
	          infile>>A[i][j];	//读入矩阵A的值
           }
          for(i=1; i<=m; i++)  
	        infile>>b[i];//读入右端项
	     
    infile.close(); 
}

LeastSquare::Function()
{

    for(i=1; i<=n; i++)   
	  for(j=1; j<=n; j++)//先将B[i][j]和G[i][i]初始化置零
         { B[i][j]=0;
           g[i]=0;
	       }	
	        
	  for(k=1; k<=m; k++)
	   {  
           for(i=1; i<=n; i++)
           for(j=1; j<=n; j++)  
		   {
	        B[i][j] += A[k][i] * A[k][j];  //计算正则矩阵B
		   }
	  }
            for(i=1; i<=n; i++)
				    for(k=1; k<=m; k++)
			{
            g[i] += A[k][i] * b[k];//计算右端项g
			}
	  
				
	//开始Cholesky分解
	if(B[1][1]<=0){
		cout << "Sorry! The data is wrong!" << endl;//首项不能小于等于零
		exit(-1);
	  }
	
	  
	B[1][1]=sqrt(B[1][1]);
	
	
	for(i=2;i<=n;i++)
		B[i][1]=B[i][1]/B[1][1];//先求第一列系数
		
		
	for(j=2;j<=n;j++)
	{        double s=0;
		    for(k=1;k<j;k++)
		      {
			     s=B[j][k]*B[j][k];
				 B[j][j]=B[j][j]-s;//求对角元素
				 
		       }
				
		      if(B[j][j]<=0)//对角元不能为小于等于零
		       {  cout << "Sorry! The data gets wrong!" << endl;
			        exit (-1);
		        }
		    
		     B[j][j]=sqrt(B[j][j]);
		     
		       for(i=j+1;i<=n;i++)
			   {     s=0;
			        for(k=1;k<j;k++)
			           {
				         s+=B[i][k]*B[j][k];

			            }
					 B[i][j]=B[i][j]-s;
					 
			         B[i][j]=B[i][j]/B[j][j];//得到某一对角元下面的一列元素
		        }
	}
	
	
	//分解结束,一下为回代求解
	    	for(i=1;i<=n;i++)
	    	    for(j=i+1;j<=n;j++)
				        B[i][j]=0;//将上三角元置零，得下三角阵L
				  
	
			g[1]=g[1]/B[1][1];//g[1]即为x1
			double temp;

			for( i=2;i<=n;i++)
				{for( j=1;j<=i-1;j++)
					{temp=0.0;
         			temp=temp+B[i][j]*g[j];
					}
         
				   g[i]=(g[i]-temp)/B[i][i];
				}

			for(j=1;j<=n;j++)
	    	    for(i=j+1;i<=n;i++)
				     B[j][i]=B[i][j];//转置L以求上三角方程组

			for(j=1;j<=n;j++)
	    	    for(i=j+1;i<=n;i++)
				        B[i][j]=0;//将下三角的元素置零

			g[n]=g[n]/B[n][n];

			for(int i=n-1;i>=1;i--)
				{for(int j=i+1;j<=n;j++)
					{temp=0.0;
					temp+=B[i][j]*g[j];
					}
				g[i]=(g[i]-temp)/B[i][i];
				}
			
			
			//回代结束，得到最小二乘解存于g[i]中

				
          SavetoFile();
		  return 1;
}

int LeastSquare::SavetoFile()
{
	ofstream outfile("results.dat");
	if(!outfile)
		{
		cout << "Sorry.The output failed!" << endl;
		return 0;
	  }
	 else{
	 	    for(int i=1; i<=n; i++)
	         outfile<<g[i]<<endl;
	        
        }outfile.close();
	return 1;
}



int main()
{   
  LeastSquare cholesky;
  cholesky.Function();
  cholesky.SavetoFile();
  getchar();
  return 1;
}