//Matrixd.cpp
#include <fstream.h>
#include <math.h>
//#include "Readfile.cpp"

Matrixd::Matrixd()
{
	rows = 3;
	cols = 3;
	int r;
	int c;
      	array = new double [rows*cols];
	for(r = 0; r < rows; r++)
	{
		for(c = 0; c < cols; c++)
		{
			array[r*cols + c] = 1.0;
		}
	}
}

Matrixd::~Matrixd()
{
	delete [] array;
}

Matrixd::Matrixd(int rowsin, int colsin, double * Matrixdin)
{
	rows = rowsin;
	cols = colsin;
	int r = 0;
	int c = 0;
	//array = Matrixdin;
        array = new double [rows*cols];
	for(r = 0; r < rows; r++)
	{
		for(c = 0; c < cols; c++)
		{
			array[r*cols + c] = Matrixdin[r*cols + c];
                        //Matrixdin[r*cols + c];
		}
	}
}

Matrixd::Matrixd(int rowsin, int colsin)		//if you know the size but not the contents yet!
{
	rows = rowsin;
	cols = colsin;
	int r = 0;
	int c = 0;
	//array = Matrixdin;
        	array = new double [rows*cols];
	for(r = 0; r < rows; r++)
	{
		for(c = 0; c < cols; c++)
		{
			array[r*cols + c] = 1.0;
                        	}
	}
}

Matrixd::Matrixd(int rowsin, double val)
{
        //this one makes a column matrix with val as each entry
	rows = rowsin;
	//cols = colsin;
	int r = 0;
	//int c = 0;
	//array = Matrixdin;
        array = new double [rows*cols];
	for(r = 0; r < rows; r++)
	{
                array[r] = val;
	}
}


Matrixd::Matrixd(int rowsin, int colsin, double val)
{
        //This constructor just fille matrix with all the same value
	rows = rowsin;
	cols = colsin;
	int r = 0;
	int c = 0;
	//array = Matrixdin;
        	array = new double [rows*cols];
	for(r = 0; r < rows; r++)
	{
		for(c = 0; c < cols; c++)
		{
			array[r*cols + c] = val;
                        	}
	}
}

Matrixd::Matrixd(int rowsin, int colsin, ifstream & in)
{
        //long n;
        rows = rowsin;
        cols = colsin;
	array = new double [rows*cols];
        double temp;
	//n = FindString(in, "Matrixd:");
	for(int r = 0; r < rows; r++)
        {
                for(int c = 0; c < cols; c++)
                {
                        temp = get_double(in);
                        array[r*cols + c] = temp;
                }
        }
}



Matrixd::Matrixd(Matrixd& m)
{
	rows = m.rows;
	cols = m.cols;
	array = new double [rows*cols];
	for(int r = 0; r < rows; r++)
	{
		for(int c = 0; c < cols; c++)
		{
			array[r*cols + c] = m.val(r,c);
                        //m.array[r*cols + c];
		}
	}
}

Matrixd Matrixd::operator+(Matrixd & Matrixdin)
{
	if(rows==Matrixdin.rows && cols==Matrixdin.cols)
	{
		double *arrayptr = new double [rows*cols];
		for(int r = 0; r < rows; r++)
		{
			for(int c = 0; c < cols; c++)
			{
				arrayptr[r*cols + c] = array[r*cols + c] +
									Matrixdin.array[r*cols + c];
			}
		}
		Matrixd m(rows, cols, arrayptr);
		delete [] arrayptr;
		return m;
	}
	else
	{
		Matrixd m;
		return m;
	}
}

Matrixd Matrixd::operator-(Matrixd & Matrixdin)
{
	if(rows==Matrixdin.rows && cols==Matrixdin.cols)
	{
		double *arrayptr = new double [rows*cols];
		for(int r = 0; r < rows; r++)
		{
			for(int c = 0; c < cols; c++)
			{
				arrayptr[r*cols + c] = array[r*cols + c] -
										Matrixdin.array[r*cols + c];
			}
		}
		Matrixd m(rows, cols, arrayptr);
		delete [] arrayptr;
		return m;
	}
	else
	{
		Matrixd m;
		return m;
	}
}

