共轭梯度法.cpp

共轭梯度法

#include<stdio.h>
#include<stdlib.h>
#include<malloc.h>
#include<iostream.h>
#include<math.h>

void main(){
    double *A;
	double *b;
	double *x;
	double *r;
	double *v;
	double *z;
	double t,c,d,sum,p,q;
	int M,n,i,j,k;
	char tag;

L:	cout<<"n=";
	cin>>n;
	
	//initial x,A,b,r,v,z
    A=(double*)malloc((n*n+1)*sizeof(double));
	x=(double*)malloc((n+1)*sizeof(double));
	b=(double*)malloc((n+1)*sizeof(double));
	r=(double*)malloc((n+1)*sizeof(double));
	v=(double*)malloc((n+1)*sizeof(double));
	z=(double*)malloc((n+1)*sizeof(double));

	//initial the values of A[i,j],b[i]
	for(i=1;i<=n;i++){

		sum=0;
		for(j=1;j<=n;j++){

			A[n*(i-1)+j]=pow((i+j+1),-1);
            sum=sum+A[n*(i-1)+j];

		}//for

		b[i]=sum/3;

	}//for

//output A,b
	cout<<"A is as follows:"<<endl;

	for(i=1;i<=n;i++){

		for(j=1;j<=n;j++) printf("%f  ",A[n*(i-1)+j]);

	cout<<endl;

	}//for
    
	cout<<endl<<"b is as follows:"<<endl;
    for(i=1;i<=n;i++) printf("%f\n",b[i]);

	//input initial value of x
	cout<<endl<<"Please input x as directions:"<<endl;
    for(i=1;i<=n;i++){

		printf("x[%d]=",i);
        cin>>x[i];

	}//for

	//r=b-Ax,v=r,c=<r,r>
	c=0;
	for(i=1;i<=n;i++){

		sum=0;
		for(j=1;j<=n;j++) sum=sum+A[n*(i-1)+j]*x[j];

		r[i]=b[i]-sum;
        v[i]=r[i];
        c=c+r[i]*r[i];

	}//for
    
	cout<<"M=";
	cin>>M;
/*	cout<<"delta=";
	cin>>delta;
	cout<<"e=";
	cin>>e;
*/
    for(k=1;k<=M;k++){

		p=0;
		for(i=1;i<=n;i++) p=p+v[i]*v[i];

		if(sqrt(p)<pow(0.00000001,n)) break;

        //z=Av,t=c/<v,z>
		q=0;
		for(i=1;i<=n;i++){

			z[i]=0;
			for(j=1;j<=n;j++) z[i]=z[i]+A[n*(i-1)+j]*v[j];
            
            q=q+v[i]*z[i];

		}//for
        t=c/q;
		
        //x=x+tv,r=r-tz,d=<r,r>
		d=0;
		for(i=1;i<=n;i++){
			x[i]=x[i]+t*v[i];
			r[i]=r[i]-t*z[i];
            d=d+r[i]*r[i];
		}//for

		if(d<=pow(0.00000001,n)) break;

		//v=r+(d/c)v,c=d
		for(i=1;i<=n;i++) v[i]=r[i]+(d/c)*v[i];
		c=d;

		//output k,x,r
		printf("step %d:\n",k);
        for(i=1;i<=n;i++){

			printf("x[%d]=%f\tr[%d]=%f\n",i,x[i],i,r[i]);

		}//for
		cout<<endl;
	}//for
	
    cout<<"Try again?<y/n>";
	cin>>tag;
	cout<<endl;

	if(tag=='y') goto L;

	free(A);
	free(b);
	free(v);
	free(x);
	free(z);
	free(r);

}//main