新建 文本文档 (2).txt

一些解决线形方程组的并行计算的程序,希望对大家有一点帮助

#include <iostream>
#include <stdlib.h>
#include <math.h>
using namespace std;
double *allocMem(int);
int GaussLineMain(double*,double *,double*,int);
int Jacobi(double*,double*,double*,double*,int,int);
int matrix(double matrixnum[],int JacobiTime,int matrixNum,int n)
{
	//int matrixNum;//矩阵的行数
	double *matrixA;//矩阵A ，初始系数矩阵
	double *matrixD;//矩阵D为A中的主对角阵
	double *matrixL;//矩阵L为A中的下三角阵
	double *matrixU;//矩阵U为A中的上三角阵
	double *B;      //为雅可比方法迭代矩阵
	double *f;      //为中间过渡矩阵
	double *x;      //为一维数组，存放结果
	double *xk;     //为一维数组，用来在迭代中使用
	double *b;      //为一维数组，存放方程组右边的系数
	
	
	int i,j,k;
//	cout<<"请输入矩阵的行数"<<endl;
//	cin>>matrixNum;
    matrixA=allocMem(matrixNum*matrixNum);
    matrixD=allocMem(matrixNum*matrixNum);
    matrixL=allocMem(matrixNum*matrixNum);
    matrixU=allocMem(matrixNum*matrixNum);
    B=allocMem(matrixNum*matrixNum);
    f=allocMem(matrixNum);
    x=allocMem(matrixNum);
    xk=allocMem(matrixNum);
    b=allocMem(matrixNum);
    
    //cout<<endl<<endl<<endl<<"<<请输入矩阵中各元素的值（为"<<matrixNum<<"*"<<matrixNum<<",共计"<<matrixNum*matrixNum<<"个元素)"<<">>:"<<endl<<endl;
////    for(i=0;i<matrixNum;i++)
////    {
////    	//cout<<"请输入矩阵中第"<<i+1<<"行的"<<matrixNum<<"个元素:";
////    	for(j=0;j<matrixNum;j++)
////    	  cin>>*(matrixA+i*matrixNum+j);
////    }
for(int h=0; h<9;h++)
{
	*(matrixA+h)=matrixnum[h];
	//cout<<*(matrixA+h)<<endl;
}
    cout<<"<<请输入方程组右边系数各元素值，共记"<<matrixNum<<"个"<<">>:"<<endl<<endl;
    for(i=0;i<matrixNum;i++)
       cin>>*(b+i);
       //下面将Ａ分裂为Ａ＝Ｄ－Ｌ－Ｕ
       
       for(i=0;i<matrixNum;i++)
         for(j=0;j<matrixNum;j++)
         *(matrixD+i*matrixNum+j)=*(matrixL+i*matrixNum+j)=*(matrixU+i*matrixNum+j)=0;
    	//cout<<"right"<<endl;
   for(i=0;i<matrixNum;i++)
     for(j=0;j<matrixNum;j++)
     {
      if(i==j&&*(matrixA+i*matrixNum+j))
      *(matrixD+i*matrixNum+j)=1/(*(matrixA+i*matrixNum+j));
      else if(i>j)*(matrixL+i*matrixNum+j)=-*(matrixA+i*matrixNum+j);
      else *(matrixU+i*matrixNum+j)=-*(matrixA+i*matrixNum+j);
      //求Ｂ矩阵中的元素
     }
      //cout<<"right"<<endl;
     for(i=0;i<matrixNum;i++)
      for(j=0;j<matrixNum;j++)
      {
      	//cout<<"right"<<endl;
      	double temp=0;
      	for(k=0;k<matrixNum;k++)
      	temp+=*(matrixD+i*matrixNum+k)*(*(matrixL+k*matrixNum+j)+*(matrixU+k*matrixNum+j));
      	*(B+i*matrixNum+j)=temp;
      	//cout<<*(B+i*matrixNum+j)<<endl;
      }
     //求f中的元素
     //cout<<"here"<<endl;
      for(i=0;i<matrixNum;i++)
      {
      	double temp=0;
      	for(j=0;j<matrixNum;j++)
      	 temp+=*(matrixD+i*matrixNum+j)*(*(b+j));
      	 *(f+i)=temp;
      	 
      	}
      	//cout<<"here"<<endl;
      	//计算Ｘ的初始向量值
      	GaussLineMain(matrixA,x,b,matrixNum);
      	//cout<<"here"<<endl;
      	//利用雅可比迭代公式求解ＸＫ的值
      	//int JacobiTime;
      	//cout<<endl<<endl<<endl<<"<<雅可比迭代开始，请输入希望迭代的次数>>:"<<endl;
      	//cin>>JacobiTime;
      	while(JacobiTime<=0)
      	{
      		cout<<"迭代次数必须大于０，请重新输入：";
      		cin>>JacobiTime;
      		
      	}
      	Jacobi(x,xk,B,f,matrixNum,JacobiTime);
      	//输出线形方程组的解
     	cout<<endl<<endl<<endl<<"<<方程组运算结果如下>>"<<endl;
      	cout.precision(20);
      	for(i=0;i<matrixNum;i++)
      	  cout<<"x"<<i+1+n*matrixNum<<"="<<*(xk+i)<<endl;
  
//   	  cout<<endl<<endl<<"求解过程结束..."<<endl<<endl;
    	  
     	  delete[] matrixA;
    	  delete[] matrixD;
     	  delete[] matrixL;
      	  delete[] matrixU;
     	  delete[] B;
    	  delete[] f;
     	  delete[] x;
     	  delete[] xk;
     	  delete[] b;
      	   
      	   return 0;
      	        	
      
}
////////////////////////////////////////////////////////
//分配内存空间函数