Matrixd Matrixd::operator*(Matrixd & Matrixdin)
{
	//Matrixd Multiplication
	if(cols==Matrixdin.rows)
	{
		int rows2 = Matrixdin.rows;
		int cols2 = Matrixdin.cols;
		double *arrayptr = new double [rows*cols2];
		int r = 0;
		int c = 0;
		int k = 0;
		double sum = 0;
		for(r = 0; r<rows; r++)
		{
			for(c = 0; c<cols2; c++)
			{
				sum = 0.0;
				//arrayptr[r*cols2 + c] = 0;
				for(k = 0; k<cols; k++)
				{
					sum += array[r*cols + k] *
								Matrixdin.array[k*cols2+c];
				}
				arrayptr[r*cols2+c] = sum;
			}
		}
		Matrixd m(rows, cols2, arrayptr);
		delete [] arrayptr;
		return m;
	}
	else
	{
		Matrixd m;
		return m;
	}
}


Matrixd & Matrixd::operator=(Matrixd & Matrixdin)
{
	rows = Matrixdin.rows;
	cols = Matrixdin.cols;
	//array = Matrixdin.array;
        double temp;

        array = new double [rows*cols];
	for(int r = 0; r < rows; r++)
	{
		for(int c = 0; c < cols; c++)
		{
                        temp = Matrixdin.val(r,c);
			array[r*cols + c] = temp;
                        //Matrixdin.array[r*cols + c];
		}
	}
	return *this;
}



void Matrixd::printMatrixd(ofstream & out)
{
	out.precision(4);
	out.setf(ios::scientific);
	int r = 0;
	int c = 0;
	for(r=0; r<rows; r++)
	{
		//out.width(12);
		//out<<energy.en[j];
		for(c=0; c<cols; c++)
		{
			out.width(12);
			out<< A(r, c);
			//array[r*cols + c];
		}
		out<<endl;
	}
}


double Matrixd::A(int r, int c)
{
	return array[r*cols + c];
}


double Matrixd::A(int r)
{
	return array[r*cols];		//for a column vector
}

double Matrixd::val(int r, int c)
{
        double temp = array[r*cols + c];
        return temp;
}

double Matrixd::val(int r)
{
        double temp = array[r*cols];
	return temp;		//for a column vector
}

void Matrixd::set(int r, int c, double it)
{
	array[r*cols + c] = it;
}

void Matrixd::set(int r, double it)
{
        array[r*cols] = it;         //for column vector
        //return 0;
}


double Matrixd::addrow(int r)
{
	double temp = 0.0;
	for(int c = 0; c < cols; c++)
	{
		temp += val(r,c);
	}
	return temp;
}

double Matrixd::addcol(int c)
{
	double temp = 0.0;
	for(int r = 0; r < rows; r++)
	{
		temp += val(r,c);
	}
	return temp;
}



void Matrixd::transpose()
{
        //transpose the Matrixd
        //no need to redo the array... same size as before
        int old_rows = rows;
        int old_cols = cols;
        rows = old_cols;
        cols = old_rows;
        Matrixd temp;
        temp = *this;
        for(int r = 0; r < rows; r++)
        {
                for(int c = 0; c < cols; c++)
                {
                     array[r*cols+c] = temp.array[c*rows+r];
                }
        }
}

