// BRYAN WEIR
// Problem Set 5
// PART 1
// 7 April 97
	
#include <iostream.h>
#include <fstream.h>
#iinclude <stdio.h>
#include <math.h>
#include <stdlib.h>

#define NULL 0

class Matrix;
class Vector;
class Solver;

class Matrix
{
  int numrows, numccols;
  double *data;
  
public: 
  Matrix()
  {
    numrows = 0;
    numcols = 0;
    data = NULL;
  }
Matrix(int r, int c)
  {
    numrows = r;
    numcols = c;
    data = new double[r*c];
  }

  void set(int r, int c, double d)
  {
    if (data == NULL)
      exit;
    if (r<1 || c<1 || r>numrows || c>numcols)
      {
	cout << "Error!!! Bad dimensions" << endl;
	exit;
      }
    data[numcols*(r-1)+(c-1)] = d; 
  }

  double get(int r, int c)
  {
    if (data == NULL)
      exit;
     if (r<1 || c<1 || r>numrows || c>numcols)
      {
	cout << "Error!!! Bad dimensions" << endl;
	exit;
      }
   return data[numcols*(r-1)+(c-1)]; 
  }

  void set_cols_rows(int r, int c)
  {
    if (numrows || numcols)
      exit;
    numrows = r;
    numcols = c;
    data = new double[r*c];
  }
  
  int getrows()
  {
    if (!numrows)
      exit;
    return numrows;
  }

  int getcols()
  {
    if (!getcols)
      exit;
    return numcols;
  }

  Matrix *copy()
  {
    int i,j;

    Matrix *A = new Matrix(this->numrows, this->numcols);
    for (i=0; i<=numrows; i++)
      for (j=0; j<=numcols; j++)
	A->set(i,j, this->get(i,j));
    return A;
  }

  void swap(int row1, int row2)
  {
    double temp;
    int i;

    if (row1<1 || row2<1 || row1>numrows || row2>numrows)
      exit;
    for (i=1; i<numcols; i++)
      {
	temp = this->get(row1, i);
	this->set(row1, i, this->get(row2, i));
	this->set(row2, i, temp);
      }
  }
  
  void printmx()
  {
    int i, j;
    
    for(i=1; i<=numrows; i++)
      {
	for (j=1; j<numcols; j++)
	  cout << this->get(i, j) << " ";
	cout << endl;
      }
  }
}; 

 
class Vector : public Matrix
{
public:
  Vector()
  {
    numcols = 1;
    numrows = r;
    data = null;
  }

  Vector(int r)
  {
    numcols = 1;
    numrows = r;
    data = new double[r];
  }

  double get(int r)
  {
    if (data == NULL)
      exit;
    if (r,1 || r>numrows)
      {
	cout << "Error!!  Bad Dimensions." << endl;
	exit;
      }
   return data[numcols*(r-1)];  
  }
  void set(int r, double d)
  {
    if (data == NULL)
      exit;
    if (r<1 || r>numrows)
      {
	cout << "Error!!  Bad Dimensions." << endl;
	exit;
      }
    data[numcols*(r-1)] = d;
  }

};

class Solver
{
  Matrix *M;
  Matrix *N;
  Vector *v;

 public:
  Solver()
    {
      M = NULL;
      N = NULL;
      v = NULL;
    }
  Solver(Matrix *that_matrix)
    {
      int i;
      
      N = that_matrix->copy();
      v = new Vector(that_matrix->getrows());
      M = NULL;

      for(i=1; i<=N->getrows(); i++)
	b->set(i,0);
    }

  Matrix *getN()
    {
      return N;
    }
  
  Vector *getv()
    {
      return v;
    }
  
  void setv(Vector *random)
    {
      int i;
      v = new Vector(random->getrows());
      for(i=1; i<random->getrows(); i++)
	v->set(i, random->get(i));
    }
  void setup()
    {
      if (N==NULL)
	exit(-1);
      if (N->getrows() != N->getcols())
	exit(-1);
      
      int i,j;
      int size = N->getrows();
      M = new Matrix(size, 2*size+2);

      for(i=1; i<=size; i++)
	for(j=1; j<=size*2+2; j++)
	  {
	    if(j<=size)
	      M->set(i, j, N->get(i,j));
	    if(j==(size+1))
	      M->set(i,j, v->get(i));
	    if((j>(size+1)) && (j<=(2*size+1)))
	      if(j-i-N->getrows()==1)
		M->set(i,j,1);
	      else
		M->set(i,j,0);
	    if(j==(2*size+2))
	      M->set(i,j,i);
	  }
    }
  
M->print();
}


