// sistemas.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "math.h"
#include "conio.h"
#include "malloc.h"

// definição de um tipo matriz
typedef struct matriz {
	int linha;
	int coluna;
	float *aij;
} matriz;

// define para acesso aos coeficientes da matriz]
#define coef(m,a,b) ((m).aij[((a)-1)*((m).coluna)+(b)-1])

// a matriz
matriz m3_3, m3_3_aumentada;

void matriz_initial(void)
{
	m3_3.linha=3;
	m3_3.coluna=3;

	m3_3.aij=(float *)malloc(sizeof(float)*9);
	m3_3.aij[0]=2;
	m3_3.aij[1]=3;
	m3_3.aij[2]=-1;
	m3_3.aij[3]=-1;
	m3_3.aij[4]=1;
	m3_3.aij[5]=4;
	m3_3.aij[6]=1;
	m3_3.aij[7]=-8;
	m3_3.aij[8]=-1;

	m3_3_aumentada.linha=3;
	m3_3_aumentada.coluna=4;

	m3_3_aumentada.aij=(float *)malloc(sizeof(float)*12);
	m3_3_aumentada.aij[0]=2;
	m3_3_aumentada.aij[1]=3;
	m3_3_aumentada.aij[2]=-1;
	m3_3_aumentada.aij[3]=0;

	m3_3_aumentada.aij[4]=-1;
	m3_3_aumentada.aij[5]=1;
	m3_3_aumentada.aij[6]=4;
	m3_3_aumentada.aij[7]=3;

	m3_3_aumentada.aij[8]=1;
	m3_3_aumentada.aij[9]=-8;
	m3_3_aumentada.aij[10]=-1;
	m3_3_aumentada.aij[11]=1;
}

float determinante(matriz *m)
{
	int sinal, i, j, k;
	float res=0;
	matriz subm;
	int dim_matriz;

	// a matriz nao e quadrada, nao tem determinante
	if(m->coluna!=m->linha)
		return 0;
	dim_matriz=m->linha;
	// matriz a uma dimensao, o determinante e o unico coeficiente
	if(dim_matriz==1)
		return (res=coef(*m,1,1));
	// matriz a dimensao 2, o determinante e diretamente calculavel
	else if(dim_matriz==2)
		return (res=coef(*m,1,1)*coef(*m,2,2)-coef(*m,1,2)*coef(*m,2,1));
	// matriz a dimensao >2, produto dos coeficientes de uma linha com os sub-determinantes
	else {

		// preparação do espaço para uma matriz de dimensão dim_matriz-1
		subm.coluna=subm.linha=dim_matriz-1;
		subm.aij=(float *)malloc(sizeof(float)*(dim_matriz-1)*(dim_matriz-1));

		// para cada coeficiente da primeira linha
		for(i=1,sinal=1;i<=dim_matriz;i++,sinal*=-1) {
			// preparaçao da submatriz, eliminando a primeira linha e i-ésima coluna
			for(j=2;j<=dim_matriz;j++)
				for(k=1;k<=dim_matriz-1;k++) {
					coef(subm,j-1,k)=coef(*m,j,(k>=i)?k+1:k);
				}
			// o resultado é acressentado do novo termo: subdeterminante pelo coeficiente e o sinal
			res+=sinal*coef(*m, 1, i)*determinante(&subm);
		}
		
		free(subm.aij);
		return res;
	}	
}


// recebe a matriz aumentada (cada linha tem n+1 elementos)
void eliminacao_gauss(matriz *m)
{
	float max, val, factor;
	int i, j, k, l_max;

	// eliminaçao das variaveis,
	// da primeira a penultima coluna de variavel (sem considerar os termos constantes)
	for(i=1;i<=m->coluna-2;i++) {
		// procure a linha do coeficiente maior
		max=0;
		for(j=i;j<=m->linha;j++)
			if((val=fabs(coef(*m,j,i)))>max) {
				max=val;
				l_max=j;
			}
		// nao tem coeficiente nao nulo, melhor parar
		if(max==0)
			return;
		// a linha do maior coeficiente nao nulo é a primeira das linhas a modificar
		// a linha i e a linha do maior coeficiente sao trocadas
		if(l_max!=i) {
			for(k=1;k<=m->coluna;k++) {
				val=coef(*m,i,k);
				coef(*m,i,k)=coef(*m,l_max,k);
				coef(*m,l_max,k)=val;
			}
		}
		// para as linhas de i ate a ultima
		for(j=i+1;j<=m->linha;j++) {
			factor=coef(*m,j,i)/coef(*m,i,i);
			coef(*m,j,i)=0;
			for(k=i+1;k<=m->coluna;k++) {
				coef(*m,j,k)-=factor*coef(*m,i,k);
			}
		}
	}
}