void Matrixd::augment(Matrixd & Matrixdin)
{
        //Just add new columns onto end of *this
        int oldcols = cols;
        int colsin = Matrixdin.cols;
        int rowsin = Matrixdin.rows;
        Matrixd tempin;
        tempin = Matrixdin;
        if(rowsin < rows)
        {
		double *arrayptr = new double [rows*colsin];
                int diff = rows - rowsin;
                int el_diff = diff * colsin;
                int element;
                for(int r2 = 0; r2 < rows; r2++)
                {
                        for(int c2 = 0; c2 < colsin; c2++)
                        {
                                element = r2*colsin+c2;
                                if(element < el_diff)
                                        arrayptr[element] = 0.0;
                                else
                                {
                                        arrayptr[element] =
                                                Matrixdin.array[element - el_diff];
                                }
                        }
                }
                Matrixd addedrows(rows, colsin, arrayptr);
                tempin = addedrows;
        }
        Matrixd temp;
        temp = *this;
        delete [] array;
        cols += colsin;
	array = new double [rows*cols];
        double tempdouble;
        int element;
        //Fill in whole "new " array
        for(int r = 0; r < rows; r++)
        {
                for(int c = 0; c < cols; c++)
                {
                        element = r*cols + c;
                        if(c<oldcols)
                        {
                                tempdouble = temp.array[r*oldcols+c];
                                array[element] = tempdouble;
                        }
                        else
                        {
                                //tempdouble = tempin.array[r*colsin+c-colsin];
                                tempdouble = tempin.array[element - (r+1)*oldcols];
                                array[element] = tempdouble;
                        }
                }
        }
}

Matrixd Matrixd::extract(int n)
{
        double temp[10000];  //largest inverse is 100x100
        for(int r = 0; r < rows; r++)
        {
                for(int c = cols - n; c < cols; c++)
                {
                        temp[r*n + c -n] = this->val(r,c);
                }
        }
        Matrixd B(rows, n, temp);
        return B;
}



Matrixd Matrixd::ident(int n)
{
        delete [] array;
        rows = n;
        cols = n;
	array = new double [rows*cols];
        //double * im = new double[n*n];            //identity Matrixd vector
        for(int i = 0; i < n; i++)
        {
                for(int j = 0; j < n; j++)
                {
                        if(i==j)
                                array[i*n+j]= 1.0;
                        else
                                array[i*n+j]= 0.0;
                }
        }
}

Matrixd Matrixd::GJ(ofstream & out)
{
     //Gauss-Jordan Elim.
    //The coefficient Matrixd is *this
    //This version is exclusively for Matrixd inversion...
    //i.e. no argument sent.

	int n = rows;				//number of equations
        Matrixd I;
        I.ident(n);
        //Now augment A with the identity Matrixd
        Matrixd A;
        A = *this;
        A.augment(I);
        int ncols = A.cols;     //this is # of cols in new augmented Matrixd A
        //return A;

	double dum;
	double factor;
	double sum =0;
	int i = 0;
	int j = 0;
	int k = 0;
	int ii = 0;

	int * o = new int[n];		//Order Vector
	double * s = new double[n];	//Scale Vector
	//double * x1 = new double[n];	//Solution Array
	//Matrixd x(n, 1, x1);		//Solution Vector
        //double check;
	//Matrixd A;
	//Matrixd C;

	//A = *this;
	//C = cin;

	 //Order the equations
	for(i = 0; i < n; i++)
	{
		o[i] = i;
		s[i] = fabs(A.A(i,0));
		for(j = 1; j < ncols; j++)
		{
			if (fabs(A.A(i,j)) > s[i])
				s[i] = fabs(A.A(i,j));
		}
	}
        //Now s[] is largest magnitude in each row (identified as o[]

	//Elimination
	for(k = 0; k < n; k++)
	{
		//partial pivot
                //i.e. switch equations so that largest pivot is used
		int pivot = k;
		double big = fabs( A.A(o[k],k) / s[o[k]] );
                //A.A(o[k],k) is the normal pivot...
                //Now search forward through remaining rows for the best pivot
		for(ii = k+1; ii < n; ii++)
		{
			dum = fabs( A.A(o[ii],k) / s[o[ii]] );
			if(dum > big)
			{
				big = dum;
				pivot = ii;
			}
		}
		dum = o[pivot];         //get ready to swap rows
		o[pivot] = o[k];        //move the best pivot row to current row
        	o[k] = dum;             //move former pivot to former position
                                        //of best pivot row
                //Now that we've settled on best pivot row...
                //normalize it
                double temp2 = A.val(o[k], k);
                for(int iii=0; iii < ncols; iii++)
                {
                        double temp1 = A.val(o[k], iii)/ temp2;
                        //double temp3 = temp1 / temp2;
                        A.set(o[k], iii, temp1);
                }
                //out << endl;
                //A.printMatrixd(out);

		//continue elimination for all rows above and below
		for(i = 0; i < n; i++)
		{
                        if(i != k)
                        {
                                //factor = A.A(o[i],k) / A.A(o[k],k);
                                factor = A.A(o[i],k);
        			for(j = 0; j < ncols; j++)
        			{
        				A.array[o[i]*ncols+j] -= factor*A.A(o[k],j);
	        		}
		        	//check = A.array[o[i]*ncols+j];
                                //C.array[o[i]] -= factor * C.array[o[k]];
                        }
		}
                //out << endl;
                //A.printMatrixd(out);
                //now extract the Matrixd inverse from augmented Matrixd
	}
        Matrixd B;
	B=A;
        B = A.extract(n);
        return B;
}

