#include <allinc.h>
const int S=2;
void display(apmatrix <double> &);
void main()
{
	//double list[S*S]={0,2,1,2,3,4,3,1,2,3,4,3,2,4,2,1};//{1,3,0,1,4,2,1,2,0};
   //double list[S*S]={1,3,0,1,4,2,1,2,0};
   double list[S*S]={1,2,2,1};
	apmatrix <double> a(S*S,S*S+1,0);

	for(int A=1; A<=S; A++)// <-regular A, not [A]
	{
		for(int B=1; B<=S; B++)
		{
			if(B!=1)
				a[B-1+(A-1)*S-1][B+(A-1)*S-1]=1;
			if(B!=S)
				a[B+1+(A-1)*S-1][B+(A-1)*S-1]=1;

			if(A!=1)
			{
				a[B+(A-2)*S-1][B+(A-1)*S-1]=1;
				if(B!=1)
					a[B-1+(A-2)*S-1][B+(A-1)*S-1]=1;
				if(B!=S)
					a[B+1+(A-2)*S-1][B+(A-1)*S-1]=1;
			}

			if(A!=S)
			{
				a[B+(A)*S-1][B+(A-1)*S-1]=1;
				if(B!=1)
					a[B-1+(A)*S-1][B+(A-1)*S-1]=1;
				if(B!=S)
					a[B+1+(A)*S-1][B+(A-1)*S-1]=1;
			}
		}
	}
	for(int A=1; A<=S*S; A++)
	{
		a[A-1][S*S+1-1]=list[A-1];
	}

   ofstream ouf;
   ouf.open("zmatrixOUF.txt");//writes text for Mathematica to use
   ouf<<"RowReduce[{";
   for(int i=0; i<a.numrows(); i++)
   {
   	ouf<<"{";
   	for(int j=0; j<a.numcols(); j++)
      {
      	if(j<a.numcols()-1)
	      	ouf<<char(0x61+i)<<char(0x61+j);//a[i][j];
         else
         	ouf<<a[i][j];
         if(j<a.numcols()-1)
         	ouf<<",";
      }
      ouf<<"}";
      if(i<a.numrows()-1)
      	ouf<<",";
   }
   ouf<<"}]";
   ouf.close();




   display(a);

   //where RREFing begins
   for(int i=0; i<a.numrows()-1; i++)
   {
      //make sure it's not zero, swap if it is
   	if(a[i][i]==0)
      {
      	for(int j=i+1; j<a.numrows()-1; j++)
         	if(a[j][i]!=0)
            {
            	for(int k=0; k<a.numcols(); k++)
               {
               	double temp=a[j][k];
                  a[j][k]=a[i][k];
                  a[i][k]=temp;
               }
               j=a.numrows()-1;
            }
      }
      //pivot element is now non-zero
      for(int j=i+1; j<a.numrows(); j++)
      {
      	//check to see if target row is nonzero, else skip
      	if(a[j][i]!=0)
         {
         	//check if target row needs factored up a bit
            if(a[j][i]-a[i][i]*int(a[j][i]/a[i][i])!=0)
            	for(int k=0; k<a.numcols(); k++)
               	a[j][k]*=a[i][i];

            //now do math to cancel out that row.
            double factor=a[j][i]/a[i][i];
            for(int k=0; k<a.numcols(); k++)
            	a[j][k]-=a[i][k]*factor;
         }
      }
      clrscr(); display(a);
   }

   //a[8][8]=1;
   //a[8][9]=1;
   //clrscr(); display(a);

   //reverse direction
   for(int i=a.numrows()-1; i>=1; i--)
   {
   	if(a[i][i]!=0)
      {
         a[i][a.numcols()-1]/=a[i][i];
      	a[i][i]/=a[i][i];
      	for(int j=i-1; j>=0; j--)
         	if(a[j][i]!=0)
            {
            	//factor up target if necessary
            	if(a[j][i]-a[i][i]*int(a[j][i]/a[i][i])!=0)
               	for(int k=0; k<a.numcols(); k++)
                  	a[j][k]*=a[i][i];

               //cancel out
               double factor=a[j][i]/a[i][i];
               for(int k=0; k<a.numcols(); k++)
               	a[j][k]-=a[i][k]*factor;
            }
		   clrscr(); display(a);
      }
   }
   cout<<"\ndone."; getch();
}

void display(apmatrix <double> &a)
{
   for(int i=0; i<a.numrows(); i++)
   {
   	for(int j=0; j<a.numcols(); j++)
      	cout<<setw(5)<<setprecision(2)<<a[i][j]<<" ";
      cout<<endl;
   }
   getch();
}