/*
 *  LinTool.cpp
 *
 *  Created by Assistant on 8/28/06.
 *  Copyright 2006 __MyCompanyName__. All rights reserved.
 *
 */
#include <iostream>


#include "LinTool.h"
#include <stdlib.h>
#include <ctime> 

using namespace std;

#define DIM_ERROR 1

bool seedYes = true;

LT_matrix::LT_matrix(int Nr, int Nc)
{
			
	int iR, iC;
	
	Nrows = Nr;
	Ncols = Nc;
	
	A = (double**) malloc(Nrows*sizeof(double));
	
	for(iR=0; iR<Nrows; iR++)
	{
		A[iR] = (double*) malloc(Ncols*sizeof(double));
	}
	
	//init with 0.0s:
	for (iR=0; iR<Nrows ; iR++)
	{
		for (iC=0; iC<Ncols ; iC++)
		{
			A[iR][iC]=0.0;
		}
	}
}

LT_matrix::~LT_matrix()
{
	int iR;
	for(iR=0; iR<Nrows; iR++)
	{
		free(A[iR]);
	}
	free(A);
}







void LT_matrix::solve(LT_matrix* bb, LT_matrix* xx)
{
	// Gauss Elimination
	// Transformation auf Dreiecksmatrix
	
	int i,j,k,iR,iC,n;              // Zaehlvariablen
	double f;               // Faktor zur Zeilenaddition
	double** a;
	double* b;
	double* x;
	
	if (Nrows != Ncols)
	{
		cout << "LinTool: error using solve: A in Ax=b must be square\n";
		exit(DIM_ERROR);
	}
	if (Nrows != bb->Nrows || Nrows != xx->Nrows)
	{
		cout << "LinTool: error using solve: x and b in Ax=b must be column vectors of same length as size of A\n";
		exit(DIM_ERROR);
	}
	if (bb->Ncols!=1 || xx->Ncols!=1)
	{
		cout << "LinTool: error using solve: x and b in Ax=b must be column vectors\n";
		exit(DIM_ERROR);
	}
	n=Nrows;
	
	// Copy A, b
	a = (double**) malloc(n*sizeof(double));
	for(iR=0; iR<n; iR++)
	{
		a[iR] = (double*) malloc(n*sizeof(double));
		for (iC=0; iC<Ncols ; iC++)
		{
			a[iR][iC]=A[iR][iC];
		}
	}
	
	b = (double*) malloc(n*sizeof(double));
	for(iR=0; iR<n; iR++) b[iR]=bb->A[iR][0];
	
	x = (double*) malloc(n*sizeof(double)); 
	
	
	
	for (i=0; i<n-1; i++){                      // n-1 Eliminationsschritte
	//  write(a);
	for (j=i+1; j<n; j++)                    // jeweils alle folgenden Zeilen bearbeiten
	{
		f = - a[j][i] / a[i][i];               // Faktor bestimmen
		a[j][i] = 0.0;                         // dieses Element wird zu 0 gemacht
		for (k=i+1; k<n; k++)                  // zum Rest der j-ten Zeile das ...
			a[j][k] = a[j][k] + f * a[i][k];     // ... f-fache der i-ten Zeile addieren
		b[j] = b[j] + f * b[i];                // ebenso auf der rechten Seite
	}
	}
	//write(a);
	// x[i] aus Dreiecksmatrix bestimmen
	for (i=n-1; i>=0; i--)
	{
		x[i] = b[i];                             // rechte Seite
		for (j=n-1; j>i; j--)
			x[i] = x[i] - a[i][j] * x[j];          // alle bekannten x[j] einsetzen
		x[i] = x[i] / a[i][i];                   // durch Diagonalelement dividieren
	}
	
	// Copy x back
	for(iR=0; iR<n; iR++) xx->A[iR][0]=x[iR];
	
	//free memory:
	for(iR=0; iR<n; iR++) free(a[iR]);
	free(a);
	free(b);
	free(x);
}







//overload the [][] access
double* LT_matrix::operator[](int i)
{
	return A[i];
}



std::ostream& operator<<(std::ostream &s, LT_matrix &M)
{
	int iR, iC;
	s << M.Nrows <<" rows, "<< M.Ncols <<" columns\n";
	
	for (iR=0; iR<M.Nrows ; iR++)
	{
		for (iC=0; iC<M.Ncols ; iC++)
		{
			 cout << M.A[iR][iC] << "\t";
		}
		cout << endl;
	}
	return s;
}



//Sum of all matrix elements
double LT_matrix::ElemSum()
{
	int iR, iC;
	double sumRes=0.0;
	for (iR=0; iR<Nrows ; iR++)
	{
		for (iC=0; iC<Ncols ; iC++)
		{
			 sumRes += A[iR][iC];
		}
	}
	return sumRes;
}



/*
LT_matrix LT_matrix::operator+(LT_matrix B)
{
	LT_matrix C(Nrows, Ncols);
	int iR, iC;
	for (iR=0; iR<Nrows ; iR++)
	{
		for (iC=0; iC<Ncols ; iC++)
		{
			 C.A[iR][iC] = A[iR][iC] + B.A[iR][iC];
		}
	}
	return C;
}

*/

void LT_matrix::copy( LT_matrix* N)
{
	int iR, iC;
	if (Nrows != N->Nrows || Ncols != N->Ncols)
	{
		cout << "LinTool: error using copy: must have same size. Quiting programm."<< endl;
		exit(DIM_ERROR);
	}
	for (iR=0; iR<Nrows ; iR++)
	{
		for (iC=0; iC<Ncols ; iC++)
		{
			A[iR][iC] = N->A[iR][iC];
		}
	}
}

