Revision 1264 trunk/src/linad99/eigen.cpp
eigen.cpp (revision 1264)  

25  25 
ad_exit(1); 
26  26 
} 
27  27 
dmatrix m1=symmetrize(m); 
28 
int n=m1.rowsize(); 

29 
int cmin=m1.colmin(); 

30 
int rmin=m1.rowmin(); 

31  28 
m1.colshift(1); // set minimum column and row indices to 1 
32  29 
m1.rowshift(1); 
30 
int n=m1.rowsize(); 

33  31 
dvector diag(1,n); 
34  32 
dvector off_diag(1,n); 
35  33  
...  ...  
68  66 
ad_exit(1); 
69  67 
} 
70  68 
int n=m.rowsize(); 
71 
int l,k,j,i;


69 
int l,j,i; 

72  70 
double scale,hh,h,g,f; 
73  71  
74  72 
for (i=n;i>=2;i) 
...  ...  
77  75 
h=scale=0.0; 
78  76 
if (l > 1) 
79  77 
{ 
80 
for (k=1;k<=l;k++) 

78 
for (int k=1;k<=l;k++)


81  79 
scale += fabs(m[i][k]); 
82  80 
if (scale == 0.0) 
83  81 
e[i]=m[i][l]; 
84  82 
else 
85  83 
{ 
86 
for (k=1;k<=l;k++) 

84 
for (int k=1;k<=l;k++)


87  85 
{ 
88  86 
m[i][k] /= scale; 
89  87 
h += m[i][k]*m[i][k]; 
...  ...  
101  99 
m[j][i]=m[i][j]/h; 
102  100 
#endif 
103  101 
g=0.0; 
104 
for (k=1;k<=j;k++) 

102 
for (int k=1;k<=j;k++)


105  103 
g += m[j][k]*m[i][k]; 
106 
for (k=j+1;k<=l;k++) 

104 
for (int k=j+1;k<=l;k++)


107  105 
g += m[k][j]*m[i][k]; 
108  106 
e[j]=g/h; 
109  107 
f += e[j]*m[i][j]; 
...  ...  
113  111 
{ 
114  112 
f=m[i][j]; 
115  113 
e[j]=g=e[j]hh*f; 
116 
for (k=1;k<=j;k++) 

114 
for (int k=1;k<=j;k++)


117  115 
m[j][k] = (f*e[k]+g*m[i][k]); 
118  116 
} 
119  117 
} 
...  ...  
138  136 
for (j=1;j<=l;j++) 
139  137 
{ 
140  138 
g=0.0; 
141 
for (k=1;k<=l;k++) 

139 
for (int k=1;k<=l;k++)


142  140 
g += m[i][k]*m[k][j]; 
143 
for (k=1;k<=l;k++) 

141 
for (int k=1;k<=l;k++)


144  142 
m[k][j] = g*m[k][i]; 
145  143 
} 
146  144 
} 
...  ...  
182  180 
{ 
183  181 
dvector& d = (dvector&) _d; 
184  182 
dvector& e = (dvector&) _e; 
183 
#ifdef EIGEN_VECTORS 

185  184 
dmatrix& z = (dmatrix&) _z; 
185 
#endif 

186  186 
int max_iterations=30; 
187  187 
int n=d.size(); 
188  188 
max_iterations+=10*(n/100); 
189 
int m,l,iter,i,k;


189 
int m,l,iter,i; 

190  190 
double s,r,p,g,f,dd,c,b; 
191  191  
192  192 
for (i=2;i<=n;i++) e[i1]=e[i]; 
...  ...  
232  232 
g=c*rb; 
233  233 
/* Next loop can be omitted if eigenvectors not wanted */ 
234  234 
#ifdef EIGEN_VECTORS 
235 
for (k=1;k<=n;k++) 

235 
for (int k=1;k<=n;k++)


236  236 
{ 
237  237 
f=z[k][i+1]; 
238  238 
z[k][i+1]=s*z[k][i]+c*f; 
...  ...  
266  266 
int max_iterations=30; 
Also available in: Unified diff