# Un programa para calcular los tres eigenvalores 
# menores de una matrix de tamano n x n, y el correspondiente
# eigenvector (o vector singular) tal que min(li) y li > 1e-8, para i=1,2,3,
# donde li es cada uno de los eigenvalores calculados
#
# Ideas tomadas de:
# [1] G.H. Golun and C.F. Van Loan, Matrix Computations, 3rd edition
#   chapter 8.5, page 439 
# [2] W. Barth, R.S. Martin and J.H. Wilkinson, Calculation of the
#   eigenvalues of a symmetric tridiagonal matrix by the method of
#   bisection, Numerische Mathematik 9, pp. 386-393, 1967
#
#  Fraga 24/01/2014
#
# 1. Se calcula B = A^T * A
# 2. Se reduce B a la forma tridiagonal
# 3. Se hallan los extremos de los valores de
#    todos los eigenvalores y se aplica [2]
##### Ok al 25/01/2014
# 4. Se hallan los tres eigenvalores más pequeños
# 5. Se busca el eigenvalor más pequeños y mayor a 1e-8, este es l
# 6. Se calcula el eigenvector correspondiente
#    al eigenvalor más pequeño por el método
#    de iteración inversa, usando la matriz B:
#     (B-lI)x = b0
# 
#    b = [1, 1, ..., 1]^T; b0 = b/norm(b);
#    Se resuelve usando la descomposición QR de (B-lI)   
#
##### Ok al 28/01/2014
## 21/05/2014 Corregido un bug en la función biseccion()
use strict;

# Realiza la transformación Householder
sub house( $$$$ )
{
	my $A = shift;
	my $n = shift;
	my $k = shift;
	my $v = shift;

	# Cambiado: house( A( k+1:n-1, k )
    # v(0) = A( k+1, k )
    # v(1) = A( k+2, k )
	# v(3) = A( k+3, k ) 
    #
    # v(m) = A( n-1, k )
	#
	# Si v(0,...,n-1)
	# Si v(k+1,...,n-1)

	my $sigma = 0;
	for( my $i=$k+2; $i<$n; $i++ ) {
		$sigma += $A->[$i][$k] * $A->[$i][$k];
	}
	$v->[$k+1] = 1;
	for( my $i=$k+2; $i<$n; $i++ ) {
		$v->[$i] = $A->[$i][$k];
	}
	my $beta;
	if ( $sigma < 1e-10 ) {
		$beta = 0;
	}
	else {
		my $mu = sqrt( $A->[$k+1][$k]*$A->[$k+1][$k] + $sigma );
		if ( $A->[$k+1][$k] <= 0.0 ) {
			$v->[$k+1] = $A->[$k+1][$k] - $mu;
		}
		else {
			$v->[$k+1] = -$sigma/($A->[$k+1][$k] + $mu);
		}
		my $val = $v->[$k+1]*$v->[$k+1]; 
		$beta = 2.0*$val/($sigma+$val); 
		$val = $v->[$k+1];
		for( my $i=$k+1; $i<$n; $i++ ) {
			$v->[$i] /= $val;
		}
	}
	# for( my $i=$k+1; $i<$n; $i++ ) {
	# 	print "$v->[$i] ";
	# }
	# print "\nbeta= $beta\nk= $k\n";
	# print "n = $n\n";
	return( $beta );
}

# Rutina que reduce la matriz A, de tamaño n x n,
# a una matriz tridiagonal. La matriz de salida se sobreescribe en A
sub tridiagonal( $$ )
{
	my $A = shift;
	my $n = shift;

	my $v = [];
	my $p = [];
	my $w = [];
	for( my $k=0; $k<$n-2; $k++ ) {
		#                                    n = 5,  m =  n - (k+1) +1 = n-k
		# Original: house( A( k+1:n , k )  si k=1, m = 5 - 1 = 4
		# Cambiado: house( A( k+1:n-1, k ) si k=0, m = n-1-(k+1)+1 = n-k-1 = 4
		my $beta = house( $A, $n, $k, $v );

		for( my $i=$k+1; $i<$n; $i++ ) { 
			my $suma = 0;
			for( my $j=$k+1; $j<$n; $j++ ) { 
				$suma += $A->[$i][$j] * $v->[$j];
			}
			$p->[$i] = $beta*$suma;
		}
		my $val = 0;
		for( my $j=$k+1; $j<$n; $j++ ) { 
			$val += $p->[$j] * $v->[$j];
		}
		$val = $val * $beta / 2.0;

		for( my $j=$k+1; $j<$n; $j++ ) { 
			$w->[$j] = $p->[$j] - $val * $v->[$j]; 
		}

		my $suma = 0;
		for( my $j=$k+1; $j<$n; $j++ ) { 
			$val = $A->[$j][$k];
			$suma += $val * $val;
		}
		$val = sqrt( $suma );

		# si k=2
		$A->[$k+1][$k] = $val;
		$A->[$k][$k+1] = $val;

		for( my $i=$k+1; $i<$n; $i++ ) { 
			for( my $j=$k+1; $j<$n; $j++ ) { 
				$A->[$i][$j] = $A->[$i][$j] - $v->[$i]*$w->[$j] - $w->[$i]*$v->[$j];
			}
		}
	}
}

