1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
|
/******************************************************************************
// FUNCION QUE CALCULA LOS AUTOVALORES DE UNA MATRIZ (LOS DEVUELVE EN UN ARRAY)
// PARA IMPLEMENTAR LA FUNCIÓN TENER EN CUENTA LO SIGUIENTE :
// (1) la función atan(x) devuelve el arco tangente de x
//*****************************************************************************/
Array1D< real > jacobi (Array2D< real >&A, // matriz original
Array2D< real >&Autovectores, // matriz de salida con los autovectores ordenados por columna
const int NMaxIter, // número máximo de iteraciones
const real TOL, // tolerancia para detener el algoritmo
int &Niter) // variable de salida con el numero de iteraciones
{
/** COMPROBAMOS QUE LAS DIMENSIONES SON COHERENTES */
if(A.dim1()!=A.dim2() || A.dim1()!=Autovectores.dim1() || A.dim2()!=Autovectores.dim2() ){
printf("jacobi() : matrix dimensions are different\n");
return(Array1D< real >());
}
/** COMPROBAMOS QUE LA MATRIZ A ES SIMÉTRICA */
int dim=A.dim1();
for(int i=0;i<dim;i++){
for(int j=i+1;j<dim;j++){
if(A[i][j]!=A[j][i]){
printf("jacobi() : matrix is not symmetric\n");
return(Array1D< real >());
}
}
}
/** HACER ALUMNO */
/* Bloque(1):Inicialización de variables*/
// hacemos una copia de la matriz
Array2D< real > A2=A.copy();
// inicializamos la matriz de autovectores a la identidad
Autovectores=0.;
for(int i=0;i<dim;i++) Autovectores[i][i]=1.;
// hacemos las iteraciones del metodo de Jacobi
int n,i,j,p,q;
real R,C,T,CO,SI,D;
/* Fin Bloque(1)*/
/* Bloque(2):Algoritmo de Jacobi y cálculo de autovectores*/
for (Niter=0;Niter<=NMaxIter;Niter++){
p=1;
q=0;
R = fabs(A2[p][q]);
/* Bloque(2.1) Cálculo del elemento mayor no perteneciente a la diagonal*/
for (i=2;i<dim;i++){
for (j=0;j<=i-1;j++){
if (fabs(A2[i][j])>R){
R=fabs(A2[i][j]);
p=j;
q=i;
}
}
} /* Fin Bloque(2.1)*/
/* Bloque(2.2): Solución: Si el elemento más grande es menor que la tolerancia
dada se considera que el algortimo ha finalizado*/
if (R<TOL)
{
Array1D< real > a(dim);
for(int l=0;l<dim;l++) a[l]=A2[l][l];
return (a);
} /* Fin Bloque(2.2) */
/* Bloque(2.3): Cálculo del seno y coseno de la matriz de rotación*/
C=(A2[q][q]-A2[p][p])/(2*A2[p][q]);
if (C<0){ T=-C-sqrtl(1.+C*C);}
else{ T=-C+sqrtl(1.+C*C); }
CO=1./sqrtl(1.+T*T);
SI=CO*T; /* Fin Bloque 2.3) */
/* Bloque(2.4): Modificación de la matriz A*/
for (j=0;j<dim;j++){
if(j!=p && j!=q){
D=A2[p][j];
A2[j][p]=CO*D-SI*A2[q][j];
A2[p][j]=A2[j][p];
A2[j][q]=SI*D+CO*A2[q][j];
A2[q][j]=A2[j][q];
}
}
A2[p][p]=A2[p][p]-T*A2[p][q];
A2[q][q]=A2[q][q]+T*A2[p][q];
A2[p][q]=0;
A2[q][p]=0; /* Fin Bloque(2.4)*/
/* Bloque(2.5): Modificación de la matriz de autovectores*/
real b_ip,b_iq;
for(i=0;i<dim;i++){
b_ip=Autovectores[i][p];
b_iq=Autovectores[i][q];
Autovectores[i][p]=CO*b_ip - SI*b_iq;
Autovectores[i][q]=SI*b_ip + CO*b_iq;
}/* Fin Bloque(2.5)*/
} /* Fin Bloque(2)*/
printf ("Error, numero de iteraciones excedido\n");
return(Array1D< real >());
}
|
