#include <stdio.h>
#include <math.h>
/* cg.c : 1999 3 13 by Oguni */
#define N 4
main(){
  double A[N][N]={4.,-1.,0.,0.,-1.,4.,-1.,0.,0.,-1.,4.,-1.,0.,0.,-1.,4.};
  double b[N]={3.,2.,2.,3.}, x[N];
  double tol[2]={0.000001, 0.};
  int n=4, itr[2]={20, 1}, i;
  cg(A,b,x,tol,n,itr);
  printf("%d %f\n", itr[1], tol[1]);
  for (i=0; i<n; i++){
    printf("%f\n", x[i]); }
}
/* function */
cg(A,b,x,tol,n,itr)
double A[N][N], b[N], x[N], tol[2];
int n, itr[2];
{
  int k, i, ic;
  double g[N], d[N], 	q[N], qd[N], val[N];
  double u0, u, s, r, t, beta, vn, vb, bn, err;
  for (i=0; i<n; i++){
    g[i]=-b[i]; d[i]=b[i]; }
  u0=0.0;
  for (i=0; i<n; i++){
    u0+=b[i]*b[i]; }
  bn=sqrt(u0); k=0; u=u0; ic=1;
  while (ic==1){
    k++; itr[1]=k;
    mvsub(A,d,qd,n);
    for (i=0; i<n; i++){ q[i]=qd[i]; }
    r=0.0;
    for (i=0; i<n; i++){ r+=d[i]*q[i]; }
    if (r==0) printf("zero divide.\n");
    s=u/r;
    for (i=0; i<n; i++){
      x[i]+=s*d[i]; g[i]+=s*q[i]; }
    t=0.0;
    for (i=0; i<n; i++){
      t+=g[i]*q[i]; }
    beta=t/r;
    for (i=0; i<n; i++){
      d[i]=-g[i]+beta*d[i]; }
    u*=beta;
    mvsub(A,x,val,n);
    vn=0.0; vb=0.0;
    for (i=0; i<n; i++){
      vn+=val[i]*val[i]; vb+=val[i]*b[i]; }
    vn=sqrt(vn);
    if (1.0-vb/(vn*bn)<tol[0] && u/u0<=tol[0]) ic=0;
    if (k>itr[0]) ic=0;
  }
  mvsub(A,x,val,n);
  for (i=0; i<n; i++){
    val[i]=val[i]-b[i]; }
  err=0.0;
  for (i=0; i<n; i++){
    err+=val[i]*val[i]; }
  tol[1]=sqrt(err);
  return;
}
/* function */
mvsub(A,b,c,n)
double A[N][N], b[N], c[N];
int n;
{
  int i, j;
  double s;
  for (i=0; i<n; i++){
    s=0.0;
    for (j=0; j<n; j++){
      s+=A[i][j]*b[j];
    }
    c[i]=s;
  }
}