Matrixd Matrixd::GJ(Matrixd & Matrixdin)
{
     //General Gauss-Jordan Elim.
    //The coefficient Matrixd is *this

	int n = rows;				//number of equations
        //Now augment A with the input Matrixd
        Matrixd A;
        A = *this;
        A.augment(Matrixdin);
        int ncols = A.cols;     //this is # of cols in new augmented Matrixd A
        //return A;

	double dum;
	double factor;
	double sum =0;
	int i = 0;
	int j = 0;
	int k = 0;
	int ii = 0;

	int * o = new int[n];		//Order Vector
	double * s = new double[n];	//Scale Vector

        //Order the equations
	for(i = 0; i < n; i++)
	{
		o[i] = i;
		s[i] = fabs(A.A(i,0));
		for(j = 1; j < ncols; j++)
		{
			if (fabs(A.A(i,j)) > s[i])
				s[i] = fabs(A.A(i,j));
		}
	}
        //Now s[] is largest magnitude in each row (identified as o[]

	//Elimination
	for(k = 0; k < n; k++)
	{
		//partial pivot
                //i.e. switch equations so that largest pivot is used
		int pivot = k;
		double big = fabs( A.A(o[k],k) / s[o[k]] );
                //A.A(o[k],k) is the normal pivot...
                //Now search forward through remaining rows for the best pivot
		for(ii = k+1; ii < n; ii++)
		{
			dum = fabs( A.A(o[ii],k) / s[o[ii]] );
			if(dum > big)
			{
				big = dum;
				pivot = ii;
			}
		}
		dum = o[pivot];         //get ready to swap rows
		o[pivot] = o[k];        //move the best pivot row to current row
        	o[k] = dum;             //move former pivot to former position
                                        //of best pivot row
                //Now that we've settled on best pivot row...
                //normalize it
                double temp2 = A.val(o[k], k);
                if( temp2 == 0 )
                        temp2 = 1e-8;
                for(int iii=0; iii < ncols; iii++)
                {
                        double temp1 = A.val(o[k], iii)/ temp2;
                        //double temp3 = temp1 / temp2;
                        A.set(o[k], iii, temp1);
                }

		//continue elimination for all rows above and below
		for(i = 0; i < n; i++)
		{
                        if(i != k)
                        {
                                //factor = A.A(o[i],k) / A.A(o[k],k);
                                factor = A.A(o[i],k);
        			for(j = 0; j < ncols; j++)
        			{
        				A.array[o[i]*ncols+j] -= factor*A.A(o[k],j);
	        		}
                        }
		}
	}
        return A;
        double temp = 1.0;
}

Matrixd Matrixd::Inverse()
{
         //Matrixd Inverse
         //send identity Matrixd to augment for Gauss-Jordan
	int n = rows;				//number of equations
        Matrixd I;
        I.ident(n);
        Matrixd A;
        A = this->GJ(I);
        Matrixd B;
        B = A.extract(n);
        return B;
}

double Matrixd::condition()
{
        Matrixd Cnorm;   //make a copy of C
        Cnorm = *this;
        Cnorm.normalize();
        double temp = Cnorm.norm();
        Matrixd D;
        D = Cnorm.Inverse();
        double temp2 = D.norm();
        double cond = temp * temp2;
        return cond;
}

