LU.cpp

// Resolve um sistema linear pela substituição LU.
double* LU (double **A, double *b)
{
	int i, j, k, p;
	double *pivot = new double[n];
	double Amax, t, m, r, Mult;
	// 1 PASSO: Transformar a matriz A do problema em duas matrizes triangulares L e U. 	
	for (i = 0; i < n; i++)
		pivot[i] = i;
	for (j = 0; j < n-1; j++)
	{
		cout << "Iteracao " << j << endl;
		imprimeMatriz(A,n);		
		// Escolher pivot
		p = j;
		Amax = abs(A[j][j]);
		// Verifica na coluna a ser eliminada qual elemento possui o maior valor absoluto, este elemento será o pivô.		
		for (k = j+1; k < n; k++)
		{
			if (abs(A[k][j]) > Amax)
			{
				Amax = abs(A[k][j]);
				p = k;
			}
		}
		// Se (p != j) então deve-se trocar de linhas 
		if (p != j)
		{
			for (k = 0; k < n; k++)
			{
				t = A[j][k];
				A[j][k] = A[p][k];
				A[p][k] = t;	
			}
			m = pivot[j];
			pivot[j] = pivot[p];
			pivot[p] = m;
		}
		if (abs(A[j][j]) != 0)
		{
			// Eliminação de Gauss
			r = 1 / A[j][j];
			for (i = j+1; i < n; i++)
			{
				Mult = A[i][j]*r;
				A[i][j] = Mult;
				for (k = j+1; k < n; k++)
					A[i][k] = A[i][k] - Mult*A[j][k];
			}
		}
	}
	cout << "\nDecomposicao L/U" << endl;
	imprimeMatriz(A,n);
	// A matriz A agora é L\U. Elementos abaixo da diagonal principal são de L, os da diagonal principal para cima pertencem a U.
	
	// 2 PASSO: Realizar substituições sucessivas para resolver o sistema triangular inferior: Ly = b		
	double *y = new double[n];
	double soma;
	k = pivot[0];
	y[0] = b[k];
	// Percorre somente os elementos de L em A.
	for (i = 1; i < n; i++)
	{
		soma = 0;
		for (j = 0; j <= i-1; j++)
		{
			soma += A[i][j]*y[j]; 
		}
		k = pivot[i];
		y[i] = b[k] - soma;	
	}
	// 3 PASSO: Realizar substituições retroativas para resolver o sistema triangular superior: Ux = y 			
	double *x = new double[n];
	x[n-1] = y[n-1] / A[n-1][n-1];
	for (i = n-2; i >= 0; i--)
	{
		soma = 0;
		for (j = i+1; j < n; j++)
			soma += A[i][j]*x[j];
		x[i] = (y[i] - soma) / A[i][i];
	}	

	delete [] pivot;
	delete [] y;

	imprimeSolucao(x);
	return (x);
}