double *allocMem(int num)
{
	double *head;
	if((head=new double[num])==NULL)
	{
		cout<<"内存空间分配失败，程序终止运行～！"<<endl;
		exit(0);
	}
	return head;
}
////////////////////////////////////////////////////////////////
//计算Ａx=b中Ｘ的初始向量值，采用高斯列诸元素消去法

int GaussLineMain(double *A,double*x,double*b,int num)
{
	int i,j,k;
	for(i=0;i<num-1;i++)
	{
		//首先每列选主元，最大的一个
		double lineMax=fabs(*(A+i*num+i));
		int lineI=i;
		for(j=i;j<num;j++)
		  if(fabs(*(A+j*num+i))>fabs(lineMax))lineI=j;
		  //主元所在的行和当前处理行做行交换，右系数b也随之交换
		  for(j=i;j<num;j++)
		  {
		  	lineMax=*(A+i*num+j);
		  	*(A+i*num+j)=*(A+lineI*num+j);
		  	*(A+lineI*num+j)=lineMax;
		  	lineMax=*(b+i);
		  	*(b+i)=*(b+lineI);
		  	*(b+lineI)=lineMax;
		  	
		  }
		if(*(A+i*num+i)==0) continue;
		//　将Ａ变为上三角矩阵，同样b也随之变换
		for(j=i+1;j<num;j++)
		{
			double temp=-*(A+j*num+i)/(*(A+i*num+i));
			for(k=i;k<num;k++)
			{
				*(A+j*num+k)+=temp*(*(A+i*num+k));
				
				
			}
			*(b+j)+=temp*(*(b+i));
		}
		
		}
	//验证Ａx=b是否有唯一解，就是验证Ａ的行列式是否为０；
	double determinantA=1;	
	for(i=0;i<num;i++)
	  determinantA*=*(A+i*num+i);
	  if(determinantA==0)
	  {
	  	cout<<endl<<endl<<"通过计算，矩阵Ａ的行列式为｜Ａ｜＝０，既Ａ没有唯一解\n　　程序退出..."<<endl<<endl;
	  	exit(0);
	  }
	  //从最后一行开始，回代求解Ｘ的初始向量
	 for(i=num-1;i>=0;i--)
	 {
	  for(j=num-1;j>i;j--)
	    *(b+i)-=*(A+i*num+j)*(*(x+j));
	    *(x+i)=*(b+i)/(*(A+i*num+i));
	 }
	 return 0;
}
//////////////////////////////////////////////////////////////////////////
//利用雅可比迭代公式求解Ｘ的精确值
int Jacobi(double *x,double *xk,double *B,double *f,int num,int time)
  {
  	int t=1,i,j;
  	while(t<=time)
  	{
  		for(i=0;i<num;i++)
  		{
  			double temp=0;
  			for(j=0;j<num;j++)
  			 temp+=*(B+i*num+j)*(*(x+j));
  			 *(xk+i)=temp+*(f+i);
  		}
  		for(i=0;i<num;i++)
  		 *(x+i)=*(xk+i);
  		 t++;
  	}
  	return 0;
  }
  int main()
{
	int m,Linenum,jacoTime;
	int n=0;
	
	cout<<"请输入行数：（行数与列数相等 既是个方阵）"<<endl;
	cin>>Linenum;
	double matrixValue[Linenum*Linenum];
	cout<<"<<请输入矩阵中各元素的值（为"<<Linenum<<"*"<<Linenum<<",共计"<<Linenum*Linenum<<"个元素)"<<">>:"<<endl<<endl;
	for(int j=0;j<Linenum;j++)
	{
	 cout<<"请输入矩阵中第"<<j+1<<"行的"<<Linenum<<"个元素:";
	 for(int h=0;h<Linenum;h++)
	 cin>>matrixValue[j*Linenum+h];
	}
	cout<<"请输入你要迭代的次数："<<endl;
	cin>>jacoTime;
	for(m=0;m<3;m++)
	{  
		matrix(matrixValue,jacoTime,Linenum,n);
	    n++;
	}
	cout<<endl<<endl<<"求解过程结束..."<<endl<<endl;
	return 0;
}