/** @(#) matrixdouble/squaredmatrix.cpp */
/**
* Define las funciones del encabezado.
* @author Daniel Alba Cuellar, Omar Posada Villarreal
* @version 1.0, 01/03/2002
*/
#include "squaredmatrix.h"
// constructor, ~, friend-----------------------------------------------------
/** No se limpia porque se copiaran los datos. */
SquaredMatrix :: SquaredMatrix(RectangularMatrix &matrix)
: RectangularMatrix(matrix.getRows(), matrix.getRows(),
matrix.getLogic(), false) {
copyFrom(matrix);
}
/* Destructor predeterminado. */
/** Libera recursos. */
// operator-------------------------------------------------------------------
// public---------------------------------------------------------------------
/** Copia el triangulo superior actual al triangulo inferior para obtener una
* matriz simetrica rellena. No toca la diagonal.
* Logica cero (row), var(col). */
void SquaredMatrix :: upperTriangularToSymmetric() {
int i, j;
int nc = columns + logic;
for (i = 0; i < (rows - 1); ++i) {
for (j = (i + 1 + logic); j < nc; ++j) {
m[j - logic][i + logic] = m[i][j];
}
}
}
/** Forma la matriz R en la actual de la factorizacion QR,
* en base a la original A y la diagonal.
* R deberia ser Triangular superior.
* No checa dimensiones.
* Uso:
* R.createRMatrix(A, d);
* Logica uno.
* @param A Original obtenida de qrdcmp().
* @param d Diagonal */
void SquaredMatrix :: createRMatrix(SquaredMatrix &A, VectorDouble &d) {
setLogic(1);
A.setLogic(1);
d.setLogic(1);
int i, j;
clean();
// Copiar la triangular superior
for (i = 1; i <= (rows - 1) ; ++i) {
for (j = (i + 1); j <= columns; ++j) {
m[i - 1][j] = A[i][j]; //m[][] equiv R[][]
}
}
/*
for (i = 1; i <= (n - 1) ; ++i) {
for (j = (i + 1); j <= n ; ++j) {
R(i, j) = A[i][j];
}
}
*/
// Copia la diagonal
for (i = 1; i <= rows; ++i) {
m[i - 1][i] = d[i];
}
}
/** Calcula la norma max para la matriz actual.
* Logica uno.
*/
double SquaredMatrix :: norm1() {
//* Compute norm1 for M[m,n] (matrix)
setLogic(1);
double max = 0.0;
double sum = 0.0;
int i, j;
for (i = 1; i <= rows; ++i) {
sum += fabs(m[i-1][1]);
}
max = sum;
for (j = 2; j <= columns; ++j) {
sum = 0.0;
for (i = 1; i <= rows; ++i) {
sum += fabs(m[i -1][j]);
}
if (sum > max) {
max = sum;
}
}
return max;
}
/** Estima el numero de condicion "l1" de una matriz triangular superior "R",
* usando el algoritmo contenido en Cline, Moler, Sterwart, Wilkinson [1979].
*/
double SquaredMatrix :: condition() {
setLogic(1);
int n = rows;
VectorDouble x(n, 1);
VectorDouble diag(n, 1); // antes M2 Diagonal de la matriz
VectorDouble p(n, 1);
VectorDouble pm(n, 1);
double xp;
double xm;
double temp;
double tempm;
double xnorm;
int i, j;
// 1, 2.
double est;
est = norm1(); // Estimador del numero de condicion
// 3.
diag = takeDiagonal();
x[1] = 1.0 / diag[1];
// 4.
for (i = 2; i <= n; ++i) {
p[i] = m[0][i] * x[1];
}
//5.
for (j = 2; j <= n ; ++j) {
xp = (1.0 - p[j]) / diag[j];
xm = ( -1.0 - p[j]) / diag[j];
temp = fabs(xp);
tempm = fabs(xm);
// 5.5
for (i = (j + 1); i <= n; ++i) {
pm[i] += (m[j - 1][i] * xm);
tempm += fabs(pm[i]) / fabs(diag[i]);
p[i] += m[j - 1][i] * xp;
temp += fabs(p[i]) / fabs(diag[i]);
}
if (temp >= tempm) {
x[j] = xp; // ej = 1
} else {
x[j] = xm;
for (i = (j + 1); i <= n; ++i) {
p[i] = pm[i];
}
} // if
} // for
xnorm = x.norm1();
est = est / xnorm;
return est;
}
/** Calcula la cota inferior de los eigen valores de una matriz simetrica.
* <pre>
* minEigVal = min{1 <= i <= n} (aii - sum{j=1, j!=i}(|aij|) )
* </pre>
* @see Dennis. pag 60. Teorema de Gershgorin. */
double SquaredMatrix :: findMinLimitEigenValue() {
double minEig = numeric_limits<double>::max();
double sum; // == aux
int i, j;
for (i = 0; i < rows; ++i) {
sum = 0.0;
for (j = 0; j < columns; ++j) {
sum += fabs(m[i].v[j]);
}
sum = m[i].v[i] - sum - fabs(m[i].v[i]); // quitar i!=j
if (sum < minEig) {
minEig = sum;
}
}
return minEig;
}
/** Calcula la cota superior de los eigen valores de una matriz simetrica.
* <pre>
* maxEigVal = max{1 <= i <= n} (aii + sum{j=1, j!=i}(|aij|) )
* </pre>
* @see Dennis. pag 60. Teorema de Gershgorin. */
double SquaredMatrix :: findMaxLimitEigenValue() {
double maxEig = numeric_limits<double>::min();
double sum; // == aux
int i, j;
for (i = 0; i < rows; ++i) {
sum = 0.0;
for (j = 0; j < columns; ++j) {
sum += fabs(m[i].v[j]);
}
sum = m[i].v[i] + sum - fabs(m[i].v[i]); // quitar i!=j
if (sum > maxEig) {
maxEig = sum;
}
}
return maxEig;
}
/** Convierte a la matriz actual en identidad.
* Ceros excepto en la diagonal con unos. */
void SquaredMatrix :: toIdentity() {
clean();
int i;
for (i = 0; i < rows; ++i) {
m[i].v[i] = 1.0;
}
}
/** Realiza M pasos iterativos del metodo de la potencia
* para obtener el eigenvalor dominante de la matriz actual.
* Tambien incorpora un proceso de ortogonalizacion
* para encontrar los sucesivos eigenvalores.
* @author Daniel Alba Cuellar, Omar Posada Villarreal (generalizar)
* @param x Vector inicial proporcionado por el usuario.
* @param M numero de iteraciones.
* @param seq se intenta encontrar el eigenvalor numero seq. [1, n].
* @param return V Almacena en la matriz V el eigenvector generado.
* Para encontrar el de mayor valor absoluto, seq = 1.
* @param return w Eigen vector generado.
* @return Maximo eigen valor.
*/
/*
main
A = this, n = rows
// realiza el proceso iterativo para encontrar los n eigenvalores de A
for (i=1; i <= n; i++)
{
printf("eigenvector %d:\n", i);
power_d( iter, i );
// almacena en la matriz V el eigenvector generado
for (j=1; j <= n; j++)
V[j][i] = w[j];
}
*/
// Falta probarlo
double SquaredMatrix :: dominantEigenValue(VectorDouble &x, int M, int seq,
SquaredMatrix &V, VectorDouble &w)
{ int k, i, j, l;
double r; // cociente Raleigh, eigen valor, lambda
double max_y;
VectorDouble y(rows, 1);
VectorDouble v_temp(rows, 1); /* vector temporal */
double e_temp; /* escalar temporal */
//n == rows
w = x;
/* Proceso de ortogonalizacion: construye un vector inicial w a ser usado
por el metodo de la potencia */
/* w = x - suma((V(.,j)*x/V(.,j)*V(.,j))V(.,j)).
V(.,i) es el j-ésimo eigenvector de A. */
for (j=1; j < seq; j++)
{ /* extrae el j-ésimo eigenvector de la matriz V */
for (l=1; l <= rows; l++)
v_temp[l] = V[l][j];
// innerProduct <v_temp, x> / <v_temp,v_temp>
e_temp = (v_temp * x) / (v_temp * v_temp);
for (l=1; l <= rows; l++)
w[l] = w[l] - e_temp * v_temp[l];
}
/* aqui comienza el metodo de la potencia */
for (k = 1; k <= M; k++)
{ /* y = Aw */
y = (*this) * w;
/* r = w*y/w*w. esta es la aproximacion del eigenvalor por medio
del cociente de Rayleigh */
r = (w * y) / (w * w);
/* calcula la norma infinito de y */
max_y = fabs(y[1]);
for (i = 2; i <= rows; i++)
{ if (fabs(y[i]) > max_y)
{ max_y = fabs(y[i]);
}
}
/* w = y/norma(y). norma(y) es la norma infinito de y */
for (i = 1; i <= rows; i++)
{ y[i] = y[i]/max_y;
w[i] = y[i];
}
/* Proceso de ortogonalizacion: elimina de w los errores de redondeo */
/* w = x - suma((V(.,j)*x/V(.,j)*V(.,j))V(.,j)).
V(.,i) es el j-ésimo eigenvector de A. */
for (j=1; j < seq; j++)
{ /* extrae el j-ésimo eigenvector de la matriz V */
for (l=1; l <= rows; l++)
v_temp[l] = V[l][j];
e_temp = (v_temp * y)/(v_temp * v_temp);
for (l=1; l <= rows; l++)
w[l] = w[l] - e_temp * v_temp[l];
}
}
return r;
}
// private--------------------------------------------------------------------
// Fin------------------------------------------------------------------------