1 //数值计算程序-特征值和特征向量
2
3 //////////////////////////////////////////////////////////////
4 //约化对称矩阵为三对角对称矩阵
5 //利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵
6 //a-长度为n*n的数组,存放n阶实对称矩阵
7 //n-矩阵的阶数
8 //q-长度为n*n的数组,返回时存放Householder变换矩阵
9 //b-长度为n的数组,返回时存放三对角阵的主对角线元素
10 //c-长度为n的数组,返回时前n-1个元素存放次对角线元素
11 void eastrq(double a[],int n,double q[],double b[],double c[]);
12 //////////////////////////////////////////////////////////////
13 //求实对称三对角对称矩阵的全部特征值及特征向量
14 //利用变型QR方法计算实对称三对角矩阵全部特征值及特征向量
15 //n-矩阵的阶数
16 //b-长度为n的数组,返回时存放三对角阵的主对角线元素
17 //c-长度为n的数组,返回时前n-1个元素存放次对角线元素
18 //q-长度为n*n的数组,若存放单位矩阵,则返回实对称三对角矩阵的特征向量组
19 // 若存放Householder变换矩阵,则返回实对称矩阵A的特征向量组
20 //a-长度为n*n的数组,存放n阶实对称矩阵
21 int ebstq(int n,double b[],double c[],double q[],double eps,int l);
22 //////////////////////////////////////////////////////////////
23 //约化实矩阵为赫申伯格(Hessen berg)矩阵
24 //利用初等相似变换将n阶实矩阵约化为上H矩阵
25 //a-长度为n*n的数组,存放n阶实矩阵,返回时存放上H矩阵
26 //n-矩阵的阶数
27 void echbg(double a[],int n);
28 //////////////////////////////////////////////////////////////
29 //求赫申伯格(Hessen berg)矩阵的全部特征值
30 //返回值小于0表示超过迭代jt次仍未达到精度要求
31 //返回值大于0表示正常返回
32 //利用带原点位移的双重步QR方法求上H矩阵的全部特征值
33 //a-长度为n*n的数组,存放上H矩阵
34 //n-矩阵的阶数
35 //u-长度为n的数组,返回n个特征值的实部
36 //v-长度为n的数组,返回n个特征值的虚部
37 //eps-控制精度要求
38 //jt-整型变量,控制最大迭代次数
39 int edqr(double a[],int n,double u[],double v[],double eps,int jt);
40 //////////////////////////////////////////////////////////////
41 //求实对称矩阵的特征值及特征向量的雅格比法
42 //利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量
43 //返回值小于0表示超过迭代jt次仍未达到精度要求
44 //返回值大于0表示正常返回
45 //a-长度为n*n的数组,存放实对称矩阵,返回时对角线存放n个特征值
46 //n-矩阵的阶数
47 //u-长度为n*n的数组,返回特征向量(按列存储)
48 //eps-控制精度要求
49 //jt-整型变量,控制最大迭代次数
50 int eejcb(double a[],int n,double v[],double eps,int jt);
51 //////////////////////////////////////////////////////////////
52
53
54
55 选自<<徐世良数值计算程序集(C)>>
56
57 每个程序都加上了适当地注释,陆陆续续干了几个月才整理出来的啊。
58
59 今天都给贴出来了
60
61 #include "stdio.h"
62 #include "math.h"
63 //约化对称矩阵为三对角对称矩阵
64 //利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵
65 //a-长度为n*n的数组,存放n阶实对称矩阵
66 //n-矩阵的阶数
67 //q-长度为n*n的数组,返回时存放Householder变换矩阵
68 //b-长度为n的数组,返回时存放三对角阵的主对角线元素
69 //c-长度为n的数组,返回时前n-1个元素存放次对角线元素
70 void eastrq(double a[],int n,double q[],double b[],double c[])
71 {
72 int i,j,k,u,v;
73 double h,f,g,h2;
74 for (i=0; i<=n-1; i++)
75 {
76 for (j=0; j<=n-1; j++)
77 {
78 u=i*n+j; q[u]=a[u];
79 }
80 }
81 for (i=n-1; i>=1; i--)
82 {
83 h=0.0;
84 if (i>1)
85 {
86 for (k=0; k<=i-1; k++)
87 {
88 u=i*n+k;
89 h=h+q[u]*q[u];
90 }
91 }
92
93 if (h+1.0==1.0)
94 {
95 c[i-1]=0.0;
96 if (i==1)
97 {
98 c[i-1]=q[i*n+i-1];
99 }
100 b[i]=0.0;
101 }
102 else
103 {
104 c[i-1]=sqrt(h);
105 u=i*n+i-1;
106 if (q[u]>0.0)
107 {
108 c[i-1]=-c[i-1];
109 }
110 h=h-q[u]*c[i-1];
111 q[u]=q[u]-c[i-1];
112 f=0.0;
113 for (j=0; j<=i-1; j++)
114 {
115 q[j*n+i]=q[i*n+j]/h;
116 g=0.0;
117 for (k=0; k<=j; k++)
118 {
119 g=g+q[j*n+k]*q[i*n+k];
120 }
121 if (j+1<=i-1)
122 {
123 for (k=j+1; k<=i-1; k++)
124 {
125 g=g+q[k*n+j]*q[i*n+k];
126 }
127 }
128 c[j-1]=g/h;
129 f=f+g*q[j*n+i];
130 }
131 h2=f/(h+h);
132 for (j=0; j<=i-1; j++)
133 {
134 f=q[i*n+j];
135 g=c[j-1]-h2*f;
136 c[j-1]=g;
137 for (k=0; k<=j; k++)
138 {
139 u=j*n+k;
140 q[u]=q[u]-f*c[k-1]-g*q[i*n+k];
141 }
142 }
143 b[i]=h;
144 }
145 }
146 b[0]=0.0;
147 for (i=0; i<=n-1; i++)
148 {
149 if ((b[i]!=0.0)&&(i-1>=0))
150 {
151 for (j=0; j<=i-1; j++)
152 {
153 g=0.0;
154 for (k=0; k<=i-1; k++)
155 {
156 g=g+q[i*n+k]*q[k*n+j];
157 }
158 for (k=0; k<=i-1; k++)
159 {
160 u=k*n+j;
161 q[u]=q[u]-g*q[k*n+i];
162 }
163 }
164 }
165 u=i*n+i;
166 b[i]=q[u];
167 q[u]=1.0;
168 if (i-1>=0)
169 {
170 for (j=0; j<=i-1; j++)
171 {
172 q[i*n+j]=0.0;
173 q[j*n+i]=0.0;
174 }
175 }
176 }
177 return;
178
179
180
181
182 //求实对称三对角对称矩阵的全部特征值及特征向量
183 //利用变型QR方法计算实对称三对角矩阵全部特征值及特征向量
184 //n-矩阵的阶数
185 //b-长度为n的数组,返回时存放三对角阵的主对角线元素
186 //c-长度为n的数组,返回时前n-1个元素存放次对角线元素
187 //q-长度为n*n的数组,若存放单位矩阵,则返回实对称三对角矩阵的特征向量组
188 // 若存放Householder变换矩阵,则返回实对称矩阵A的特征向量组
189 //a-长度为n*n的数组,存放n阶实对称矩阵
190 int ebstq(int n,double b[],double c[],double q[],double eps,int l)
191 {
192 int i,j,k,m,it,u,v;
193 double d,f,h,g,p,r,e,s;
194 c[n-1]=0.0;
195 d=0.0;
196 f=0.0;
197 for (j=0; j<=n-1; j++)
198 {
199 it=0;
200 h=eps*(fabs(b[j])+fabs(c[j]));
201 if (h>d)
202 {
203 d=h;
204 }
205 m=j;
206 while ((m<=n-1)&&(fabs(c[m])>d))
207 {
208 m=m+1;
209 }
210 if (m!=j)
211 {
212 do
213 {
214 if (it==l)
215 {
216 printf("fail\n");
217 return(-1);
218 }
219 it=it+1;
220 g=b[j];
221 p=(b[j+1]-g)/(2.0*c[j]);
222 r=sqrt(p*p+1.0);
223 if (p>=0.0)
224 {
225 b[j]=c[j]/(p+r);
226 }
227 else
228 {
229 b[j]=c[j]/(p-r);
230 }
231 h=g-b[j];
232 for (i=j+1; i<=n-1; i++)
233 {
234 b[i]=b[i]-h;
235 }
236 f=f+h;
237 p=b[m];
238 e=1.0;
239 s=0.0;
240 for (i=m-1; i>=j; i--)
241 {
242 g=e*c[i];
243 h=e*p;
244 if (fabs(p)>=fabs(c[i]))
245 {
246 e=c[i]/p;
247 r=sqrt(e*e+1.0);
248 c[i+1]=s*p*r;
249 s=e/r;
250 e=1.0/r;
251 }
252 else
253 {
254 e=p/c[i];
255 r=sqrt(e*e+1.0);
256 c[i+1]=s*c[i]*r;
257 s=1.0/r;
258 e=e/r;
259 }
260 p=e*b[i]-s*g;
261 b[i+1]=h+s*(e*g+s*b[i]);
262 for (k=0; k<=n-1; k++)
263 {
264 u=k*n+i+1;
265 v=u-1;
266 h=q[u];
267 q[u]=s*q[v]+e*h;
268 q[v]=e*q[v]-s*h;
269 }
270 }
271 c[j]=s*p;
272 b[j]=e*p;
273 }
274 while (fabs(c[j])>d);
275 }
276 b[j]=b[j]+f;
277 }
278 for (i=0; i<=n-1; i++)
279 {
280 k=i; p=b[i];
281 if (i+1<=n-1)
282 {
283 j=i+1;
284 while ((j<=n-1)&&(b[j]<=p))
285 {
286 k=j;
287 p=b[j];
288 j=j+1;
289 }
290 }
291 if (k!=i)
292 {
293 b[k]=b[i];
294 b[i]=p;
295 for (j=0; j<=n-1; j++)
296 {
297 u=j*n+i;
298 v=j*n+k;
299 p=q[u];
300 q[u]=q[v];
301 q[v]=p;
302 }
303 }
304 }
305 return(1);
306 }
307
308 //约化实矩阵为赫申伯格(Hessen berg)矩阵
309 //利用初等相似变换将n阶实矩阵约化为上H矩阵
310 //a-长度为n*n的数组,存放n阶实矩阵,返回时存放上H矩阵
311 //n-矩阵的阶数
312 void echbg(double a[],int n)
313 { int i,j,k,u,v;
314 double d,t;
315 for (k=1; k<=n-2; k++)
316 {
317 d=0.0;
318 for (j=k; j<=n-1; j++)
319 {
320 u=j*n+k-1;
321 t=a[u];
322 if (fabs(t)>fabs(d))
323 {
324 d=t;
325 i=j;
326 }
327 }
328 if (fabs(d)+1.0!=1.0)
329 {
330 if (i!=k)
331 {
332 for (j=k-1; j<=n-1; j++)
333 {
334 u=i*n+j;
335 v=k*n+j;
336 t=a[u];
337 a[u]=a[v];
338 a[v]=t;
339 }
340 for (j=0; j<=n-1; j++)
341 {
342 u=j*n+i;
343 v=j*n+k;
344 t=a[u];
345 a[u]=a[v];
346 a[v]=t;
347 }
348 }
349 for (i=k+1; i<=n-1; i++)
350 {
351 u=i*n+k-1;
352 t=a[u]/d;
353 a[u]=0.0;
354 for (j=k; j<=n-1; j++)
355 {
356 v=i*n+j;
357 a[v]=a[v]-t*a[k*n+j];
358 }
359 for (j=0; j<=n-1; j++)
360 {
361 v=j*n+k;
362 a[v]=a[v]+t*a[j*n+i];
363 }
364 }
365 }
366 }
367 return;
368 }
369
370
371
372 //求赫申伯格(Hessen berg)矩阵的全部特征值
373 //利用带原点位移的双重步QR方法求上H矩阵的全部特征值
374 //返回值小于0表示超过迭代jt次仍未达到精度要求
375 //返回值大于0表示正常返回
376 //a-长度为n*n的数组,存放上H矩阵
377 //n-矩阵的阶数
378 //u-长度为n的数组,返回n个特征值的实部
379 //v-长度为n的数组,返回n个特征值的虚部
380 //eps-控制精度要求
381 //jt-整型变量,控制最大迭代次数
382 int edqr(double a[],int n,double u[],double v[],double eps,int jt)
383 {
384 int m,it,i,j,k,l,ii,jj,kk,ll;
385 double b,c,w,g,xy,p,q,r,x,s,e,f,z,y;
386 it=0;
387 m=n;
388 while (m!=0)
389 {
390 l=m-1;
391 while ((l>0)&&(fabs(a[l*n+l-1])>eps*(fabs(a[(l-1)*n+l-1])+fabs(a[l*n+l]))))
392 {
393 l=l-1;
394 }
395 ii=(m-1)*n+m-1;
396 jj=(m-1)*n+m-2;
397 kk=(m-2)*n+m-1;
398 ll=(m-2)*n+m-2;
399 if (l==m-1)
400 {
401 u[m-1]=a[(m-1)*n+m-1];
402 v[m-1]=0.0;
403 m=m-1; it=0;
404 }
405 else if (l==m-2)
406 {
407 b=-(a[ii]+a[ll]);
408 c=a[ii]*a[ll]-a[jj]*a[kk];
409 w=b*b-4.0*c;
410 y=sqrt(fabs(w));
411 if (w>0.0)
412 {
413 xy=1.0;
414 if (b<0.0)
415 {
416 xy=-1.0;
417 }
418 u[m-1]=(-b-xy*y)/2.0;
419 u[m-2]=c/u[m-1];
420 v[m-1]=0.0; v[m-2]=0.0;
421 }
422 else
423 {
424 u[m-1]=-b/2.0;
425 u[m-2]=u[m-1];
426 v[m-1]=y/2.0;
427 v[m-2]=-v[m-1];
428 }
429 m=m-2;
430 it=0;
431 }
432 else
433 {
434 if (it>=jt)
435 {
436 printf("fail\n");
437 return(-1);
438 }
439 it=it+1;
440 for (j=l+2; j<=m-1; j++)
441 {
442 a[j*n+j-2]=0.0;
443 }
444 for (j=l+3; j<=m-1; j++)
445 {
446 a[j*n+j-3]=0.0;
447 }
448 for (k=l; k<=m-2; k++)
449 {
450 if (k!=l)
451 {
452 p=a[k*n+k-1];
453 q=a[(k+1)*n+k-1];
454 r=0.0;
455 if (k!=m-2)
456 {
457 r=a[(k+2)*n+k-1];
458 }
459 }
460 else
461 {
462 x=a[ii]+a[ll];
463 y=a[ll]*a[ii]-a[kk]*a[jj];
464 ii=l*n+l;
465 jj=l*n+l+1;
466 kk=(l+1)*n+l;
467 ll=(l+1)*n+l+1;
468 p=a[ii]*(a[ii]-x)+a[jj]*a[kk]+y;
469 q=a[kk]*(a[ii]+a[ll]-x);
470 r=a[kk]*a[(l+2)*n+l+1];
471 }
472 if ((fabs(p)+fabs(q)+fabs(r))!=0.0)
473 {
474 xy=1.0;
475 if (p<0.0)
476 {
477 xy=-1.0;
478 }
479 s=xy*sqrt(p*p+q*q+r*r);
480 if (k!=l)
481 {
482 a[k*n+k-1]=-s;
483 }
484 e=-q/s;
485 f=-r/s;
486 x=-p/s;
487 y=-x-f*r/(p+s);
488 g=e*r/(p+s);
489 z=-x-e*q/(p+s);
490 for (j=k; j<=m-1; j++)
491 {
492 ii=k*n+j;
493 jj=(k+1)*n+j;
494 p=x*a[ii]+e*a[jj];
495 q=e*a[ii]+y*a[jj];
496 r=f*a[ii]+g*a[jj];
497 if (k!=m-2)
498 {
499 kk=(k+2)*n+j;
500 p=p+f*a[kk];
501 q=q+g*a[kk];
502 r=r+z*a[kk];
503 a[kk]=r;
504 }
505 a[jj]=q;
506 a[ii]=p;
507 }
508 j=k+3;
509 if (j>=m-1)
510 {
511 j=m-1;
512 }
513 for (i=l; i<=j; i++)
514 {
515 ii=i*n+k;
516 jj=i*n+k+1;
517 p=x*a[ii]+e*a[jj];
518 q=e*a[ii]+y*a[jj];
519 r=f*a[ii]+g*a[jj];
520 if (k!=m-2)
521 {
522 kk=i*n+k+2;
523 p=p+f*a[kk];
524 q=q+g*a[kk];
525 r=r+z*a[kk];
526 a[kk]=r;
527 }
528 a[jj]=q;
529 a[ii]=p;
530 }
531 }
532 }
533 }
534 }
535 return(1);
536 }
537
538
539
540 //求实对称矩阵的特征值及特征向量的雅格比法
541 //利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量
542 //返回值小于0表示超过迭代jt次仍未达到精度要求
543 //返回值大于0表示正常返回
544 //a-长度为n*n的数组,存放实对称矩阵,返回时对角线存放n个特征值
545 //n-矩阵的阶数
546 //u-长度为n*n的数组,返回特征向量(按列存储)
547 //eps-控制精度要求
548 //jt-整型变量,控制最大迭代次数
549 int eejcb(double a[],int n,double v[],double eps,int jt)
550 {
551 int i,j,p,q,u,w,t,s,l;
552 double fm,cn,sn,omega,x,y,d;
553 l=1;
554 for (i=0; i<=n-1; i++)
555 {
556 v[i*n+i]=1.0;
557 for (j=0; j<=n-1; j++)
558 {
559 if (i!=j)
560 {
561 v[i*n+j]=0.0;
562 }
563 }
564 }
565 while (1==1)
566 {
567 fm=0.0;
568 for (i=0; i<=n-1; i++)
569 {
570 for (j=0; j<=n-1; j++)
571 {
572 d=fabs(a[i*n+j]);
573 if ((i!=j)&&(d>fm))
574 {
575 fm=d;
576 p=i;
577 q=j;
578 }
579 }
580 }
581 if (fm<eps)
582 {
583 return(1);
584 }
585 if (l>jt)
586 {
587 return(-1);
588 }
589 l=l+1;
590 u=p*n+q;
591 w=p*n+p;
592 t=q*n+p;
593 s=q*n+q;
594 x=-a[u];
595 y=(a[s]-a[w])/2.0;
596 omega=x/sqrt(x*x+y*y);
597 if (y<0.0)
598 {
599 omega=-omega;
600 }
601 sn=1.0+sqrt(1.0-omega*omega);
602 sn=omega/sqrt(2.0*sn);
603 cn=sqrt(1.0-sn*sn);
604 fm=a[w];
605 a[w]=fm*cn*cn+a[s]*sn*sn+a[u]*omega;
606 a[s]=fm*sn*sn+a[s]*cn*cn-a[u]*omega;
607 a[u]=0.0;
608 a[t]=0.0;
609 for (j=0; j<=n-1; j++)
610 {
611 if ((j!=p)&&(j!=q))
612 {
613 u=p*n+j;
614 w=q*n+j;
615 fm=a[u];
616 a[u]=fm*cn+a[w]*sn;
617 a[w]=-fm*sn+a[w]*cn;
618 }
619 }
620 for (i=0; i<=n-1; i++)
621 {
622 if ((i!=p)&&(i!=q))
623 {
624 u=i*n+p;
625 w=i*n+q;
626 fm=a[u];
627 a[u]=fm*cn+a[w]*sn;
628 a[w]=-fm*sn+a[w]*cn;
629 }
630 }
631 for (i=0; i<=n-1; i++)
632 {
633 u=i*n+p;
634 w=i*n+q;
635 fm=v[u];
636 v[u]=fm*cn+v[w]*sn;
637 v[w]=-fm*sn+v[w]*cn;
638 }
639 }
640 return(1);
641 }