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 "); 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 "); 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 }
转自：http://www.cnblogs.com/longdouhzt/archive/2012/04/22/2465191.html