int _stklen=100000;

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

#define MAXI 200
#define MAXJ 200

float m[MAXI][MAXJ];
float m2[MAXI][MAXJ];
float b[MAXI], x[MAXI];
float b2[MAXI];
float p;
int di, dj;
int n;

void inits();
void readmats();
void checkmats();
void printmats();
void printx();
void invertmats();


void inits() {
  int i,j;
  di=dj=0;
  x[0] = 0; x[1] = 0;
  return;
}

void readmats() {
   int i, j;
   float t;
   di = dj = 4;
   m[0][0] = 7.0; m[0][1] = 2.0; m[0][2] = 1.0;
   m[1][0] = 2.0; m[1][1] = 4.0; m[1][2] = 1.0;
   m[2][0] = 1.0; m[2][1] = -1.0; m[2][2] = 60.0;
   b[0] = 2.0; b[1] = 5.0; b[2] = 0.3;
  randomize();
/*  n = 5; */
  di = dj = n;
  for (i=0; i<n; i++)
     for (j=0; j<n; j++) {
        t = (float) random(10);
        m[i][j] = t;
        m2[i][j] = t;
     }
  for (i=0; i<n; i++) {
      t = (float) random(10);
      b[i] = t;
      b2[i] = t;
   }
   return;
}

/* void checkmats2()
{
   int i,j;
   float *p1, t,d1,d2, e;

   for (i=0; i<di; i++) {
      t = 0.0;
      for (j=0; j<dj; j++) {
         p1 = &m2[i][j];
         d1 = *p1;
         p1 = x[j];
         d2 = *p1;
         d1 *= d2;
      }
      e= b2[i];
      t -= e;
      d2 = b2[i];
      printf("%f %f\n",(double) t,(double) d2);
   }
   return;
} */

void checkmats()
{
   int i,j;
   float *p, d1, d2, d3;
   for (i=0; i<di; i++) {
      d3 = 0.0;
      for (j=0; j<dj; j++) {
         p = &m2[i][j];
         d1 = *p;
         p = &x[j];
         d2 = *p;
         d3 += d1*d2;
         /* printf(" %f",(double) d1); */
      }
      d3 -= b2[i];
      printf(" %f %f\n",(double) d3, (double)  b2[i]);
   }
   return;
}

void printmats(s) 
char s[];
{
   int i,j;
   if (s!=NULL && s[0]!='\0') printf("%s\n",s);
   for (i=0; i<di; i++) {
      for (j=0; j<dj; j++)
         printf(" %f",(double) m[i][j]);
      printf("   %00.2f\n",(double)b[i]);
   }
   return;
}

void printx(s) 
char s[];
{
   int i,j;
   if (s!=NULL && s[0]!='\0') printf("%s\n",s);
   for (i=0; i<di; i++) {
      printf(" %f\n",(double) x[i]);
   }
   return;
}


void invertmats() {
   int n,i,j;
   float *y, p1, d1;

   for (n=0; n<di-1; n++) {      /* With all diagonal elements. */
      for (i=di-1; i>n; i--) {   /* Make zeros starting at the lowest row. */
         y = &m[n][n];           /* The pointer stuff is because of a definate */
         p1 = *y;                /* compiler bug, it actually hangs my whole */
         y = &m[i][n];           /* when doing a 2d array assignent to a var */
         d1 = *y;                /* bad, bad. This works, though */
         d1 /= p1;               /* This is the factor to multiply the eq with */
         for (j=n; j<dj; j++) {  /* Now multiply all relevant row elements */
            y = &m[n][j];
            p1=*y;
            m[i][j] -= d1*p1;    /* subtract the higher equation */
         }
         /*   m[i][j] = m[i][j] - p * m[n][j]; */

         b[i] = b[i] - d1 * b[n]; /* and don't forget the right hand side */
      }
   if (0) printmats("  ");              /* do a matrix print, for small mats */
   }
                                 /* Now the sys is upper tri */
   for (n=di-1; n>0; n--) {      /* sweep clear the upper tri */
      for (i=n-1; i>=0; i--) {   /* vertical stripes of zeros */
         y = &m[n][n];
         p1 = *y;
         y = &m[i][n];
         d1 = *y;
         d1 /= p1;
         m[i][n] -= d1*p1;       /* make use of the existing zeros */
         b[i] = b[i] - d1 * b[n]; /* and the right hand side of course */
      }
   if (0) printmats("  ");
   }
                                 /* Now we have a diagonal matrix */
   for (i=0; i<di; i++) {
      y = &m[i][i];
      p1 = *y;
      x[i] = b[i] / p1;       /* write solution vector elements */
   }
   return;
}

int main(argc,argv)
int argc;
char *argv[];
{
   inits();
   if (argc == 2)
      if (sscanf(argv[1],"%d",&n) != 1)
          return(-1);
   readmats(n);
   printmats("read:");
   printf("Starting invert... \n");
   invertmats();
/*   printmats("out:"); */
   checkmats();
   printx("x:");
   return 0;
}