// a matriz quadrada a triangularizar e devolve a matriz de transformação dos termos constantes
void eliminacao_gauss_multisolucao(matriz *m, matriz *m_b)
{
	int dim_matriz, l_max, i, j, k;
	float val, max, factor;

	if(m->coluna!=m->linha)
		return;
	dim_matriz=m->coluna;

	// criação da matriz de transformação dos termos constantes, 
	// e inicialização coma matriz identidade
	m_b->coluna=m_b->linha=dim_matriz;
	m_b->aij=(float *)malloc(sizeof(float)*dim_matriz*dim_matriz);
	for(i=1;i<=dim_matriz;i++)
		for(j=1;j<=dim_matriz;j++)
			if(i==j)
				coef(*m_b,i,j)=1;
			else
				coef(*m_b,i,j)=0;

	// eliminaçao das variaveis,
	// da primeira a penultima coluna de variavel
	for(i=1;i<=dim_matriz-1;i++) {
		// procure a linha do coeficiente maior
		max=0;
		for(j=i;j<=m->linha;j++)
			if((val=fabs(coef(*m,j,i)))>max) {
				max=val;
				l_max=j;
			}
		// nao tem coeficiente nao nulo, melhor parar
		if(max==0)
			return;
		// a linha do maior coeficiente nao nulo é a primeira das linhas a modificar
		// a linha i e a linha do maior coeficiente sao trocadas
		// mesma coisa para a matriz de transformação
		if(l_max!=i) {
			for(k=1;k<=dim_matriz;k++) {
				val=coef(*m,i,k);
				coef(*m,i,k)=coef(*m,l_max,k);
				coef(*m,l_max,k)=val;

				val=coef(*m_b,i,k);
				coef(*m_b,i,k)=coef(*m_b,l_max,k);
				coef(*m_b,l_max,k)=val;
			}
		}
		// para as linhas de i ate a ultima
		for(j=i+1;j<=m->linha;j++) {
			factor=coef(*m,j,i)/coef(*m,i,i);
			coef(*m,j,i)=0;
			for(k=i+1;k<=dim_matriz;k++) {
				coef(*m,j,k)-=factor*coef(*m,i,k);
			}
			// modificação dos termos da matriz de transformação
			for(k=1;k<=dim_matriz;k++) {
				coef(*m_b,j,k)-=factor*coef(*m_b,i,k);
			}
		}
	}
	
}

// recebe a matriz trianguçar aumentada (cada linha tem n+1 elementos)
float *solucao_sistema_triangular(matriz *m)
{
	float *s, soma;
	int i, j;

	s=(float *)malloc(sizeof(float)*(m->linha));
	for(i=m->linha;i>=1;i--) {
		// inicializaçao da soma com o termo constante
		soma=coef(*m,i,m->coluna);
		// determinação
		for(j=i+1;j<=m->coluna-1;j++)
			soma-=coef(*m,i,j)*s[j-1];
		s[i-1]=soma/coef(*m,i,i);
	}
	return s;
}

void print_matriz(matriz *m)
{
	int i,j;

	for(i=1;i<=m->linha;i++) {
		for(j=1;j<=m->coluna;j++)
			cprintf("%.2f, ", coef(*m,i,j));
		cprintf("\r\n");
	}
}

int main(int argc, char* argv[])
{
	matriz m_trans;
	int i;
	float *solucao;

	matriz_initial();

	print_matriz(&m3_3);
	cprintf("determinante: %f\r\n", determinante(&m3_3));
	cprintf("\r\n");

	print_matriz(&m3_3_aumentada);
	eliminacao_gauss(&m3_3_aumentada);
	cprintf("depois da eliminacao gaussiana\r\n");
	print_matriz(&m3_3_aumentada);
	solucao=solucao_sistema_triangular(&m3_3_aumentada);
	cprintf("solucao:\r\n");
	for(i=0;i<m3_3_aumentada.linha;i++)
		cprintf("%.2f, ", solucao[i]);
	cprintf("\r\n");

	eliminacao_gauss_multisolucao(&m3_3, &m_trans);
	cprintf("depois da eliminacao gaussiana, com produção da matriz de transformação\r\n");
	print_matriz(&m3_3);
	cprintf("matriz de transformação\r\n");
	print_matriz(&m_trans);

	while(!getch());
	return 0;
}