#include <stdio.h>
#include <math.h>
/* invert.c : 1999 3 8 by Oguni */
#define N 3
main(){
  int n=3, itr[2]={20, 0};
  int i;
  double A[N][N]={3.,2.,1.,2.,3.,2.,1.,2.,3.};
  double x[N+1]={1.,0.,0.,0.}, tol=0.00001;
  invert(A,x,tol,n,itr);
  printf("%d %f\n", itr[1], x[n]);
}
/* function invert */
invert(A,x,tol,n,itr)
double A[N][N], x[N+1], tol;
int n, itr[2];
{
  double L[N][N], U[N][N], y[N], s[N], c, d, xmin;
  int m, i;
  lu(A,L,U,tol,n);
  for (m=1; m<=itr[0]; m++){
    for (i=0; i<n; i++){ s[i]=x[i]; }
    subs(L,U,s,y,tol,n);
    c=0; d=0;
    for (i=0; i<n; i++){
      c+=y[i]*x[i];  d+=y[i]*y[i]; }
    c=c/d;
    if (fabs(xmin-c)<tol) return;
    xmin=c; x[n]=xmin; itr[1]=m;
    for (i=0; i<n; i++){ x[i]=y[i]/sqrt(d); }
  }
  itr[1]=m;
  return;
}
/* function lu */
lu(A,L,U,tol,n)
double A[N][N], L[N][N], U[N][N], tol;
int n;
{
  double akk, aik;
  int i, j, k;
  for (k=0; k<n-1; k++){
    if (fabs(A[k][k])<tol) {
      printf("The pivot is too small at step # %d", k); }
    akk=1.0/A[k][k];
    for (i=k+1; i<n; i++){
      aik=-A[i][k]*akk;
      for (j=k+1; j<n; j++){ A[i][j]+=aik*A[k][j]; }
    }
    for (j=k+1; j<n; j++){ A[k][j]*=akk; }
  }
  for (i=0; i<n-1; i++){
    for (j=0; j<=i; j++){ L[i][j]=A[i][j]; }
    for (j=i+1; j<n; j++){ U[i][j]=A[i][j]; }
    U[i][i]=1.0;
  }
  U[n-1][n-1]=1.0;
  for (j=0; j<n; j++){ L[n-1][j]=A[n-1][j]; }
}
/* function subs */
subs(L,U,b,x,tol,n)
double L[N][N], U[N][N], b[N], x[N], tol;
int n;
{
  double y[N], s;
  int i, j, k;
  y[0]=b[0]/L[0][0];
  for (i=1; i<=n-1; i++){
    if (fabs(L[i][i])<tol) {
      printf("The pivot is too small at step #%d", i); }
    s=0.0;
    for (j=0; j<=i-1; j++){
      s+=L[i][j]*y[j]; }
    y[i]=(b[i]-s)/L[i][i];
  }
  x[n-1]=y[n-1]/U[n-1][n-1];
  for (k=n-2; k>=0; k--){
    s=0.0;
    for (j=k+1; j<n; j++){ s+=U[k][j]*x[j]; }
    x[k]=(y[k]-s)/U[k][k];
  }
}