void solve(double *det, Matrix *inv)
{
  int size=A->getrows();
  double factor;
  int a=1;
  int i, j, k, l, m, n, largest;
  *det = 1;
  
  for(i=1; i<=size; i++)
    {
      largest = i;
      
      for(j=1; j<=size; j++)
	if(absolut(M->get(j,i)) > absolut(M->get(largest, i)))
	  largest = j;
    }
  
  if(largest!=i)
    {
      M->swap(largest, i);
      a = -a;
    }
  
  for(k=i+1; k<=size; k++)
    {
      factor = M->get(k,i)/M->get(i,i);
      for(l=i; l<=2*size+1; l++)
	{
	  M->set(k,l,M->get(k,l) - factor*(M->get(i,l)));
	}
    }
  
  M->print();
  
  for (n=1; n<=size; n++)
    (*det) = (*det)*(M->get(n,n));    
  (*det) = (*det)*a;
  
  if(*det==0)
    return;
  
  for(i=size; i>=1; i--)
    {
      for(k=i-1; k>=1; k--)
	{
	  factor = M->get(k,i)/M->get(i,i);
	  for(m=i; m<=(2*size+1); m++)
	    {
	      M->set(k, m, M->get(k,m) - factor*(M->get(i, m)));
	    }
	}
    }
  M->print();
  cout << endl;
  
  for(i=1; i<=size; i++)
    {
      for(m=(2*size+1); m>=i; m--)
	M->set(i, m, M->get(i,m)/M->get(i,i));
    }

  M->print();
  cout << endl;

  (*inv) = *(new Matrix(size, size));
  for(i=1; i<=size; i++)
    for(j=1; j<=size; j++)
      (*inv).set(i, j, M->get(i, j+size+1));

  inv->print();

};

Matrix *Matrix::inverse()
{
  Solver S(this);
  S.setup();
  double det;
  Matrix *inv = new Matrix(this->getrows(), this->getcols());
  S.solve(&det, inv);
  return inv;
}

double Matrix::determ()
{
  Solver S(this);
  S.setup();
  double det;
  Matrix inv(this->getrows(), this->getcols());
  S.solve(&det, &inv);
  return det;
}


double absolut(double n)
{
  if (n<0)
    return (-n);
  else
    return n;
}




Matrix operator+ (Matrix M1 , Matrix M2)
{
  int i, j;

  if ((M1->checkrows() != M2->checkrows ()) ||
      (M1->checkcols() != M2->checkcols()))
    {
      cout << "BZZZZZZZZ!!!  Error, unmatched crap, no addition goin'on here!"
	   << endl;
      exit(-1);
    }
  Matrix *temp = new Matrix(M1.checkrows(), M1.checkcols());

  for (int i=0; i<M1.checkrows(); i++)
    {
      for (int j=0; j<M1.checkcols(); j++)
	{
	  temp->set(i, j, M1.get(i, j) + M2.get(i, j));
	}
    }
  return (*temp);
};


Matrix operator- (Matrix M1 , Matrix M2)
{
 int i, j;

  if ((M1->checkrows() != M2->checkrows ()) ||
      (M1->checkcols() != M2->checkcols()))
    {
      cout <<"BZZZZZZZZ!!  Error, unmatched crap, no subtraction goin'on here!"
	   << endl;
      exit(-1);
    }
  Matrix *temp = new Matrix(M1.checkrows(), M1.checkcols());
  
  for(int i=0; i<M1.checkrows(); i++)
    {
      for(int j=0; j<M1.checkcols; j++)
	{
	  temp->set(i,j, M1.get(i,j) - M2.get(i,j));
	}
    }
  return (*temp);
};


Matrix operator* (double multiplier, Matrix M2)
{
  int i,j;

  Matrix *temp = new Matrix(M1.checkrows(), M1.checkcols());

  for(i=0; i<m1.checkrows(); i++)
    {
      for (j=0; j<M1.checkcols(); j++)
	{
	  temp->set(i,j,M1.get(i,j,)*factor);
	}
    }
  return (*temp);
};


Matrix operator- (Matrix M1)
{
  return (-1*M1);
};



void input_file(char *file, int freedom)
{
  ifstream crap ("ps5_q1.in");

  crap >> freedom >> L.load1 >> L.load2 >> L.load3 >> L.load4 >> S.stiffy11
       >> S.stiffy12 >> S.stiffy13 >> S.stiffy14 >> S.stiffy21 >> S.stiffy22 
       >> S.stiffy23 >> S.stiffy24 >> S.stiffy31 >> S.stiffy32 >> S.stiffy33 
       >> S.stiffy34 >> S.stiffy41 >> S.stiffy42 >> S.stiffy43 >> S.stiffy44;
  
  crap.close();
};

void main()
{
  int freedom;

};