/*
Aminur Rahman
Class:  Math 631 Graduate Linear Algebra
Instructor:  Dr. Johnathan Luke

This program has been created to implement
operations on matrices and implement
Gaussian elimination to reduce matrices to
reduced row echelon form.  I will soon write
code in Java to implement these methods, and
put it up on my website as a java applet.


Please do not copy and submit as your own,
for this goes against NJIT's honor code.
Thank you.
*/

#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include <string.h>
#include <fstream.h>
#include <math.h>

double **matrix_alloc(int, int);
void *smalloc(unsigned);
void matrix_write(FILE *,double **, int, int);
void scalar_matrix_mult(double, double **, double **, int, int);
void matrix_add(double **, double **, double **, int, int);
void matrix_mult(double **A, int rA, int cA, double **B, int rB, int cB, double **C,int rC, int cC);
void add_mult(double s, double **A, int i, int j, int c);
void exchange(double **A, int i, int j);
void row_mult(double s, double **A, int i, int c);
void GE(double **A, int rA, int cA);

main()
{
	double **A, **B, **C;
    int rA,rB,rC;
    int cA,cB,cC;

    rA = rC = 3;
	cA = rB = 3;
	cB = cC = 5;

    A = matrix_alloc(rA,cA);
    B = matrix_alloc(rB,cB);
    C = matrix_alloc(rC,cC);

    A[0][0] = 1.0;
	A[0][1] = 2.0;
	A[0][2] = 4.0;
	//A[0][3] = 4.0;
	//A[0][4] = 5.0;
	A[1][0] = 2.0;
	A[1][1] = 4.0;
	A[1][2] = 2.0;
	//A[1][3] = 9.0;
	//A[1][4] = 10.0;
	A[2][0] = 4.0;
	A[2][1] = 2.0;
	A[2][2] = 1.0;
	//A[2][3] = 14.0;
	//A[2][4] = 15.0;
	//A[3][0] = 16.0;
	//A[3][1] = 17.0;
	//A[3][2] = 18.0;
	//A[3][3] = 19.0;
	//A[3][4] = 20.0;
	//A[4][0] = 21.0;
	//A[4][1] = 22.0;
	//A[4][2] = 23.0;
	//A[4][3] = 24.0;
	//A[4][4] = 25.0;


    matrix_write(stdout,A,rA,cA);
    //matrix_write(stdout,B,rB,cB);
    //scalar_matrix_mult(12.0,A,C,rA,cA);
    //matrix_write(stdout,C,rC,cC);
    //matrix_add(A,B,C,rA,cA);
	//matrix_mult(A,rA,cA,B,rB,cB,C,rC,cC);
	//matrix_write(stdout,C,rC,cC);
	//add_mult(2.0,A,0,1,cA);
	//matrix_write(stdout,A,rA,cA);
	//exchange(C,0,1);
	//matrix_write(stdout,C,rC,cC);
	//row_mult(2,B,1,cB);
    //matrix_write(stdout,B,rB,cB);
	GE(A,rA,cA);
	matrix_write(stdout,A,rA,cA);
	return 0;
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void GE(double **A, int rA, int cA){
	/*  Implements Gaussian Elimination on matrix A  */

	int j=0;//rows of A

	//puts it into row echelon form
	for(int i=0; i<cA; i++){//columns of A

		double *colA = new double[cA];
		double *temp = new double[cA];//the zero vector
		int t=j;
		for(int k=0; t<rA; k++){
			temp[k]=0;
			colA[k] = A[t][i];
			t++;
		}

		if(colA==temp)
			continue;//if the column is equivalent to the zero vector then iterates loop
		

		else{

			if(A[j][i] != 0){//elimination
				for(int k=j+1; k<rA; k++){
					double s = -(A[k][i])/(A[j][i]);
					add_mult(s,A,k,j,cA);//eliminates k using j
					double r = 1/(A[j][i]);
					row_mult(r,A,j,cA);
				}
				for(int l=j-1; l>-1; l--){
						double s = -(A[l][i])/(A[j][i]);
						add_mult(s,A,l,j,cA);//eliminates l using k
					}
				j++;//going to next row once one column is finished
			}

			else if((i != (cA-1)) && (j != (rA-1))){//exchange
				for(int k=j+1; k<rA; k++){
					if(A[k][i] != 0)
						exchange(A,j,k);
				}
				if(A[j][i] != 0)
					i--;
			}
		}
	}
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

double **matrix_alloc(int r, int c)
{
    /* Allocates data space and row pointers for a rxc matrix  */

    double **A;
    int i,j;

    A = (double **) smalloc(r*sizeof(double *));

    for(i=0; i<r; i++)  A[i] = (double *) smalloc(c*sizeof(double));

    for(i=0; i<r; i++)
    {
        for(j=0; j<c; j++)  A[i][j] = 0.0;
    }

    return(A);
}


void matrix_write(FILE *fp, double **A, int r, int c)
{
    /*  Writes A in ascii format to file pointed to by fp  */
		ofstream myfile;
		myfile.open("Math631Proj2.txt",ios::app);
    int i,j;

    for(i=0; i<r; i++)
    {
        for(j=0; j<c; j++){  fprintf(fp,"%8.3lg ",A[i][j]);
		myfile << A[i][j] << " ";
		}
		fprintf(fp,"\n");
        myfile << "\n";
    } 
	fprintf(fp,"\n");
    myfile << "\n";   
	myfile.close();
}


void scalar_matrix_mult(double alpha, double **A, double **B, int r, int c)
{
    /*  Returns alpha*A in B  */

    int i,j;

    for(i=0; i<r; i++)
    { 
        for(j=0; j<c; j++)  B[i][j] = alpha*A[i][j];
    }
}

void matrix_add(double **A, double **B, double **C, int r, int c)
{
    /*  Returns the sum of A and B in C  */

    int i, j;

    for(i=0; i<r; i++)
    {
        for(j=0; j<c; j++)  C[i][j] = A[i][j] + B[i][j];
    }
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
void matrix_mult(double **A, int rA, int cA, double **B, int rB, int cB, double **C,int rC, int cC)
{
    /* Returns the matrix product of A and B in C  */
	if((cA != rB) || (rA != rC) || (cB != cC)){
		cout<<"Error:  Please Check Dimensions.  Thank You.\n";
		exit(1);
	}
	else{
		for(int i=0;i<rA;i++){//rows of A and rows of C
			for(int j=0;j<cB;j++){//cols of B and cols of C
				for(int k=0;k<cA;k++){//cols of A and rows of B
				C[i][j] += (A[i][k])*(B[k][j]);
				}
			}
		}
	}
}

void add_mult(double s, double **A, int i, int j, int c){
	for(int l=0; l<c; l++){
		A[i][l] += s*A[j][l];
	}
}
void exchange(double **A, int i, int j){
	double *rp;
	rp=A[i];
	A[i]=A[j];
	A[j]=rp;
}
void row_mult(double s, double **A, int i, int c){
	for(int l=0; l<c; l++){
		A[i][l]=s*A[i][l];
	}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void *smalloc(unsigned b)
{
    void *p;

    if((p=malloc(b))==NULL)
     {
	 fprintf(stderr,"Memory allocation error\n");
	 exit(1);
     }

     return(p);
}