unit FormaEc;

interface
uses
	xMatDefs, variable,FilUSR, horrores, cprocs, CYTipos;

procedure Chequeo;
procedure dimencione;
procedure ObtenerEcuaciones;

implementation

{$IFDEF testFormaec }
procedure muestra(piv,red:integer);
var k,j:word;
begin
	writeln(lst,'------------');
	writeln(lst,'piv = ', piv:8,'redundancia = ',red:8);
	for k:=piv to n do
	begin
		writeln(lst);
		for j:=piv to nm do
			if abs(m[k,j])>1e-20 then write(lst,'x')
			else write(lst,'o');
	end;
	writeln(lst);
end;
{$ENDIF}

procedure limpie; forward;


procedure borre(var x: TVectFila);
var k:integer;
begin
	for k:=1 to nm do x[k]:=0
end;


procedure compile;
begin
	chek:=false;
	nve:=0;
	nd:=0;
	limpie;
	cirt1;
	nvecirt1:=nve;
	cirt2;
end;

procedure dimencione;
var
	k:integer;
begin
	for k:=1 to qnx do indx[k]:=k;
	v[0]:=0;
	nm:=qnm;
	nx:=1;
	n:=qn;
	compile;
	inicio;{---}
	nx:=nve;
	ndiodos:=nd;
	n:=nn+nx;
	nm:=n+nx+nindep+1;
end;

procedure chequeo;
begin
	chek:=true;
	nd:=0;nve:=nvecirt1;
	cirt2
end;


procedure limpie;
var
	k, j:integer;
begin
	for k:=1 to n do
		borre(m[k]);
{ 20.2.98 cambio para inicializar la lista para cada ecuaci˘n
	m s que para cada nodo }
	for j:=0 to n do{	for j:=0 to nn do }
		arbv[j]:=-1;
	for k:=1 to nx do
	begin
		x[k]:=x0[indx[k]];
		indx[k]:=k
	end;
	x0:=x;
end;





procedure CopiarDerivada(var f:TVectEstReal);
var j:word;
begin
for j:= 1 to nx do
	f[nn+j]:=f[nn+j+nx];
end;




procedure escaler(var  a: TMatSistema);

const
 Cero = 1.0e-20; (* Estudiar un poco donde poner la definicion *)

var
 piv,pm:word;
 temp:NReal;
 pc1,pf1:word;


procedure ElimineRedundantes;
var
 pfr:integer;  (* indice a la primer fila redundancia *)
 pcr :integer; (* indice a primer columna de redundancia *)
 k,j:integer;

begin
	pfr:= nn+ nx-redundancia+1;
	pcr:=nn+nx+nx-redundancia+1;

   { Limpiamos los "ceros sucios" dejados por escaler 1}
	for k:= nn+2 to pfr-1 do
		for j:= nn+1 to k-1 do a[k,j]:=0;

   { Derivamos la relación entre las redundantes }
	for k:=pfr to pf1 do
	begin
		for j:= pf1 to pf1+nx do
		begin
			a[k,j-nx]:=a[k,j];
			a[k,j]:=0;
		end;

		{ la derivada de las entradas constantes es CERO }
		a[k,pf1+nx+1]:=0;

      { Detectamos si es necesario derivar las entradas }
		for j:= nn+nx+nx+2 to pc1 do
			if Not EsCero( a[k,j] ) then error('FormaEc.EliminarRedundantes: SERIA NECESARIO DERIVAR UNA ENTRADA');
	end;

	redundancia:=0;

   { retroceso del pivote }
	piv:= nn+1;
end;


function PruebaPivote:boolean;
var
	Estado:(Probando,Falso,Bueno);
begin
	Estado:=Probando;
	repeat
		if piv> pf1-redundancia then
		begin
			PruebaPivote:= false;
			exit;
		end;

      pm:=PointMaxCol(a,piv,pf1,piv);
		if piv<>pm then
			IntercaFila(a[piv],a[pm],piv,pc1);

		temp:=a[piv,piv];
      if abs(temp) < Cero then
			if piv > nn then
			begin
				Inc(redundancia);
				if piv<=pf1-redundancia then
				begin
					IntercaCol(a,piv,nn+nx-redundancia+1,1,pf1);
					IntercaCol(a,piv+nx,nn+nx+nx-redundancia+1,1,pf1);
					pm:=indx[piv-nn];
					indx[piv-nn]:=indx[nx-redundancia+1];
					indx[nx-redundancia+1]:=pm;
				end
				else Estado:=falso
			end
			else error('FormaEc.PruebaPivote: CIRCUITO INCONEXO')  {cir. inconexo}
		else Estado:=Bueno;
	until Estado<>Probando;
	PruebaPivote:=Estado=Bueno;
end;


var
	salida:boolean;

begin
{----------Escaler1------------------}
	piv:=1;
	redundancia:=0;

	pc1:=nn+2*nx+nindep+1;
	pf1:=nn+nx;

	salida:=false;

	{ Mientras se pueda realizar la escalerización directa }
	while PruebaPivote do
	begin
		Eliminacion1(a, piv, pf1, pc1);
		inc(piv)
	end;

	{ Si se produjo redundancia durante la escalerización directa
	tratar la matriz para eliminar la redundancia }
	if redundancia > 0 then
	begin
		ElimineRedundantes;
		{ Mientras se pueda realizar la escalerización directa }
		while PruebaPivote do
		begin
			Eliminacion1(a, piv, pf1, pc1);
			inc(piv)
		end;
	end;

	if redundancia >0 then error('FormaEc.escaler: NO PUEDO RESOLVER LA REDUNDANCIA.'); { no puedo resolver la redundancia }

	{ inicio de ESCA2 }
	for piv:= pf1 downto 2 do
	begin

		if abs(a[piv,piv]) > Cero then
			DivFila(a[piv],nn+nx+1,pc1,a[piv,piv])
		else
			error('FormaEc.escaler, abs(piv)= 0 en ESCA2'); { imposible }

		for pm:= 1 to piv-1 do
			if abs(a[pm,piv]) > Cero then
				SumFilaFac(a[pm],a[piv],nn+nx+1,pc1,-a[pm,piv]);
	end;

	DivFila(a[1],nn+nx+1,pc1,a[1,1]);

	{ Fin de ESCA2 }


end;{fin escaler}


procedure ObtenerEcuaciones;
begin
	Compile;
	Escaler(m);
end;


begin
end.(* FormaEc *)


