1 /* pdnmesh - a 2D finite element solver
2 Copyright (C) 2001-2005 Sarod Yatawatta <sarod@users.sf.net>
3 This program is free software; you can redistribute it and/or modify
4 it under the terms of the GNU General Public License as published by
5 the Free Software Foundation; either version 2 of the License, or
6 (at your option) any later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16 $Id: eig.c,v 1.12 2005/04/18 08:29:59 sarod Exp $
17 */
18
19 #ifdef HAVE_CONFIG_H
20 #include <config.h>
21 #endif
22
23
24 /* aigenvalue routines for symmetric matrices, storage is reduced */
25 #include <stdio.h>
26 #include <math.h> /* for sqrt */
27 #include <stdlib.h> /* for calloc */
28 #include <string.h> /* for memcpy */
29 #include "types.h"
30
31
32 MY_DOUBLE *eig_array; /* array to store solved eigenvalues */
33
34
35 /* macros to access elements in matrices */
36 /* all these macros assume (i,j) within [0,m-1] */
37 /* Upper Hessenberg form
38 * is represented by shifting all rows to the left so no zeros
39 * are present */
40 /* get the (i,j) element from upper Hessenberg A- m by m */
41 #define eH(A,i,j,m) \
42 ((((i)==0)||((i)==1))?A[(i)][(j)]:\
43 ((j)<(i)-1?0:A[(i)][(j)-(i)+1]))
44 /* tridiagonal form
45 * is represented by shifting all rows to the left until no zero appears */
46 /* maximum row length is 3 */
47 #define eT(A,i,j,m) \
48 ((((j)>(i)+1)||((j)<(i)-1))?0:\
49 (((i)==0)?A[(i)][(j)]:A[(i)][(j)-(i)+1]))
50 /* upper tridiagonal form */
51 /* is represented by shifting all row to left until no zero appears */
52 /* maximum row length is 3 */
53 #define eUT(A,i,j,m) \
54 ((((j)<(i))||((j)>(i)+2))?0:\
55 ((i)==0?A[(i)][(j)]:A[(i)][(j)-(i)+1]))
56
57
58 static void
convert_to_hessenberg(MY_DOUBLE ** A,MY_DOUBLE ** B,int m)59 convert_to_hessenberg(MY_DOUBLE **A, MY_DOUBLE **B, int m)
60 {
61 /* convert matrix A (m by m), symmetric to
62 * tridiagonal form B.
63 * structure of A - m by m square
64 * structure if B - tri diagonal
65 * both matrices modified
66 */
67 int i,j,k;
68
69 MY_DOUBLE *v,*u; /* vectors */
70 MY_DOUBLE norm;
71
72
73 v =(MY_DOUBLE*)calloc(m-1,sizeof(MY_DOUBLE));
74 if ( v == 0 ) {
75 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
76 exit(1);
77 }
78 u =(MY_DOUBLE*)calloc(m,sizeof(MY_DOUBLE));
79 if ( u == 0 ) {
80 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
81 exit(1);
82 }
83
84
85 /* main loop */
86 for (k=0; k < m; k++ ) {
87 /* copy vector */
88 norm =0;
89 for ( i=0; i<m-k-1; i++) {
90 v[i] = A[i+k+1][k];
91 /* also find the norm */
92 norm += v[i]*v[i];
93 }
94
95 /* store the first element of v temporarily */
96 u[0] = v[0];
97 /* form the reflector */
98 if ( v[0] > 0 )
99 v[0] += sqrt(norm);
100 else
101 v[0] -= sqrt(norm);
102
103 /* find the new norm */
104 norm -= u[0]*u[0];
105 norm += v[0]*v[0];
106 norm = sqrt(norm);
107 /* make vector v unit norm */
108 for (i=0; i<m-k-1; i++)
109 v[i] /= norm;
110
111 /* now form the product v'*A[k+1:m-1][k:m-1] */
112 /* and store in u */
113 for (i=0; i<m-k; i++ ) {
114 u[i] = 0;
115 for ( j=0; j <m-k-1; j++)
116 u[i] += v[j]*A[k+1+j][k+i];
117 }
118
119 /* now update A */
120 /* A = A -2v*u' */
121 /* note: v = m-k-1 by 1 */
122 /* u = m-k by 1 */
123 for (i=0; i < m-k-1; i++)
124 for (j=0; j<m-k; j++ )
125 A[k+1+i][k+j] -= 2*v[i]*u[j];
126
127 /* form the product A[0:m-1][k+1:m-1]*v */
128 /* size m by 1 */
129 for (i=0; i<m; i++) {
130 u[i] = 0;
131 for (j=0; j<m-k-1;j++)
132 u[i] += A[i][j+k+1]*v[j];
133 }
134
135 /* update A */
136 /* A = A - 2*u*v'*/
137 /* v - m-k-1 by 1 */
138 /* u - m by 1 */
139 for (i=0; i<m; i++)
140 for (j=0; j<m-k-1; j++)
141 A[i][j+k+1] -= 2*u[i]*v[j];
142 }
143
144 /* copy the result to tri diagonal B */
145 B[0][0]=A[0][0]; B[0][1]=A[0][1];
146 for (i=1;i<m-1;i++) {
147 B[i][0]=A[i][i-1];
148 B[i][1]=A[i][i];
149 B[i][2]=A[i][i+1];
150 }
151 B[m-1][0]=A[m-1][m-2]; B[m-1][1]=A[m-1][m-1];
152
153 free(u);
154 free(v);
155 }
156
157 static int
matrix_mult(MY_DOUBLE ** A,MY_DOUBLE ** B,MY_DOUBLE ** C,int m)158 matrix_mult(MY_DOUBLE **A, MY_DOUBLE **B,MY_DOUBLE **C,int m)
159 {
160 /* multiplies the two matrices */
161 /* C = A*B */
162 /* A - upper tri diagonal - main, super, super-super diagonals
163 * B - upper Hessenberg
164 * C - Result is tridiagonal, will be copied back to A
165 */
166 int i,j,k;
167
168 /* first find rows 1 to m-3, where we have 3 elements in rows of A */
169 for (i=1; i<m-2; i++) {
170 for (j=i-1;j<i+2;j++) {
171 C[i][j-i+1]=0;
172 for (k=i;k<=i+2;k++) /* only three terms */
173 C[i][j-i+1] +=eUT(A,i,k,m)*eH(B,k,j,m);
174 }
175 }
176 /* now the first row, only two elements in row of C */
177 i=0;
178 for (j=i;j<i+2;j++) {
179 C[i][j-i]=0;
180 for (k=i;k<=i+2;k++) /* only three terms */
181 C[i][j-i] +=eUT(A,i,k,m)*eH(B,k,j,m);
182 }
183 /* now the m-2 row, only 2 elements in row of A */
184 i=m-2;
185 for (j=i-1;j<i+2;j++) {
186 C[i][j-i+1]=0;
187 for (k=i;k<=i+1;k++) /* only two terms */
188 C[i][j-i+1] +=eUT(A,i,k,m)*eH(B,k,j,m);
189 }
190 /* now the m-1 row, only 1 element in row of A, 2 elements in row of C */
191 i=m-1;
192 for (j=i-1;j<i+1;j++) {
193 /* only one term */
194 C[i][j-i+1] =eUT(A,i,i,m)*eH(B,i,j,m);
195 }
196 /* copy C back to A */
197 A[0][0]=C[0][0]; A[0][1]=C[0][1];
198 for (i=1;i<m-1;i++) {
199 A[i][0]=C[i][0];
200 A[i][1]=C[i][1];
201 A[i][2]=C[i][2];
202 }
203 i=m-1;
204 A[i][0]=C[i][0];A[i][1]=C[i][1];
205 return(0);
206 }
207
208 static int
qr_decomp(MY_DOUBLE ** A,MY_DOUBLE ** Q,int m)209 qr_decomp(MY_DOUBLE **A, MY_DOUBLE **Q, int m)
210 {
211 /* perform a QR decompositon of A, A=Q.R
212 * R will replace A
213 * Q will store Q
214 */
215 /* A - tridiagonal
216 * Q - Upper Hessenberg
217 * R - Upper tridiagonal
218 * so A (tridiagonal) will change to A (Upper tridiagonal)
219 * hence we have 4 diagonals in A
220 */
221
222 /* not optimized, especially in callculating Q, work
223 * directly using A instead of vectors u and v
224 * e being a canonical basis vector is not considered */
225 /* when calculating R, no need to calculate last column (1 element only) */
226
227 MY_DOUBLE *v,*u,**scratch; /* vectors */
228 MY_DOUBLE norm;
229 int i,k,l;
230
231 v =(MY_DOUBLE*)calloc(m,sizeof(MY_DOUBLE));
232 if ( v == 0 ) {
233 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
234 exit(1);
235 }
236 u =(MY_DOUBLE*)calloc(m,sizeof(MY_DOUBLE));
237 if ( u == 0 ) {
238 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
239 exit(1);
240 }
241 scratch= (MY_DOUBLE **) calloc(2, sizeof(MY_DOUBLE *));
242 if ( scratch == 0 ) {
243 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
244 exit(1);
245 }
246 scratch[0] = ( MY_DOUBLE * ) calloc( (size_t)m, sizeof(MY_DOUBLE));
247 if ( scratch[0] == 0 ) {
248 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
249 exit(1);
250 }
251 scratch[1] = ( MY_DOUBLE * ) calloc( (size_t)m, sizeof(MY_DOUBLE));
252 if ( scratch[1] == 0 ) {
253 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
254 exit(1);
255 }
256
257
258 for (k=0; k<m-1; k++) {
259 /* x=A_[k:m,k] */
260 norm=0;
261 /* make u=0 */
262 memset((void*)u,0,sizeof(MY_DOUBLE)*(size_t)(m-k));
263 u[0]=eT(A,k,k,m); norm+=u[0]*u[0];
264 u[1]=eT(A,k+1,k,m); norm+=u[1]*u[1];
265 /*replaces code: for (i=0; i<m-k; i++) {
266 u[i]=eT(A,i+k,k,m);
267 norm+=u[i]*u[i];
268 }*/
269
270 memset((void*)v,0,sizeof(MY_DOUBLE)*(size_t)(m-k));
271 /* v_k=sign(x[0])||x|| e_1 +x */
272 v[0]=u[0];
273 v[0]+=(u[0]>=0?sqrt(norm):-sqrt(norm));
274 /* again calculate norm of v, not u */
275 norm +=v[0]*v[0]-u[0]*u[0];
276 norm=sqrt(norm);
277 v[0]/=norm;
278 /* v_k=v_k/||v_k|| */
279 v[1]=u[1]/norm;
280 /* all other elements are zero */
281 #ifdef DEBUG
282 printf("qr: Householder vector for col %d\n",k);
283 for (i=0;i<m-k;i++) {
284 printf("v[%d]="MDF"\n",i,v[i]);
285 }
286 #endif
287
288 /* A_[k:m,k:m]=A_[k:m,k:m]-2v_k(v_k'A_[k:m,k:m]) */
289 /* first form the product (v_k'A_[k:m,k:m]) */
290 /* and store it in u */
291 /* product vector has only 3 non zero elements, max */
292 memset((void*)u,0,sizeof(MY_DOUBLE)*(size_t)(m-k));
293 if ( k !=m-2) {
294 u[0]=v[0]*eT(A,k,k,m)+v[1]*eT(A,k+1,k,m);
295 u[1]=v[0]*eT(A,k,k+1,m)+v[1]*eT(A,k+1,k+1,m);
296 u[2]=v[1]*eT(A,k+1,k+2,m);
297 } else { /* k==m-2 */
298 u[0]=v[0]*eT(A,k,k,m)+v[1]*eT(A,k+1,k,m);
299 u[1]=v[0]*eT(A,k,k+1,m)+v[1]*eT(A,k+1,k+1,m);
300 }
301 #ifdef DEBUG
302 printf("qr: product vector for col %d\n",k);
303 for (i=0;i<m-k;i++) {
304 printf("u[%d]="MDF"\n",i,u[i]);
305 }
306 #endif
307
308 /*replaces code: for (i=0; i< m-k; i++) {
309 u[i]=0;
310 for (j=0; j<m-k; j++) {
311 u[i] +=v[j]*A[j+k][i+k];
312 }
313 }*/
314
315 /* A_[k:m,k:m]=A_[k:m,k:m]-2v_k( u' ) */
316 if (k!=m-2) {
317 /* only 6 elements of the matrix is affected */
318 /*A[k][k]-=2*v[0]*u[0];
319 A[k][k+1]-=2*v[0]*u[1];
320 A[k][k+2]=-2*v[0]*u[2]; *//* should be equality */
321
322 /*A[k+1][k]-=2*v[1]*u[0]; this is identically zero
323 A[k+1][k+1]-=2*v[1]*u[1];
324 A[k+1][k+2]-=2*v[1]*u[2]; */
325 /* ugly representation of above code */
326 if ( k==0 ) {
327 A[k][0]-=2*v[0]*u[0];
328 A[k][1]-=2*v[0]*u[1];
329 A[k][2]=-2*v[0]*u[2];
330
331 A[k+1][0]-=2*v[1]*u[0];
332 A[k+1][1]-=2*v[1]*u[1];
333 A[k+1][2]-=2*v[1]*u[2];
334 } else { /* k == m-2 and k==m-1 will never happen here */
335 A[k][1]-=2*v[0]*u[0];
336 A[k][2]-=2*v[0]*u[1];
337 A[k][3]=-2*v[0]*u[2];
338
339 A[k+1][0]-=2*v[1]*u[0];
340 A[k+1][1]-=2*v[1]*u[1];
341 A[k+1][2]-=2*v[1]*u[2];
342 }
343 } else { /* k==m-2 */
344 /* A[k][k]-=2*v[0]*u[0];
345 A[k][k+1]-=2*v[0]*u[1];
346 A[k+1][k]-=2*v[1]*u[0]; this is identically zero
347 A[k+1][k+1]-=2*v[1]*u[1]; */
348 A[k][1]-=2*v[0]*u[0];
349 A[k][2]-=2*v[0]*u[1];
350
351 A[k+1][0]-=2*v[1]*u[0];
352 A[k+1][1]-=2*v[1]*u[1];
353 }
354
355 /* replaces code: for (i=0; i< m-k; i++ ){
356 for (j=0; j<m-k; j++) {
357 A[i+k][j+k] -=2*v[i]*u[j];
358 }
359 } */
360
361 #ifdef DEBUG
362 for (i=0;i<m;i++) {
363 printf("| ");
364 if ( i==0 ) {
365 printf(""MDF" "MDF" "MDF"",A[i][0], A[i][1], A[i][2]);
366 for (l=3;l<m;l++)
367 printf(""MDF" ",0.0);
368 } else if (i==m-2) {
369 for (l=0;l<m-3;l++)
370 printf(""MDF" ",0.0);
371 printf(""MDF" "MDF" "MDF"",A[i][0], A[i][1], A[i][2]);
372 } else if (i==m-1) {
373 for (l=0;l<m-2;l++)
374 printf(""MDF" ",0.0);
375 printf(""MDF" "MDF"",A[i][0], A[i][1]);
376 } else {
377 for (l=0;l<i-1;l++)
378 printf(""MDF" ",0.0);
379 printf(""MDF" "MDF" "MDF" "MDF"",A[i][0], A[i][1], A[i][2], A[i][3]);
380 for (l=i+2;l<m-1;l++)
381 printf(""MDF" ",0.0);
382 }
383 printf("|\n");
384 }
385 #endif
386 /* save vector v (non zero part) in two arrays,
387 * since only 2 non zero elements
388 * to construct Q later */
389 scratch[0][k]=v[0];
390 scratch[1][k]=v[1];
391 }
392
393 /* note A */
394 /* for the k th column, only one element */
395 /* so vector is 1 */
396 /* note- we have to change sign in last column of Q now */
397 scratch[0][m-1]=1;
398 scratch[1][m-1]=0;
399
400 #ifdef DEBUG
401 printf("qr_decomp: R found\n");
402 for (i=0;i<m;i++) {
403 printf("| ");
404 if ( i==0 ) {
405 printf(""MDF" "MDF" "MDF"",A[i][0], A[i][1], A[i][2]);
406 for (l=3;l<m;l++)
407 printf(""MDF" ",0.0);
408 } else if (i==m-2) {
409 for (l=0;l<m-3;l++)
410 printf(""MDF" ",0.0);
411 printf(""MDF" "MDF" "MDF"",A[i][0], A[i][1], A[i][2]);
412 } else if (i==m-1) {
413 for (l=0;l<m-2;l++)
414 printf(""MDF" ",0.0);
415 printf(""MDF" "MDF"",A[i][0], A[i][1]);
416 } else {
417 for (l=0;l<i-1;l++)
418 printf(""MDF" ",0.0);
419 printf(""MDF" "MDF" "MDF" "MDF"",A[i][0], A[i][1], A[i][2], A[i][3]);
420 for (l=i+2;l<m-1;l++)
421 printf(""MDF" ",0.0);
422 }
423 printf("|\n");
424 }
425 printf("|| ");
426 for (i=0; i<m;i++) {
427 printf(""MDF" ",scratch[0][i]);
428 }
429 printf("||\n");
430 printf("|| ");
431 for (i=0; i<m;i++) {
432 printf(""MDF" ",scratch[1][i]);
433 }
434 printf("||\n");
435 #endif
436
437 /* now construct Q from stored values of v */
438 for (l=0;l<m;l++) {
439 /* unit vector e_l */
440 memset((void*)u,0,sizeof(MY_DOUBLE)*(size_t)(m));
441 u[l]=1;
442
443 #ifdef DEBUG
444 printf("unit vector %d\n",l);
445 for (i=0;i<m;i++) {
446 printf(""MDF" ",u[i]);
447 }
448 printf("\n");
449 #endif
450
451 for (k=m-1;k>=0; k--) {
452 /* retrieve v_k */
453 memset((void*)v,0,sizeof(MY_DOUBLE)*(size_t)(m-k));
454 v[0]=scratch[0][k];
455 v[1]=scratch[1][k];
456 #ifdef DEBUG
457 printf("using householder %d\n",k);
458 for (i=0;i<m-k;i++) {
459 printf("i=%d,"MDF" ",i,v[i]);
460 }
461 printf("\n");
462 #endif
463
464 /* form product */
465 norm=0;
466 norm+=v[0]*u[k];
467 if (k==m-1) {
468 u[k] -=2*v[0]*norm;
469 } else {
470 norm+=v[1]*u[k+1];
471 u[k] -=2*v[0]*norm;
472 u[k+1] -=2*v[1]*norm;
473 }
474 }
475 #ifdef DEBUG
476 printf("col %d of Q\n",l);
477 for (i=0;i<m;i++) {
478 printf(""MDF" ",u[i]);
479 }
480 printf("\n");
481 #endif
482 /* copy to Q, Q is upper hessenberg */
483 if ( l != m-1 ) {
484 for (i=0; i<l+2;i++) {
485 /*Q[i][l]=u[i]; */
486 if ( i<2 ) {
487 Q[i][l]=u[i];
488 } else {
489 Q[i][l-i+1]=u[i];
490 }
491 }
492 } else {
493 /* change sign of last column of Q see note A*/
494 for (i=0; i<l+1;i++) {
495 /* Q[i][l]=-u[i]; */
496 if ( i<2 ) {
497 Q[i][l]=-u[i];
498 } else {
499 Q[i][l-i+1]=-u[i];
500 }
501 }
502 }
503
504 }
505
506 #ifdef DEBUG
507 printf("Q found\n");
508 for (i=0;i<m;i++) {
509 printf("| ");
510 for (l=0;l<m;l++)
511 printf(""MDF" ",eH(Q,i,l,m));
512 printf("|\n");
513 }
514 #endif
515 free(v);
516 free(u);
517 free(scratch[0]);
518 free(scratch[1]);
519 free(scratch);
520 return(0);
521 }
522
523 static int
find_eigenvalues_tridiagonal(MY_DOUBLE ** H,MY_DOUBLE * v,int m)524 find_eigenvalues_tridiagonal(MY_DOUBLE **H,MY_DOUBLE *v,int m)
525 {
526 /* finds the eigenvalues of the matrix A (in tridiagonal form) and
527 * stores them in array v
528 * A will be modified
529 */
530
531 MY_DOUBLE **Q,**C,mu;
532 int i,k,pivot;
533 MY_DOUBLE *vt;
534 #ifdef DEBUG
535 int j;
536 #endif
537 #ifdef DEBUG
538 printf("find_eig_tri: solving for size m=%d\n",m);
539 printf("H=\n");
540 printf(MDF" "MDF"\n",H[0][0],H[0][1]);
541 for (i=1;i<m-1;i++) {
542 printf(MDF" "MDF" "MDF"\n",H[i][0],H[i][1],H[i][2]);
543 }
544 printf(MDF" "MDF"\n",H[m-1][0],H[m-1][1]);
545 printf("fing_eig_tri:=================\n");
546 for (i=0;i<m;i++) {
547 printf("find_eig_tri: eigenvalue array %d="MDF"\n",i,v[i]);
548 }
549 #endif
550
551 /* FIXME - if m is 2, we used closed form formula to find
552 * the eigenvalues */
553 if ( m==2 ) {
554 #ifdef DEBUG
555 printf("find_eig_tri: analytical solution, m=2\n");
556 #endif
557 mu=(H[0][0]-H[1][1])*(H[0][0]-H[1][1])+4*H[0][1]*H[1][0];
558 if ( mu>=0 ) {
559 mu=sqrtf((float)mu);
560 }
561 v[0]=0.5*(H[0][0]+H[1][1]+mu);
562 v[1]=0.5*(H[0][0]+H[1][1]-mu);
563 if ( ABS(v[0])<ABS(v[1]) ) {
564 mu=v[0];
565 v[0]=v[1];
566 v[1]=mu;
567 }
568 #ifdef DEBUG
569 for (i=0;i<m;i++) {
570 printf("find_eig_tri: eigenvalue array %d="MDF"\n",i,v[i]);
571 }
572 #endif
573 return(0);
574 }
575
576
577 /* Q is an Upper Hessenberg matrix */
578 Q= (MY_DOUBLE **) calloc(m, sizeof(MY_DOUBLE *));
579 if ( Q == 0 ) {
580 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
581 exit(1);
582 }
583 Q[0] = ( MY_DOUBLE * ) calloc( (size_t)m, sizeof(MY_DOUBLE));
584 if ( Q[0] == 0 ) {
585 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
586 exit(1);
587 }
588 for (i=1;i<m;i++) {
589 Q[i] = ( MY_DOUBLE * ) calloc( (size_t)(m-i+1), sizeof(MY_DOUBLE));
590 if ( Q[i] == 0 ) {
591 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
592 exit(1);
593 }
594 }
595 /* C is a temporary matrix, tri diagonal */
596 C= (MY_DOUBLE **) calloc(m, sizeof(MY_DOUBLE *));
597 if ( C == 0 ) {
598 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
599 exit(1);
600 }
601 for (i = 1; i < m-1; ++i) {
602 C[i] = ( MY_DOUBLE * ) calloc( (size_t)3, sizeof(MY_DOUBLE));
603 if ( C[i] == 0 ) {
604 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
605 exit(1);
606 }
607 }
608 C[0] = ( MY_DOUBLE * ) calloc( (size_t)2, sizeof(MY_DOUBLE));
609 if ( C[0] == 0 ) {
610 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
611 exit(1);
612 }
613 C[m-1] = ( MY_DOUBLE * ) calloc( (size_t)2, sizeof(MY_DOUBLE));
614 if ( C[m-1] == 0 ) {
615 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
616 exit(1);
617 }
618
619 /* FIXME improvements
620 * 1) shifts
621 * 2) divide to submatrices whenever off diagonal becomes almost 0
622 */
623 for (k=0; k< 8; k++) {
624 /* find new shift */
625 mu=H[m-1][1]; /* last diagonal element, but we use special storage */
626 /* substract the shift from H */
627 /* H = H-mu.I */
628 for (i=1;i<m; i++) {
629 /* H[i][i]-=mu; *
630 * but special storage for H */
631 H[i][1]-=mu;
632 }
633 H[0][0]-=mu;
634 #ifdef DEBUG
635 printf("find_eig_tri:iteration %d mu="MDF" \n",k,mu);
636 #endif
637 qr_decomp(H,Q,m); /* H=Q.H */
638 #ifdef DEBUG
639 printf("R after QR=\n");
640 printf(MDF" "MDF"\n",H[0][0],H[0][1]);
641 for (i=1;i<m-1;i++) {
642 printf(MDF" "MDF" "MDF"\n",H[i][0],H[i][1],H[i][2]);
643 }
644 printf(MDF" "MDF"\n",H[m-1][0],H[m-1][1]);
645 printf("=================\n");
646 printf("Q=\n");
647 for (i=0;i<m;i++) {
648 printf("| ");
649 for (j=0;j<m;j++)
650 printf(""MDF" ",eH(Q,i,j,m));
651 printf("|\n");
652 }
653 #endif
654
655 matrix_mult(H,Q,C,m); /* H=H*Q */
656
657 /* add the shift back to H */
658 /* H = H + mu. I */
659 for (i=1;i<m; i++) {
660 H[i][1]+=mu;
661 }
662 H[0][0]+=mu;
663
664 #ifdef DEBUG
665 printf("RQ=\n");
666 printf(MDF" "MDF"\n",H[0][0],H[0][1]);
667 for (i=1;i<m-1;i++) {
668 printf(MDF" "MDF" "MDF"\n",H[i][0],H[i][1],H[i][2]);
669 }
670 printf(MDF" "MDF"\n",H[m-1][0],H[m-1][1]);
671 printf("=================\n");
672 #endif
673
674 /* examine any off diagonal elements are zero
675 * if nearly zero, make them zero */
676 for (i=1;i<m-1; i++) {
677 if(ABS(H[i][0])<=TOL)
678 H[i][0]=0;
679 if(ABS(H[i][2])<=TOL)
680 H[i][2]=0;
681 }
682 if ( ABS(H[0][1])<TOL) H[0][1]=0;
683 if ( ABS(H[m-1][0])<TOL) H[m-1][0]=0;
684
685 /* now see if any row has both off diagonal elements zero.
686 * if so the diagonal element in that row is an eigenvalue */
687 /* since we are shifting with the least eigenvalue,
688 * the chances are this is the first we find. Hence
689 * check the last row */
690 pivot=0;
691 for (i=1;i<m;i++) {
692 if((H[i][0]==0)) {
693 /* we assume H to be symmetric here
694 * so if lower diagonal has a zero,
695 * upper diagonal also has a zero */
696 pivot=i; /* row */
697 }
698 }
699 /* pivot can be in [1,m-1] only */
700 /* now if pivot != 0, that means we have
701 * found an eigenvalue. Hence we split the matrix H
702 * into two at that point and call find_eigenvalue
703 * recursively on the two sub matrices */
704 if (pivot==1) {
705 /* first element is an eigenvalue */
706 v[0]=H[0][0];
707 #ifdef DEBUG
708 printf("find_eig_tri: splitting lower eigenvalue %d="MDF"\n",pivot,v[0]);
709 #endif
710 /* new sub matrix size m-1 */
711 /* move elements up one row */
712 H[0][0]=H[pivot][1];H[0][1]=H[pivot][2];
713 for (i=1;i<m-1-pivot;i++) {
714 H[i][0]=H[i+pivot][0];
715 H[i][1]=H[i+pivot][1];
716 H[i][2]=H[i+pivot][2];
717 }
718 H[m-1-pivot][0]=H[m-1][0];
719 H[m-1-pivot][1]=H[m-1][1];
720
721 /* eigenvector array */
722 vt=(MY_DOUBLE *) calloc(m-1, sizeof(MY_DOUBLE));
723 if ( vt == 0 ) {
724 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
725 exit(1);
726 }
727 #ifdef DEBUG
728 printf("find_eig_tri:entering recursive more case 1\n");
729 #endif
730 for(i=0;i<m;i++) {
731 free(Q[i]);
732 free(C[i]);
733 }
734 free(Q);free(C);
735
736 /* recursive call */
737 find_eigenvalues_tridiagonal(H,vt,m-1);
738
739 /* copy back result */
740 for(i=1;i<m;i++)
741 v[i]=vt[i-1];
742 #ifdef DEBUG
743 for (i=0;i<m;i++) {
744 printf("find_eig_tri: eigenvalue array after case 1 %d="MDF"\n",i,v[i]);
745 }
746 #endif
747
748 /* free memory */
749 free(vt);
750
751 /* break loop */
752 return(0);
753 } else if (pivot==m-1) {
754 /* last element is an eigenvalue */
755 v[m-1]=H[m-1][1];
756 #ifdef DEBUG
757 printf("find_eig_tri: splitting upper eigenvalue %d="MDF"\n",pivot,v[m-1]);
758 #endif
759 vt=(MY_DOUBLE *) calloc(m-1, sizeof(MY_DOUBLE));
760 if ( vt == 0 ) {
761 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
762 exit(1);
763 }
764 #ifdef DEBUG
765 printf("find_eig_tri:entering recursive more case 2\n");
766 #endif
767 for(i=0;i<m;i++) {
768 free(Q[i]);
769 free(C[i]);
770 }
771 free(Q);free(C);
772
773 /* recursive call */
774 /* now only the first m-1 rows of H is used */
775 find_eigenvalues_tridiagonal(H,vt,m-1);
776
777 /* copy back result */
778 for(i=0;i<m-1;i++)
779 v[i]=vt[i];
780 #ifdef DEBUG
781 for (i=0;i<m;i++) {
782 printf("find_eig_tri: eigenvalue array after case 2 %d="MDF"\n",i,v[i]);
783 }
784 #endif
785
786 free(vt);
787
788 /* break loop */
789 return(0);
790 } else if (pivot) {
791 /* no eigenvalue found, but we
792 * can subdivide the matrix into two */
793 /* first the upper matrix */
794 /* rows 0 to pivot-1, size pivot */
795 /* allocate memory */
796 #ifdef DEBUG
797 printf("find_eig_tri: splitting to two no eigenvalue, pivot at %d\n",pivot);
798 #endif
799 /* new sub matrix size pivot */
800
801 /* eigenvector array */
802 vt=(MY_DOUBLE *) calloc(pivot, sizeof(MY_DOUBLE));
803 if ( vt == 0 ) {
804 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
805 exit(1);
806 }
807 #ifdef DEBUG
808 printf("find_eig_tri:entering recursive more case 3 a\n");
809 #endif
810
811 for(i=0;i<m;i++) {
812 free(Q[i]);
813 free(C[i]);
814 }
815 free(Q);free(C);
816
817 /* recursive call */
818 find_eigenvalues_tridiagonal(H,vt,pivot);
819
820 /* copy back result */
821 for(i=0;i<pivot;i++)
822 v[i]=vt[i];
823 /* free memory */
824 free(vt);
825
826 /* now the lower matrix */
827 /* rows pivot to m-1 */
828 /* new sub matrix size m-pivot */
829 /* move rows upwards by pivot */
830 H[0][0]=H[pivot][1];H[0][1]=H[pivot][2];
831 for (i=1;i<m-1-pivot;i++) {
832 H[i][0]=H[i+pivot][0];
833 H[i][1]=H[i+pivot][1];
834 H[i][2]=H[i+pivot][2];
835 }
836 H[m-1-pivot][0]=H[m-1][0];
837 H[m-1-pivot][1]=H[m-1][1];
838
839
840 /* eigenvector array */
841 vt=(MY_DOUBLE *) calloc(m-pivot, sizeof(MY_DOUBLE));
842 if ( vt == 0 ) {
843 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
844 exit(1);
845 }
846 #ifdef DEBUG
847 printf("find_eig_tri:entering recursive more case 3 b\n");
848 #endif
849
850
851 /* recursive call */
852 find_eigenvalues_tridiagonal(H,vt,m-pivot);
853
854 /* copy back result */
855 for(i=0;i<m-pivot;i++)
856 v[i+pivot]=vt[i];
857 free(vt);
858
859 /* break loop */
860 return(0);
861 }
862 }
863
864 /* if we are here, no matrix subdivision was done
865 * and we exited the loop simply due to
866 * exceeding the number of iterations. need to copy the eigenvalues
867 * we have estimated so far */
868 for (i=0;i<m;i++) {
869 #ifdef DEBUG
870 printf("find_eig_tri: eigenvalues estimated %d="MDF"\n",i,eT(H,i,i,m));
871 #endif
872 v[i]=eT(H,i,i,m);
873 }
874 for(i=0;i<m;i++) {
875 free(Q[i]);
876 free(C[i]);
877 }
878 free(Q);free(C);
879
880
881 return(0);
882
883 }
884
885 static int
find_eigenvalues(MY_DOUBLE ** A,MY_DOUBLE * v,int m)886 find_eigenvalues(MY_DOUBLE **A,MY_DOUBLE *v,int m)
887 {
888 /* finds the eigenvalues of the matrix A and stores them
889 * in array v
890 * A will be modified
891 */
892
893 MY_DOUBLE **H,**Q,**C;
894 int i;
895 #ifdef DEBUG
896 int j;
897 #endif
898 /* setup matrix H 4 diagonal */
899 H= (MY_DOUBLE **) calloc(m, sizeof(MY_DOUBLE *));
900 if ( H == 0 ) {
901 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
902 exit(1);
903 }
904 for (i = 1; i < m-2; ++i) {
905 H[i] = ( MY_DOUBLE * ) calloc( (size_t)4, sizeof(MY_DOUBLE));
906 if ( H[i] == 0 ) {
907 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
908 exit(1);
909 }
910 }
911 H[0] = ( MY_DOUBLE * ) calloc( (size_t)3, sizeof(MY_DOUBLE));
912 if ( H[0] == 0 ) {
913 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
914 exit(1);
915 }
916 H[m-2] = ( MY_DOUBLE * ) calloc( (size_t)3, sizeof(MY_DOUBLE));
917 if ( H[m-2] == 0 ) {
918 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
919 exit(1);
920 }
921 H[m-1] = ( MY_DOUBLE * ) calloc( (size_t)2, sizeof(MY_DOUBLE));
922 if ( H[m-1] == 0 ) {
923 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
924 exit(1);
925 }
926
927 /* Q is an Upper Hessenberg matrix */
928 Q= (MY_DOUBLE **) calloc(m, sizeof(MY_DOUBLE *));
929 if ( Q == 0 ) {
930 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
931 exit(1);
932 }
933 Q[0] = ( MY_DOUBLE * ) calloc( (size_t)m, sizeof(MY_DOUBLE));
934 if ( Q[0] == 0 ) {
935 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
936 exit(1);
937 }
938 for (i=1;i<m;i++) {
939 Q[i] = ( MY_DOUBLE * ) calloc( (size_t)(m-i+1), sizeof(MY_DOUBLE));
940 if ( Q[i] == 0 ) {
941 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
942 exit(1);
943 }
944 }
945 /* C is a temporary matrix, tri diagonal */
946 C= (MY_DOUBLE **) calloc(m, sizeof(MY_DOUBLE *));
947 if ( C == 0 ) {
948 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
949 exit(1);
950 }
951 for (i = 1; i < m-1; ++i) {
952 C[i] = ( MY_DOUBLE * ) calloc( (size_t)3, sizeof(MY_DOUBLE));
953 if ( C[i] == 0 ) {
954 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
955 exit(1);
956 }
957 }
958 C[0] = ( MY_DOUBLE * ) calloc( (size_t)2, sizeof(MY_DOUBLE));
959 if ( C[0] == 0 ) {
960 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
961 exit(1);
962 }
963 C[m-1] = ( MY_DOUBLE * ) calloc( (size_t)2, sizeof(MY_DOUBLE));
964 if ( C[m-1] == 0 ) {
965 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
966 exit(1);
967 }
968
969 convert_to_hessenberg(A,H,m);
970 #ifdef DEBUG
971 printf("A, Hessenberg=\n");
972 for (i=0;i<m;i++) {
973 printf("| ");
974 for (j=0;j<m;j++)
975 printf(""MDF" ",eT(H,i,j,m));
976 printf("|\n");
977 }
978 printf("A =\n");
979 for (i=0;i<m;i++) {
980 printf("| ");
981 for (j=0;j<m;j++)
982 printf(""MDF" ",A[i][j]);
983 printf("|\n");
984 }
985 #endif
986
987 /* call the recursive routine */
988 find_eigenvalues_tridiagonal(H,v,m);
989 #ifdef DEBUG
990 printf("find_eig: eigenvalues found\n");
991 for(i=0;i<m;i++)
992 printf("find_eig: eigenvalue %d="MDF"\n",i,v[i]);
993 #endif
994 for(i=0;i<m;i++) {
995 free(H[i]);
996 free(Q[i]);
997 free(C[i]);
998 }
999 free(H);free(Q);free(C);
1000 return(0);
1001
1002 }
1003
1004
1005 static int
lu_solve(MY_DOUBLE ** A,MY_DOUBLE * x,MY_DOUBLE * y,int m)1006 lu_solve(MY_DOUBLE **A, MY_DOUBLE *x, MY_DOUBLE *y, int m)
1007 {
1008 /* solves Ax=y by LU decomposition with pivoting */
1009 /* A,x, y will be changed */
1010 MY_DOUBLE max;
1011 int i,j,k,*v;
1012 /* decompose A=L. U in place */
1013 /* row vector to store pivot data instead of P */
1014 v =(int*)calloc(m,sizeof(int));
1015 if ( v == 0 ) {
1016 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1017 exit(1);
1018 }
1019 /* U = A, P=I; L=I initially */
1020 for (i=0;i<m;i++) {
1021 v[i]=i; /* initially no permutation */
1022 }
1023
1024 for (k=0; k<m-1;k++) {
1025 /* select i (>=k) with max |U_ik| */
1026 i=k;max=ABS(A[v[i]][k]);
1027 for (j=k;j<m;j++){
1028 if (ABS(A[v[j]][k]) >max){
1029 i=j;max=ABS(A[v[i]][k]);
1030 }
1031 }
1032
1033 /* interchange rows */
1034 /* U_k,k:m-1 <-> U_i,k:m-1 and
1035 * L_k,0:k-1 <-> L_i,0:k-1
1036 */
1037 /* instead we interchange the row index */
1038 j=v[k];
1039 v[k]=v[i];
1040 v[i]=j;
1041
1042 /* now U_k,k should have the max value from column k
1043 * if U_k,k is zero matrix is singular
1044 * hence substitude a small value*/
1045 if (A[v[k]][k]==0){
1046 printf("***: warning: matrix is singular\n");
1047 A[v[k]][k]=0.0000001;
1048 }
1049
1050 for (j=k+1;j<m;j++) {
1051 /* L_j,k=U_j,k/U_k,k */
1052 max=A[v[j]][k]/A[v[k]][k];
1053 /*A[v[j]][k]=A[v[j]][k]/A[v[k]][k]; */
1054 /* U_j,k:m-1=U_j,k:m-1-L_j,k*U_k,k:m-1 */
1055 for (i=k;i<m;i++) {
1056 A[v[j]][i]-=max*A[v[k]][i];
1057 }
1058 A[v[j]][k]=max;
1059 }
1060 }
1061 /* now we have A = L. U */
1062 /* solve system Ax=y */
1063 /* forward elimination */
1064 for (i=0; i<m; i++) {
1065 max=0;
1066 for (j=0;j<i;j++) {
1067 max+=A[v[i]][j]*y[v[j]];
1068 }
1069 y[v[i]]-=max;
1070 }
1071 /* back substitution */
1072 for (i=m-1;i>=0;i--) {
1073 max=0;
1074 for (j=i+1;j<m;j++) {
1075 max+=A[v[i]][j]*x[v[j]];
1076 }
1077 x[v[i]]=(y[v[i]]-max)/A[v[i]][i];
1078 }
1079
1080 /* put solution in correct place */
1081 for (i=0; i<m; i++) {
1082 y[i]=x[v[i]];
1083 }
1084 for (i=0;i<m;i++) {
1085 x[i]=y[i];
1086 }
1087 free(v);
1088 return(0);
1089 }
1090
1091
1092
1093
1094
1095 /* used in sorting eigenvalues, comparison */
1096 static int
compare_doubles(const void * a,const void * b)1097 compare_doubles(const void *a,const void *b)
1098 {
1099 const MY_DOUBLE *da=(const MY_DOUBLE *)a;
1100 const MY_DOUBLE *db=(const MY_DOUBLE *)b;
1101
1102 return((ABS(*da)<=ABS(*db)?1:-1));
1103 }
1104
1105 /* finds all eigenvalues of the matrix a */
1106 static void
find_all_eigenvalues(MY_DOUBLE ** a,MY_DOUBLE * g,MY_INT m)1107 find_all_eigenvalues(MY_DOUBLE **a, MY_DOUBLE *g, MY_INT m)
1108 {
1109 /* finds all eigenvalues of a m by m symmetric matrix a
1110 * a - m by m matrix
1111 * and returns them using an array g - size m by 1
1112 */
1113 MY_DOUBLE **B;
1114
1115 int i;
1116 #ifdef DEBUG
1117 int k,j;
1118 #endif
1119
1120 /* temporary matrix */
1121 B= (MY_DOUBLE **) calloc(m, sizeof(MY_DOUBLE *));
1122 if ( B == 0 ) {
1123 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1124 exit(1);
1125 }
1126 for (i = 0; i < m; ++i) {
1127 B[i] = ( MY_DOUBLE * ) calloc( (size_t)m, sizeof(MY_DOUBLE));
1128 if ( B[i] == 0 ) {
1129 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1130 exit(1);
1131 }
1132 memcpy((void*)B[i],(void*)a[i],(size_t)(m)*sizeof(MY_DOUBLE));
1133 }
1134
1135
1136 /* display */
1137 #ifdef DEBUG
1138 printf("A =\n");
1139 for (i=0;i<m;i++) {
1140 printf("| ");
1141 for (j=0;j<m;j++)
1142 printf(""MDF" ",a[i][j]);
1143 printf("|\n");
1144 }
1145 #endif
1146
1147 /* find eigenvalues */
1148 find_eigenvalues(B,g,m);
1149
1150 /* sort the eigenvalues */
1151 qsort(g,m,sizeof(MY_DOUBLE),compare_doubles);
1152 #ifdef DEBUG
1153 printf("find_eigenvalues_all: eigenvalues found\n");
1154 for (i=0;i<m;i++) {
1155 printf(""MDF" ",g[i]);
1156 }
1157 printf("\n");
1158 #endif
1159
1160
1161 for(i=0;i<m;i++) {
1162 free(B[i]);
1163 }
1164 free(B);
1165
1166 }
1167
1168 /* find the eigenvector(s) given an eigenvalue(s) */
1169 void
find_eigenvectors(MY_DOUBLE ** a,MY_DOUBLE ** x,MY_DOUBLE * ev,MY_INT m,MY_INT n)1170 find_eigenvectors(MY_DOUBLE **a, MY_DOUBLE **x, MY_DOUBLE *ev, MY_INT m, MY_INT n)
1171 {
1172 /* finds the eigenvector of a m by m symmetric matrix a
1173 * a - m by m matrix
1174 * x - eigenvector matrix size (m by n)
1175 * ev - eigenvalue array size n by 1
1176 * n - how many eigenvectors to find starting from the lowest eigenvalue
1177 * m - matrix size of a
1178 */
1179 MY_DOUBLE **B,*y,*v,*vec,norm;
1180
1181 int i,j,en,k;
1182 MY_DOUBLE rayleigh;
1183
1184 /* temporary matrix */
1185 B= (MY_DOUBLE **) calloc(m, sizeof(MY_DOUBLE *));
1186 if ( B == 0 ) {
1187 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1188 exit(1);
1189 }
1190 for (i = 0; i < m; ++i) {
1191 B[i] = ( MY_DOUBLE * ) calloc( (size_t)m, sizeof(MY_DOUBLE));
1192 if ( B[i] == 0 ) {
1193 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1194 exit(1);
1195 }
1196 memcpy((void*)B[i],(void*)a[i],(size_t)(m)*sizeof(MY_DOUBLE));
1197 }
1198
1199 /* allocate memory for array of eigenvalues */
1200 v= (MY_DOUBLE *) calloc(m, sizeof(MY_DOUBLE));
1201 if ( v == 0 ) {
1202 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1203 exit(1);
1204 }
1205
1206
1207 /* display */
1208 #ifdef DEBUG
1209 printf("A =\n");
1210 for (i=0;i<m;i++) {
1211 printf("| ");
1212 for (j=0;j<m;j++)
1213 printf(""MDF" ",a[i][j]);
1214 printf("|\n");
1215 }
1216 #endif
1217
1218 /* find all eigenvalues */
1219 find_all_eigenvalues(a,v,m);
1220
1221 #ifdef DEBUG
1222 printf("eigenvalues found\n");
1223 for (i=0;i<m;i++) {
1224 printf(""MDF" ",v[i]);
1225 }
1226 printf("\n");
1227 #endif
1228
1229 /* do an inverse iteration to find the eigenvector */
1230 y =(MY_DOUBLE*)calloc(m,sizeof(MY_DOUBLE));
1231 if ( y == 0 ) {
1232 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1233 exit(1);
1234 }
1235 vec =(MY_DOUBLE*)calloc(m,sizeof(MY_DOUBLE));
1236 if ( vec == 0 ) {
1237 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1238 exit(1);
1239 }
1240 for (en=m-1;en>=m-n;en--) {
1241 /* initial guess */
1242 memset((void*)vec,0,(size_t)m*sizeof(MY_DOUBLE));
1243 vec[0]=1; /* start with unit vector */
1244
1245 j=0;
1246 while(j<EIGENVECTOR_ITERATION_LIMIT) {
1247 /* inverse iteration */
1248 /* copy A to B */
1249 /* B=A-mu.I */
1250 for (i=0; i<m; i++) {
1251 memcpy((void*)B[i],(void*)a[i],(size_t)(m)*sizeof(MY_DOUBLE));
1252 B[i][i] -=v[en]; /* last eigenvalue */
1253 y[i]=vec[i];
1254 }
1255 #ifdef DEBUG
1256 printf("solving system\n");
1257 for (i=0;i<m;i++) {
1258 printf("| ");
1259 for (k=0;k<m;k++)
1260 printf(""MDF" ",B[i][k]);
1261 printf("|| "MDF"\n",y[i]);
1262 }
1263 #endif
1264
1265 lu_solve(B, vec, y, m);
1266 /* normalize x */
1267 norm=0;
1268 for(i=0;i<m;i++) {
1269 norm+=vec[i]*vec[i];
1270 }
1271 norm=sqrt(norm);
1272 for (i=0; i<m;i++) {
1273 /*vec[i]/=sqrt(norm); */
1274 vec[i]=vec[i]/norm;
1275 }
1276 #ifdef DEBUG
1277 printf("find_eigvec: eigenvector estimate: for eigenvalue "MDF" iteration %d\n",v[en],j);
1278 for (i=0;i<m;i++) {
1279 printf(""MDF" ",vec[i]);
1280 }
1281 printf("\n");
1282 #endif
1283 /* calculate the Rayleigh coefficient to find convergence */
1284 /* store a*v in x */
1285 for(i=0; i<m; i++) {
1286 x[i][m-1-en]=0;
1287 for(k=0; k <m; k++) {
1288 x[i][m-1-en] +=a[i][k]*vec[k];
1289 }
1290 }
1291 /* multiply v^T * a *v */
1292 rayleigh=0;
1293 for(i=0; i<m; i++)
1294 rayleigh+=x[i][m-1-en]*v[i];
1295 if(ABS(rayleigh-v[en]) <= TOL)
1296 j=EIGENVECTOR_ITERATION_LIMIT; /* break loop */
1297 j++;
1298 }
1299
1300 for (i=0;i<m;i++) {
1301 x[i][m-1-en]=vec[i];
1302 #ifdef DEBUG
1303 printf(MDF" ",x[i][m-1-en]);
1304 #endif
1305 }
1306 #ifdef DEBUG
1307 printf("\n");
1308 printf("find_eigvec: eigenvalue %d estimate: "MDF" \n",en, v[en]);
1309 #endif
1310
1311 }
1312 i=0;
1313 for (en=m-1;en>=m-n;en--) {
1314 /* copy eigenvalue to output */
1315 ev[i++]=v[en];
1316 }
1317
1318 for(i=0;i<m;i++) {
1319 free(B[i]);
1320 }
1321 free(B);
1322 free(y);
1323 free(v);
1324 free(vec);
1325 }
1326