void  LT_matrix::add(LT_matrix* M, LT_matrix* N)
{
	int iR, iC;
	if (Nrows != M->Nrows || Nrows != N->Nrows || Ncols != M->Ncols || Ncols != N->Ncols)
	{
		cout << "LinTool: error using add: all 3 matrix dimensions must agree. Quiting programm."<< endl;
		exit(DIM_ERROR);
	}
	
	for (iR=0; iR<Nrows ; iR++)
	{
		for (iC=0; iC<Ncols ; iC++)
		{
			A[iR][iC] = M->A[iR][iC] + N->A[iR][iC];
		}
	}
}

void  LT_matrix::sub(LT_matrix* M, LT_matrix* N)
{
	int iR, iC;
	if (Nrows != M->Nrows || Nrows != N->Nrows || Ncols != M->Ncols || Ncols != N->Ncols)
	{
		cout << "LinTool: error using sub: all 3 matrix dimensions must agree. Quiting programm."<< endl;
		exit(DIM_ERROR);
	}
	
	for (iR=0; iR<Nrows ; iR++)
	{
		for (iC=0; iC<Ncols ; iC++)
		{
			A[iR][iC] = M->A[iR][iC] - N->A[iR][iC];
		}
	}
}

void  LT_matrix::mult_elemwize(LT_matrix* M, LT_matrix* N)
{
	int iR, iC;
	if (Nrows != M->Nrows || Nrows != N->Nrows || Ncols != M->Ncols || Ncols != N->Ncols)
	{
		cout << "LinTool: error using mult_elemwize: all 3 matrix dimensions must agree. Quiting programm."<< endl;
		exit(DIM_ERROR);
	}
	
	for (iR=0; iR<Nrows ; iR++)
	{
		for (iC=0; iC<Ncols ; iC++)
		{
			A[iR][iC] = M->A[iR][iC] * N->A[iR][iC];
		}
	}
}

void  LT_matrix::div_elemwize(LT_matrix* M, LT_matrix* N)
{
	int iR, iC;
	if (Nrows != M->Nrows || Nrows != N->Nrows || Ncols != M->Ncols || Ncols != N->Ncols)
	{
		cout << "LinTool: error using div_elemwize: all 3 matrix dimensions must agree. Quiting programm."<< endl;
		exit(DIM_ERROR);
	}
	
	for (iR=0; iR<Nrows ; iR++)
	{
		for (iC=0; iC<Ncols ; iC++)
		{
			A[iR][iC] = M->A[iR][iC] / N->A[iR][iC];
		}
	}
}




void LT_matrix::mult( LT_matrix* M, LT_matrix* N)
{
	int iR, iC, k;
	double **P, **Q;
	if (Nrows != M->Nrows || Ncols != N->Ncols || M->Ncols != N->Nrows)
	{
		cout << "LinTool: error using mult: all 3 matrix dimensions must agree in A = B*C. Quiting programm."<< endl;
		exit(DIM_ERROR);
	}
	
	P=M->A;
	Q=N->A;
	
	for (iR=0; iR<Nrows ; iR++)
	{
		for (iC=0; iC<Ncols ; iC++)
		{
			A[iR][iC] = 0.0; //M->A[iR][iC] / N->A[iR][iC];
			for (k=0; k<M->Ncols; k++) A[iR][iC] += P[iR][k] * Q[k][iC];
		}
	}
}

void LT_matrix::mult(double t)
{
	int iR, iC;
	for (iR=0; iR<Nrows ; iR++)
	{
		for (iC=0; iC<Ncols ; iC++)
		{
			A[iR][iC] *= t;
		}
	}
}



LT_matrix* LT_matrix::transp(void)
{
	int iR, iC;
	
	double** At = (double**) malloc(Ncols*sizeof(double));
	
	for(iR=0; iR<Ncols; iR++)
	{
		At[iR] = (double*) malloc(Nrows*sizeof(double));
	}
	
	//copy:
	for (iR=0; iR<Nrows ; iR++)
	{
		for (iC=0; iC<Ncols ; iC++)
		{
			At[iC][iR] = A[iR][iC];
		}
	}
	
	//free:
	for(iR=0; iR<Nrows; iR++)
	{
		free(A[iR]);
	}
	free(A);
	
	//exchange:
	A=At;
	iR=Nrows;
	Nrows=Ncols;
	Ncols=iR;
	
	return this;
}


//Fill marix with randoms
void LT_matrix::rand_fill(double min, double max)
{
	int iR, iC;
	double mrand = RAND_MAX;
	double val;
	
	if (seedYes)
	{
		srand((unsigned)time(0));	// seed the generator with system time
		seedYes = false;
	}
	
	for (iR=0; iR<Nrows ; iR++)
	{
		for (iC=0; iC<Ncols ; iC++)
		{
			val = double(rand())/mrand + double(rand())/(mrand*mrand) + double(rand())/(mrand*mrand*mrand);
			val *= (max-min);
			val -= min;
			A[iR][iC] = val;
		}
	}
}


bool LT_matrix::is_this(LT_matrix* M)
{
	return (M == this);
}







//news:

//Find maximal value in matrix
double LT_matrix::max()
{
	int iR, iC;
	double max = A[0][0];
	
	for (iR=0; iR<Nrows ; iR++)
	{
		for (iC=0; iC<Ncols ; iC++)
		{
			if (A[iR][iC] > max) max = A[iR][iC];
		}
	}
	return max;
}

//Find minimum value in matrix
double LT_matrix::min()
{
	int iR, iC;
	double min = A[0][0];
	
	for (iR=0; iR<Nrows ; iR++)
	{
		for (iC=0; iC<Ncols ; iC++)
		{
			if (A[iR][iC] < min) min = A[iR][iC];
		}
	}
	return min;
}