double Matrixd::norm()
{
        //find uniform row norm
        double sum;
        double maxsum = 0.0;
        for(int r = 0; r < rows; r++)
        {
                sum = 0.0;
                for(int c = 0; c < cols; c++)
                {
                        sum += fabs(this->val(r,c));
                }
                if(sum > maxsum)
                        maxsum = sum;
        }
        return maxsum;
}

void Matrixd::normalize()
{
        //normalize the Matrixd so that largest element in each row is 1
        double big;
        double currvalue;
        for(int r = 0; r < rows; r++)
        {
                big = 0;
                //big is largest value for each row
                for(int c = 0; c < cols; c++)
                {
                        currvalue = fabs(this->val(r,c));
                        if(currvalue > big)
                                big = currvalue;
                }
                //Now divide all elements by largest value
                for(int c = 0; c < cols; c++)
                {
                        if(big == 0)
                                big = 1e-10;

                        this->set(r,c,this->val(r,c)/big);
                }
        }
}

Matrixd Matrixd::Gauss(Matrixd & cin)
{

    //The coefficient Matrixd is *this

	int n = rows;				//number of equations

	double dum;
	double factor;
	double sum =0;
	int i = 0;
	int j = 0;
	int k = 0;
	int ii = 0;

	int * o = new int[n];		//Order Vector
	double * s = new double[n];	//Scale Vector
	double * x1 = new double[n];	//Solution Array
	Matrixd x(n, 1, x1);			//Solution Vector
	Matrixd A;
	Matrixd C;

	A = *this;
	C = cin;

	 //Order the equations
	for(i = 0; i < n; i++)
	{
		o[i] = i;
		s[i] = fabs(A.A(i,0));
		for(j = 1; j < n; j++)
		{
			if (fabs(A.A(i,j)) > s[i])
				s[i] = fabs(A.A(i,j));
		}
	}

	//Elimination
	for(k = 0; k < n-1; k++)
	{
		//partial pivot
		int pivot = k;
		double big = fabs( A.A(o[k],k) / s[o[k]] );
		for(ii = k+1; ii < n; ii++)
		{
			dum = fabs( A.A(o[ii],k) / s[o[ii]] );
			if(dum > big)
			{
				big = dum;
				pivot = ii;
			}
		}
		dum = o[pivot];
		o[pivot] = o[k];
		o[k] = dum;

		//continue elimination
		for(i = k+1; i < n; i++)
		{
			factor = A.A(o[i],k) / A.A(o[k],k);
			for(j = k+1; j < n; j++)
			{
				A.array[o[i]*cols+j] -= factor*A.A(o[k],j);
			}
			C.array[o[i]] -= factor * C.array[o[k]];
		}
	}

	//Back Substitute
        if(A.A(o[n-1],n-1) != 0)
                x.array[n-1] = C.array[o[n-1]] / A.A(o[n-1],n-1);
        else
                x.array[n-1] = 1e7;
	for(i = n-2; i >= 0; i--)
	{
		sum = 0;
		for(j = i+1; j < n; j++)
		{
			sum += A.A(o[i],j) * x.array[j];
		}
		x.array[i] = (C.array[o[i]] - sum) / A.A(o[i],i);
	}

	delete [] x1;
	delete [] o;
	delete [] s;
	return x;
}

