#include <stdio.h>
#include <math.h>
/* minsnln.c : 1999 5 3 by Oguni */
#define N 2
main(){
  double x[N]={-2., 5.}, A[N][N], an[N];
  double tol=0.0001, val;
  int n=2, itr[2]={10,1};
  minsnln(A,an,x,tol,n,itr);
  printf("Minimization of simultaneous non-linear equations.\n");
  printf("%f  %f \n", x[0], x[1]);
  printf("%d\n", itr[1]);
  val=100.0*(x[1]-x[0]*x[0])*(x[1]-x[0]*x[0])+(1.-x[0])*(1.-x[0]);
  printf("Substtution check %f\n", val);
}
/* function */
nabla(A,x,an,n)
double A[N][N], x[N], an[N];
int n;
{
  an[0]=-400.0*x[0]*(x[1]-x[0]*x[0])+2.0*(x[0]-1.);
  an[1]=200.0*(x[1]-x[0]*x[0]);
  A[0][0]=1200.0*x[0]*x[0]-400.0*x[1]+2.0;
  A[1][0]=-400.0*x[0];
  A[0][1]=-400.0*x[0];
  A[1][1]=200.0;
  return;
}
/* minsnln */
minsnln(A,an,x,tol,n,itr)
double A[N][N], an[N], x[N], tol;
int n, itr[2];
{
  int i,j,k;
  double tol1=0.000000001, r, s, w[N];
  for (i=1; i<=itr[0]; i++) {
    nabla(A,x,an,n);
    an[0]=-an[0];  an[1]=-an[1];
    gauss(A, an, w, tol1, n);
    for (j=0; j<n; j++){
      x[j]+=w[j];
    }
    s=0.0;
    for (k=0; k<n; k++){
      s+=w[k]*w[k]; }
    if (sqrt(fabs(s))<tol){ break; }
  }
  itr[1]=i;
  return;
}
/* function */
gauss(A, b, x, tol, n)
double A[N][N], b[N], x[N], tol;
int n;
{
  int i, j, k;
  double akk, aik, s;
  for (k=0; k<n-1; k++){
    if (fabs(A[k][k])<tol){
      printf("The pivot is too small at step # %d", k); return; }
    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]; }
      b[i]+=aik*b[k];
    }
    for (j=k+1; j<n; j++){
      A[k][j]*=akk; }
    b[k]*=akk;
  }
  x[n-1]=b[n-1]/A[n-1][n-1];
  for (k=n-2; k>=0; k--){
    s=0.0;
    for (j=k+1; j<n; j++){ s+=A[k][j]*x[j]; }
    x[k]=b[k]-s;
  }
}