# Rutina que multiplica A^T * A
# Solo es necesario calcular la porción superior (o inferior)
# de la matriz dado que el resultado es una matriz simétrica
# Args: S = E^T * E, y n, el tamaño de las matrices es n x n
sub multiplicaATA( $$$ )
{
	my $S = shift;
	my $E = shift;
	my $n = shift;

	for( my $i=0; $i<$n; $i++ ) {
		for( my $j=$i; $j<$n; $j++ ) {
			my $suma = 0;
			for( my $k=0; $k<$n; $k++ ) {
				$suma += $E->[$k][$i] * $E->[$k][$j];
			}
			$S->[$i][$j] = $suma;
			$S->[$j][$i] = $suma;
		}
	}
}

# Rutina que lee la matriz de entrada almacenada en un archivo
# Args: nombre_archivo referencia_a_la_matriz_de_salida
sub lee_archivo( $$ ) 
{
	my $arch = shift;
	my $M = shift;

	open( my $ARCH, '<', $arch ) or die "No se puede abrir archivo $arch\n";

	my $n = <$ARCH>;	
	chop( $n );

	my $i=0;
	while( <$ARCH> ) {
		my @d = split;
		if ( $#d > 1 ) {
			for( my $j=0; $j<$n; $j++ ) {
				$M->[$i][$j] = $d[$j];
			}
			$i++;
		}
	}
	close( $ARCH );
	return( $n );
}

# Imprime una matrix A de tamaño n x n
sub imprime( $$ )
{
	my $A = shift;
	my $n = shift;

	for( my $i=0; $i<$n; $i++ ) {
		for( my $j=0; $j<$n; $j++ ) {
			print "$A->[$i][$j] ";
		}
		print "\n";
	}
	print "\n";
}

# Rutina que encuentra los eigenvalores de una matriz tridiagonal
# Args: M, la matriz tridiagonal de entrada
#       n, el tamaño de la matriz de n x n
#      eigen, referencia para almacenar los eigenvalores
#      m1, m2, los eigenvalores deseados, l_m1 < l_m1+1 < l_m1+2 < ... < l_m2
sub biseccion( $$$$$ )
{
	my $M = shift;
	my $n = shift;
	my $eigen = shift;
	my $m1 = shift;
	my $m2 = shift;

	my @beta;
	$beta[0] = 0;
	for( my $i=$n-1; $i>=1; $i-- ){
		$beta[$i] = $M->[$i][$i-1]*$M->[$i][$i-1];
	}
	# Calculo of xmin y xmax
	my $k;
	$k = $n-1;
	my $xmin;
	my $xmax;
	$xmin = $M->[$k][$k] - abs( $M->[$k][$k-1] );
	$xmax = $M->[$k][$k] + abs( $M->[$k][$k-1] );

	for( my $i=$n-2; $i>0; $i-- ){
		my $h = abs($M->[$i][$i-1]) + abs($M->[$i+1][$i]);
		my $v = $M->[$i][$i] + $h;
		if ( $v > $xmax ) {
			$xmax = $v;
		}
		$v = $M->[$i][$i] - $h;
		if ( $v < $xmin ) {
			$xmin = $v;
		}
	}
	print "$xmin $xmax\n";
	my $DBL_EPSILON = 2.22045E-16;
	my $eps1=1e-8;
	my $eps2;
	if( $xmin + $xmax > 0 ) {
		$eps2 = $DBL_EPSILON*$xmax;
	}
	else {
		$eps2 = $DBL_EPSILON*(-$xmin);
	}
	if( $eps1 <= 0 ) {
		$eps1 = $eps2;
	}
	$eps2 = 0.5*$eps1 + 7.0*$eps2;

	my $x0 = $xmax;
	my @wu;
	for( my $i=$m1; $i<=$m2; $i++ ) {
		$eigen->[$i] = $xmax;
		$wu[$i] = $xmin;
	}
	my $z = 0;

	# Ciclo para calcular el k-ésimo eigenvalor
	for( $k=$m2; $k>=$m1; $k-- ) {
		my $xu = $xmin;
		CONTINUA: for( my $i=$k; $i>=$m1; $i-- ) {
			if( $xu < $wu[$i] ) {
				$xu = $wu[$i];
				last CONTINUA;
			}
		}
		if( $x0 > $eigen->[$k] ) {
			$x0 = $eigen->[$k];
		}
	
		while(  $x0-$xu > 2.0*$DBL_EPSILON*(abs($xu)+abs($x0))+$eps1 ) {
			$z++;

			my $x1 = ($xu+$x0)/2.0;
			my $a = 0;
			my $q = 1;
			for( my $i=0; $i<$n; $i++ ) {
				my $v;
				if ( abs($q) > 1e-10 ) {
					$v = $beta[$i]/$q;
				}
				else {
					if ( $i==0 ){
						$v = 0;
					}
					else{
						$v = abs( $M->[$i-1][$i]/$DBL_EPSILON );
					}
				}
				$q = $M->[$i][$i] - $x1 - $v; 
				if( $q<0 ) {
					$a++;
				}
			}
			if( $a<$k ) {
				if( $a<$m1 ) {
					$xu = $wu[$m1] = $x1;
				}
				else{
					$xu = $wu[$a+1] = $x1;
					if( $eigen->[$a] > $x1 ){
						$eigen->[$a] = $x1;
					}
				}
			}
			else {
				$x0 = $x1;
			}
		}
		$eigen->[$k] = ($x0+$xu)/2.0; 
	}
	$eigen->[0] = $z;
}