double Matrixd::test_iter_refine(Matrixd & bin, ofstream & out)
{
	Matrixd a;
	Matrixd b;
	a = *this;
	b = bin;
        Matrixd ainv;
        ainv = a.Inverse();
        /*
        a.printMatrixd(out);
        out << endl;
        b.printMatrixd(out);
        out << endl;
        */
        //ainv.printMatrixd(out);
        Matrixd x1;
        Matrixd x2;
        x1 = a.Inverse()*b;         //from a * x = b
        x2 = a.Gauss(b);         //from a * x = b
        /*
        out << "x1 " << endl;
        x1.printMatrixd(out);
        out << endl;
        out << "x2 " << endl;
        x2.printMatrixd(out);
        out << endl;
        */
        Matrixd r1;              //residual
        Matrixd r2;
        r1 = b - a * x1;
        /*
        out << "r1 " << endl;
        r1.printMatrixd(out);
        out << endl;
        */
        r2 = b - a * x2;
        /*
        out << "r2 " << endl;
        r2.printMatrixd(out);
        out << endl;
        */
        Matrixd e1;              //the error
        Matrixd e2;              //the error
        e1 = ainv * r1;
        e2 = a.Gauss(r2);
        /*
        out << "e1 " << endl;
        e1.printMatrixd(out);
        out << endl;
        out << "e2 " << endl;
        e2.printMatrixd(out);
        out << endl;
        */
        double enorm1;
        double xnorm1;
        double num1 = 0;
        enorm1 = e1.norm();
        xnorm1 = x1.norm();
        num1 = enorm1 / xnorm1;      //if num < 0.5 will converge
        double enorm2;
        double xnorm2;
        double num2 = 0;
        enorm2 = e2.norm();
        xnorm2 = x2.norm();
        num2 = enorm2 / xnorm2;      //if num < 0.5 will converge
        return num2;
}

Matrixd Matrixd::iterative_refinement(Matrixd & cin, ofstream & out)
{
        //iterative refinement for ill-behaved matrices
        int n = rows;                   //number of equations
	int i = 0;
	int j = 0;
	int k = 0;

	double * x1 = new double[n];	//Solution Array
        //Matrixd xtilde_start;
        //Matrixd a1;
        //a1 = *this;
        //xtilde_start = a1.Gauss(cin);
        double dum = 0;
        for(i = 0; i < n; i++)
                x1[i] = dum;

	Matrixd oldx(n, 1, x1);		//Solution Vector
	Matrixd xtilde;
        //(n, 1, x1);	//Approximate Solution Vector
        //xtilde = xtilde_start;
        Matrixd deltax(n, 1, x1);        //intermediate solution difference
        Matrixd ctilde(n, 1, x1);        //intermediate RHS difference
        Matrixd E(n, 1, x1);             //the RHS
	Matrixd a;
	Matrixd c;

	a = *this;
	c = cin;

        //xtilde = a.Inverse()*c;         //from a * x = c
        xtilde = a.Gauss(c);

        double relerr;
        double errstop = 0.01;
        double err = 1;
        double olderr = 1.1 * err;
        double cond = 1.0;
        relerr = errstop * 1.1;

        /*
        out << "a " << endl;
        a.printMatrixd(out);
        out << endl;
        out << "c " << endl;
        c.printMatrixd(out);
        out << endl;
        out << "Starting xtilde " << endl;
        xtilde.printMatrixd(out);
        out << endl;
        */
        double denom = 1.0;
        int counter = 0;
        int maxiter = 1e4;
        while(relerr > errstop && err > 0.00001 && counter < maxiter)
        {
                counter++;
                //Start each loop by calculating current ctilde's
                //xtilde is either all 0's (for first time) or the approx
                //solution from the last iteration
                ctilde = a*xtilde;      // c~ = a x~
                /*
                out << endl;
                out << "ctilde for iteration " << counter << endl;
                ctilde.printMatrixd(out);
                out << endl;
                out << "E for iteration " << counter << endl;
                */
                E = c - ctilde;         // a * deltax = c - c~ = E
                //E.printMatrixd(out);

                err = 0;
                for(i = 0; i < n; i++)
                {
                        if(c.val(i) == 0)
                                denom = 1e-8;
                        else
                                denom = c.val(i);
                        err+=E.val(i)*E.val(i)/denom/denom;
                }
                err = sqrt(err/(double)n);
                /*
                out << endl;
                out << "err for iteration " << counter << " = " << err;
                out << endl;
                */
                if(err != 0)
                {
                        relerr = fabs((err - olderr)/err);
                        /*
                        out << endl;
                        out << "relerr for iteration " << counter << " = " << relerr;
                        out << endl;
                        */
                        //Now solve for new and improved xtilde
                        //deltax = a.Inverse() * E;
                        deltax = a.Gauss(E);
                        /*
                        out << endl;
                        out << "deltax for iteration " << counter;
                        out << endl;
                        */
                        //deltax.printMatrixd(out);
                        //cond = deltax.condition();
                        oldx = xtilde;
                        xtilde = oldx + deltax;
                        /*
                        out << endl;
                        out << "xtilde for iteration " << counter+1;
                        out << endl;
                        xtilde.printMatrixd(out);
                        */
                }
                else
                        relerr = 0.9 * errstop;
                olderr = err;
        }

        return xtilde;
}

