gongetidu.cpp

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

#include <iostream.h>
#include <assert.h>
#include <math.h>



void gongetidu(double **a,double *b,double * x,const int n,const double eps)
{
	int i,j;
	int num;
	double alphak,betak,add,err,h,l,t;
	add=alphak=betak=0;
	
	double *p=new double [n];
	assert(p!=NULL);
	double *r=new double [n];
	assert(r!=NULL);
	
	for (i=0;i<n;i++)
	{
		p[i]=r[i]=b[i];
	}
	num=1;
	do 
	{
		//计算ALPHAK；
		h=l=0;
		for (i=0;i<n;i++)
		{
			t=0;
			for (j=0;j<n;j++)
			{
				t+=a[i][j]*p[j];
			}
			l+=p[i]*t;
			h+=r[i]*p[i];
		}
		alphak=h/l;
		/////////////////////////////////////////////		
		for (i=0;i<n;i++)
		{
			x[i]+=alphak*p[i];
		}
		
		
		for (i=0;i<n;i++)
		{
			h=0;
			for (j=0;j<n;j++)
			{
				h+=a[i][j]*x[j];
			}
			r[i]=b[i]-h;
		}
		//计算BETAK；
		h=l=0;
		for (i=0;i<n;i++)
		{	
			t=0;
			for (j=0;j<n;j++)
			{
				t+=a[i][j]*p[j];
			}
			h+=r[i]*t;
			l+=p[i]*t;
		}
		betak=-h/l;
		
		for (i=0;i<n;i++)
		{
			p[i]=r[i]+betak*p[i];
		}
		
		
		//计算误差
		err=0;
		for(i=0;i<n;i++) err+=r[i]*r[i];
		// 		if(num%10==0)
		{cout<<num<<"\t"<<sqrt(err)<<endl;}
		num++;
	} while(sqrt(err)>eps);
	
	delete [] p;
	p=NULL;
	delete [] r;
	r=NULL;
}


void ssorpcg(double **a, double *b, double *x,const int n,const double omiga,const double eps)
{
	int i,j,num=0;
	double delta,add,add2,add3,alphak,betak,temp;
	delta=add=add2=add3=0;


	double *r=new double[n];
	assert(r!=NULL);
	double *h=new double[n];
	assert(h!=NULL);
	double *p=new double[n];
	assert(p!=NULL);

	for (i=0;i<n;i++)
	{
		r[i]=h[i]=p[i]=0;
	}
// x[]为零向量

//r=b-Ax0
	for (i=0;i<n;i++)
	{
		r[i]=b[i];
	}
//Sovle M*h=r
	 h[0]=r[0];
	 for (i=1;i<n;i++)
	 {
		 add=0;
		 for (j=0;j<=i-1;j++)
		 {
			 add+=(omiga*a[i][j]/a[j][j])*h[j];
		 }
		 h[i]=r[i]-add;
	 }
	 h[n-1]=h[n-1]/a[n-1][n-1];
	 for (i=n-2;i>=0;i--)
	 {
		 add=0;
		 for (j=i+1;j<n;j++)
		 {
			 add+=-omiga*a[i][j]*h[j];
		 }
		 h[i]=(h[i]-add)/a[i][i];
	 }


//p=h
	 for (i=0;i<n;i++)
	 {
		 p[i]=h[i];
	 }

//|r|^2
	 delta=0;
	 for (i=0;i<n;i++)
	 {
		 delta+=r[i]*r[i];
	 }

while(sqrt(delta)>eps)
	 {
		 cout<<num++<<"\t"<<sqrt(delta)<<endl;

// ji suan AlphaK
		 add=add2=0;
		 for (i=0;i<n;i++)
		 {
			 add3=0;
			 for (j=0;j<n;j++)
			 {
				 add3+=a[i][j]*p[j];
			 }
			 add+=p[i]*add3;
			 add2+=r[i]*h[i];
		 }
		alphak=add2/add;


        temp=add2;

//x[k+1] and r[k+1]
		for (i=0;i<n;i++)
		{
			x[i]=x[i]+alphak*p[i];
			add=0;
			for (j=0;j<n;j++)
			{
				add+=a[i][j]*p[j];
			}
// 			cout<<add<<endl;
			r[i]=r[i]-alphak*add;
		}
		delta=0;
		for (i=0;i<n;i++)
		{
			delta+=r[i]*r[i];
		}

//h[k+1]=M^(-1)r[k+1]
		h[0]=r[0];
		for (i=1;i<n;i++)
		{
			add=0;
			for (j=0;j<=i-1;j++)
			{
				add+=omiga*a[i][j]/a[j][j]*h[j];
			}
			h[i]=r[i]-add;
		}
		h[n-1]=h[n-1]/a[n-1][n-1];
		for (i=n-2;i>=0;i--)
		{
			add=0;
			for (j=i+1;j<n;j++)
			{
				add+=-omiga*a[i][j]*h[j];
			}
			h[i]=(h[i]-add)/a[i][i];
		}

//betaK
		add=0;

		for (i=0;i<n;i++)
		{
			add=r[i]*h[i];
		}

		betak=add/temp;
//p(k+1)
		for (i=0;i<n;i++)
		{
			p[i]=h[i]+betak*p[i];
		}
//|r|^2
		delta=0;
		for (i=0;i<n;i++)
		{
			delta+=r[i]*r[i];
		}
	 }


	 delete [] r;
	 r=NULL;
	 delete [] h;
	 h=NULL;
	 delete [] p;
	 p=NULL;
}