#include<stdio.h>
#include<math.h>

#define N 3

void cg(n, rep, a, x, b, eps)
int n, rep;
double a[][N], x[], b[];
double eps;
{
	int i, j, l;
	double alpha, beta;
	double Ap[N], r[N], p[N];
	double rp, pAp, rr, rrp;
	double xx, norm, normm;
	
	for(i=0; i<n; i++){
		r[i] = b[i];
		for(j=0; j<n; j++)
			r[i] -= a[i][j]*x[j];
		p[i] = r[i];
	}
	
	for(l=0; l<rep; l++){
		rp = 0.;
		pAp = 0.;
		rr = 0.;
		for(i=0; i<n; i++){
			Ap[i] = 0.;
			for(j=0; j<n; j++)
				Ap[i] += a[i][j]*p[j];
			rp += r[i]*p[i];
			pAp += p[i]*Ap[i];
			rr += r[i]*r[i];
		}
		alpha = rp/pAp;
		rrp = 0.;
		normm = 0.;
		for(i=0; i<n; i++){
			norm = x[i];
			x[i] += alpha*p[i];
			norm -= x[i];
			if(fabs(norm) > normm) normm = fabs(norm);
			r[i] -= alpha*Ap[i];
			rrp += r[i]*r[i];
		}
		if(normm < eps){
			printf("%i\n", l);
			printf("convergence!!\n");
			return;
		}
		beta = rrp/rr;
		for(i=0; i<n; i++)
			p[i] = r[i]+beta*p[i];
	}
	
	printf("no convergence!!\n");
	return;	
}

void main(void){
	int i, j;
	double A[N][N]={{1, 1, 0},{1, 2, 0},{0, 0, 1}};
	double x[N];
	double b[N]={2, 3, 1};
	double eps;

	eps = pow(10, -7);
	cg(N, 100, A, x, b, eps);
	
	for(i=0; i<N; i++)
		printf("%lf\n", x[i]);
		
	return;
}