Matrixd Matrixd::diag_dom(ofstream & out)
{
        Matrixd diag(rows, cols, 0.0);
        Matrixi row(rows, 1, 999);
        Matrixd dom_matrix;
        double temp;
        dom_matrix = *this;
        for(int r = 0; r < rows; r++)
        {
                for(int c = 0; c < cols; c++)
                {
                        temp = testrow(r,c);
                        diag.set(r,c,temp);
                }
        }
        out << endl;
        diag.printMatrixd(out);
        //Now we have diag matrix to help decide which equation should go where
        for(int i = 0; i < rows; i++)
        {
                double temp = 0;
                int index = 0;
                for( int r = 0; r < rows; r++)
                {
                        if(testforuse(row, r) == 1)
                        {
                                if(diag.val(r, i) > temp)
                                {
                                        temp = diag.val(r, i);
                                        index = r;
                                }
                                row.set(i, index);
                        }
                }
        }
        out << endl;
        row.printMatrixi(out);
        dom_matrix.setorder(row);
        //put rows in the right spots
        for(int r = 0; r < rows; r++)
        {
                for(int c = 0; c < cols; c++)
                {
                        dom_matrix.set(r, c, this->val(row.val(r), c));
                }
        }
        return dom_matrix;
}

double Matrixd::testrow(int r, int cin)
{
        double sum = 0;
        for(int c = 0; c < cols; c++)
        {
                if( c != cin)
                        sum += fabs(this->val(r,c));
        }
        double temp;
        double dom = 0;
        temp = this->val(r,cin) / sum;
        dom = temp;
        /*
        if(temp > 1)
                dom = temp;
        else
                dom = 0.0;
        */
        return dom;
}

int Matrixd::testforuse(Matrixi & row, int r)
{
        int flag = 1;
        for(int i = 0; i < rows; i++)
        {
                if(row.val(i,1) == r)
                        flag = 0;
        }
        return flag;
}

void Matrixd::setorder(Matrixi & order)
{
        order_vector = order;
}

Matrixd Matrixd::Gauss_Seidel(Matrixd & cin, ofstream & out)
{
        Matrixd x(rows, 0.0);        //Start off with all zeros
        Matrixd xold(rows, 0.0);
        int flag = 0;
        int maxiter = 1000;
        int counter = 0;
        double maxerr = 0.00001;
        double relerr = 1.1 * maxerr;
        double temp = 0.0;
        double temp2 = 0.0;
        Matrixd aorig;
        aorig = *this;
        //first let's make sure that a is diagonally dominant
        Matrixd a;
        a = diag_dom(out);
        //now a is diagonally dominant
        while(flag == 0 && counter < maxiter)
        {
                flag = 1;
                counter++;
                for(int i = 0; i < rows; i++)
                {
                        xold.set(i, x.val(i));
                        temp = 0.0;
                        for(int c = 0; c < cols; c++)
                        {
                                if( c != i)
                                        temp += a.val(i,c)*x.val(c);
                        }
                        temp2 = (cin.val(i) - temp) / a.val(i,i);
                        x.set(i, temp2);
                        relerr = fabs((x.val(i) - xold.val(i)) / x.val(i));
                        if(relerr > maxerr)
                                flag = 0;
                }
        }
        //now need to reorder x before returning it
        Matrixd xorig;
        xorig = x;      //set them equal then only make changes where you need to
        for(int i = 0; i < rows; i++)
        {
                xorig.set(x.order_vector.val(i), x.val(i));
        }
        return xorig;

}