main.cpp

这个程序利用共轭梯度方法

//用共轭梯度法求解一个线性方程组
#include <iostream>
#include <assert.h>
#include <math.h>
#include <fstream>
#include "gongetidu.h"
using namespace std;

const int kk=10;
const double eps=1e-10;
const double omiga=1.0;





void main()
{
	int i,j;//循环变量

	ofstream outf( "cg.txt", ios_base::out );
	if (!outf)
	{
		cerr<<"can not open file"<<endl;
		exit( -1 );
	}
	ofstream outf2( "pcg.txt", ios_base::out );
	if (!outf)
	{
		cerr<<"can not open file"<<endl;
		exit( -1 );
	}
	ofstream ot( "ot.txt", ios_base::out );
	if (!outf)
	{
		cerr<<"can not open file"<<endl;
		exit( -1 );
	}


	
	
	//为矩阵，右端项，解分配内存空间	
	double **a=new double*[kk*kk];
	assert(a!=0);
	for(i=0;i<kk*kk;i++)
	{
		a[i]=new double[kk*kk];
		assert(a[i]!=0);
		
		for (j=0;j<kk*kk;j++)
		{
			a[i][j]=0;
		}
	}
	double *b= new double[kk*kk];
	assert(b!=0);
	double *x= new double[kk*kk];
	assert(x!=0);
	
	
	
	//为数组赋值
	for(i=0;i<kk;i++)
	{
		
		for(int j=0;j<kk;j++)
		{
			if(j==0)
			{
				a[i*kk][i*kk]=4;a[i*kk][i*kk+1]=-1;
			}
			else if(j==kk-1)
			{
				a[i*kk+kk-1][i*kk+kk-1]=4;a[i*kk+kk-1][i*kk+kk-2]=-1;
			}
			else
			{
				a[i*kk+j][i*kk+j-1]=-1;
				a[i*kk+j][i*kk+j]=4;
				a[i*kk+j][i*kk+j+1]=-1;
			}
			if(i!=kk-1)
				a[i*kk+j][(i+1)*kk+j]=-1;
			if(i!=0)
				a[i*kk+j][(i-1)*kk+j]=-1;
		}
		//为右端项赋值，并使x[]归0.
	}
	for (i=0;i<kk*kk;i++)
	{
		b[i]=1;
		x[i]=0;
	}
	for (i=0;i<kk*kk;i++)
	{
		ot<<x[i]<<endl;
	}

	
	gongetidu(a,b,x,kk*kk,eps);
	for (i=0;i<kk*kk;i++)
	{
		outf<<x[i]<<endl;
	}


	for (i=0;i<kk*kk;i++)
	{
		b[i]=1;
		x[i]=0;
	}

 	ssorpcg(a,b,x,kk*kk,omiga,eps);

	for (i=0;i<kk*kk;i++)
	{
		outf2<<x[i]<<endl;
	}
	
	//释放内存空间	
	for(i=0;i<kk*kk;i++)
	{
		delete [] a[i];
		a[i]=NULL;
	}
	delete [] a;
	a=NULL;
	delete [] b;
	b=NULL;	
	delete [] x;
	x=NULL;
}





//迭代初试向量为零向量












 








double product(double *a,double *b,const int n)
{
	int i;
	double add=0;
	for (i=0;i<n;i++)
	{
		add+=a[i]*b[i];
	}
	return add;
}

double product(double *a,double *b,double **m,const int n)
{
	int i,j;
	double add,add2;
	add=0;
	for (i=0;i<n;i++)
	{
		add2=0;
		for (j=0;j<n;j++)
		{
			add2+=m[i][j]*b[j];
		}
		add+=add2*a[i];
	}
	return add;
}