# Subrutina que calcula descomposición QR
# A = QR, A es de tamaño m x n
# No retorna el valor de  Q
sub descomposicionQR( $$$$ )
{
	my $A = shift;
	my $m = shift;
	my $n = shift;
	my $R = shift;
  
	my $vec = []; 
	for( my $k=0; $k<$n; $k++ ) {

		my $sum = 0;
		for( my $i=0; $i<$m; $i++ ) {
			$sum += $A->[$i][$k] * $A->[$i][$k];
		}
		$sum = sqrt( $sum );
		$R->[$k][$k] = $sum;
		for( my $i=0; $i<$m; $i++ ) {
			$vec->[$i] = $A->[$i][$k] / $sum;
		}
		for( my $j=$k+1; $j<$n; $j++ ) {
			$sum = 0;
			for( my $i=0; $i<$m; $i++ ) {
				$sum += $vec->[$i] * $A->[$i][$j];
			}
			$R->[$k][$j] = $sum;

			for( my $i=0; $i<$m; $i++ ) {
				$A->[$i][$j] -= $vec->[$i] * $sum;
			}
		}
	}
}

# Si se tiene el sistema Ax = b, A es de tamaño n x n, 
# se realiza la matriz aumentada [A|b] y se realiza
# la descomposición QR, se obtiene Q[ R | Q^T b ]
# Esta rutina obtiene x a partir de la matriz Rp = [ R | Q^T b ]
sub resuelve( $$$ )
{
	my $Rp = shift;
	my $n = shift;
	my $x = shift;

	my $j = $n-1;

	$x->[$j] = $Rp->[$j][$n] / $Rp->[$j][$j];

	for( my $i=$n-2; $i>=0; $i-- ) {
		my $suma = 0;
		for( $j=$i+1; $j<$n; $j++ ) {
			$suma += $Rp->[$i][$j] * $x->[$j];
		}
		$x->[$i] = ( $Rp->[$i][$n] - $suma )/$Rp->[$i][$i];
	}

}


die "Args: archivo\n" if $#ARGV != 0;

my $arch = $ARGV[0];

my $A = [];
my $B = [];
my $n = lee_archivo( $arch, $A );

imprime( $A, $n );

multiplicaATA( $B, $A, $n );
imprime( $B, $n );

# Hago una copia de la matriz original
my $Bcopia = [];
for( my $i=0; $i<$n; $i++ ) {
	for( my $j=0; $j<$n; $j++ ) {
		$Bcopia->[$i][$j] = $B->[$i][$j];
	}
}

tridiagonal( $B, $n );
imprime( $B, $n );

# Se calculan los tres eigenvalores más pequenos.
my $eigen = [];

biseccion( $B, $n, $eigen, 1, 3 );
print "Iteraciones: $eigen->[0]\n";
for( my $i=1; $i<=3; $i++ ) {
	print "eigen $i: $eigen->[$i]\n";
}

# Selecciona el más pequeño que sea mayor a 1e+8
my $cual = 0;
AQUI: for( my $i=1; $i<=3; $i++ ) {
	if( abs($eigen->[$i]) > 1e-8 ){
		$cual = $i;
		last AQUI;
	}
}
die "Los eigenvalores seleccionados son todos iguales a cero\n" if $cual == 0;

# Se calcula el eigenvector asociado al eigenvalor encontrado
# Se usa la matriz original y la descomposición QR

print "e: $eigen->[$cual]\n";

for( my $i=0; $i<$n; $i++ ) {
	$Bcopia->[$i][$n] = 0.5;
	$Bcopia->[$i][$i] -= $eigen->[$cual];
}
my $R = [];
descomposicionQR( $Bcopia, $n, $n+1, $R );

for( my $i=0; $i<$n; $i++ ) {
	for( my $j=0; $j<=$n; $j++ ) {
		print "$i,$j $R->[$i][$j] ";
	}
	print "\n";
}

my $evec = [];
resuelve( $R, $n, $evec );

print "r: ";
for( my $i=0; $i<$n; $i++ ) {
	print "$evec->[$i] ";
}
print "\n";

my $suma=0;
for( my $i=0; $i<$n; $i++ ) {
	$suma += $evec->[$i]*$evec->[$i];
}
$suma = sqrt( $suma );

print "rnorm: ";
for( my $i=0; $i<$n; $i++ ) {
	$evec->[$i] /= $suma;
	print "$evec->[$i] ";
}
print "\n";
