/*
 * Biblioteca de Funcoes para resolver sistemas Ax=B,
 * fazer pivoteamento simultâneo e Eliminação Simultanea em A e B
 *----------------------------------------------------------------
 * Autores:
 * - Gustavo Sverzut Barbieri <ra008849@ic.unicamp.br> - Eng. Computação - Unicamp 
 * - Ivens Prates Teles Alves <ra008908@ic.unicamp.br> - Eng. Computação - Unicamp
 */


#include "elim_gauss.h"

/*
 *Funcao: PivoteamentoParcial
 *Descrição: Faz o pivoteamento parcial das Matrizes A e B simultaneamente 
 *Autor: Gustavo Sverzut Barbieri <ra008849@ic.unicamp.br>
 *Data: 2001-10-01
 *Parametros:
 *  Matriz *A: matriz quadrada(nxn)
 *  Matriz *B: matriz linear(nx1)
 *Retorno: nada
 */
void PivoteamentoParcial(Matriz *A, Matriz *B) {
  double pivo, aux;
  int i,j,k,r;

  for (k=0; k < (A->nc -1); k++) { 
    /* Acha o pivo */
    pivo = abs(A->info[k][k]);
    r=k;
    for (i=(k+1); i< A->nc; i++) {
      if (abs(A->info[i][k]) > pivo) {
	  pivo = abs(A->info[i][k]);
	  r=i;
      }
    }
    if (r!=k) { /* Troca linhas  de 'A' e de 'B'*/
      for (j=0; j < A->nc; j++) {
	aux = A->info[k][j];
	A->info[k][j] = A->info[r][j];
	A->info[r][j] = aux;
      }
      aux = B->info[k][0];
      B->info[k][0] = B->info[r][0];
      B->info[r][0] = aux;
    }
  }
}


/*
 *Funcao: EliminacaoGauss
 *Descrição: Faz a Eliminação de Gauss em A e B simultaneamente
 *Autor: Gustavo Sverzut Barbieri <ra008849@ic.unicamp.br>
 *Data: 2001-10-01
 *Parametros:
 *  Matriz *A: matriz quadrada(nxn)
 *  Matriz *B: matriz linear(nx1)
 *Retorno: nada
 */
void EliminacaoGauss(Matriz *A, Matriz *B) {
  int i,j,k;
  double mult;
  /* Faz a eliminacao */
  for (k=0; k < (A->nl-1); k++) {
    for (i=(k+1); i < A->nl; i++) {
      mult = A->info[i][k] / A->info[k][k];
      A->info[i][k]=0;
      for (j=k+1; j < A->nl; j++)
	A->info[i][j] = A->info[i][j] - mult * A->info[k][j];
      B->info[i][0] = B->info[i][0] - mult * B->info[k][0];
    }
  }
}


/*
 *Funcao: ResolveSistemaTriangularSuperior
 *Descrição: Resolve sistema Ax=B, sendo A triangular superior
 *Autor: Gustavo Sverzut Barbieri <ra008849@ic.unicamp.br>
 *Data: 2001-10-01
 *Parametros:
 *  Matriz *A: matriz quadrada(nxn) triangular superior
 *  Matriz *B: matriz linear(nx1)
 *  char nome[]: nome a ser atribuido à matriz-solução
 *Retorno: Apontador para matriz-solução do sistema
 */
Matriz *ResolveSistemaTriangularSuperior(Matriz *A, Matriz *B, char nome[]) {
  Matriz *X = NULL;
  int k,j;
  double s=0;
  X = matCria(B->nl,1,nome);
  
  X->info[X->nl-1][0] = B->info[B->nl-1][0] / A->info[A->nl-1][A->nc-1];
  
  for (k = (X->nl - 2); k >= 0; k--) {
    for (j=(k+1), s=0; j < X->nl; j++)
      s = s + A->info[k][j]*X->info[j][0];
    X->info[k][0] = (B->info[k][0] - s)/A->info[k][k];
  }
  return X;
}


/*
 *Funcao: Residuo
 *Descrição: Retorna o residuo na solução X em função de A e B
 *Autor: Ivens Prates Telles Alves <ra008908@ic.unicamp.br>
 *Data: 2001-10-01
 *Parametros:
 *  Matriz *A: matriz quadrada(nxn)
 *  Matriz *X: matriz linear(nx1)
 *  Matriz *B: matriz linear(nx1)
 *Retorno: resíduo na solução
 */
double Residuo(Matriz *A, Matriz *X, Matriz *B) {
  int i,j;
  double residuo=0, som=0;

  for (i=0;i<A->nl;i++) {
    som = 0;
    for (j=0;j<A->nc;j++)
      som += (A->info[i][j] * X->info[j][0]);
    som -= B->info[i][0];
    residuo += som*som;
  }
  return(sqrt(residuo));
}