(求助C语言实现Householder变换一般实矩阵为上Hessenberg矩阵的算法)

aj6249

<># include "math.h"
/ l+ [1 K+ u7 ]4 }, e1 U#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
      int n;
   double a[],p[],c[],u[];
{   
+ P7 I$ ~, s' u: c' u   int i,j,k,v;
   double sum,asum;
for(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
       for(i=0;i&lt;n;++i)//初始化u[]全为零
; t8 S  _2 D, u2 O  f3 o     u=0.0;</P>
<>
+ s1 X) a' h6 y1 f! d      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
      {
: e& K$ F( D% }  _/ P       k=i*n+j;
' l+ R2 E. ]! `: R* b0 f# J    sum+=a[k]*a[k];
) W  G3 G8 T  M+ M; Z   }
      asum=sqrt(sum);</P>
<>      for(i=j+1;i&lt;n;++i)
   {
6 V7 o9 \1 y1 |* P& ]$ F   
4 c& w' C3 O; S% U0 |2 }    if(i==j+1)
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
% i, u" J- o8 I% i3 O    else
     u=a[i*n+j];
   }
      </P>
3 W4 N  X8 I& T5 f; a, ^: S<>
   sum=0.0;  //实现P
   for(i=0;i&lt;n;++i)
+ R  g# X. O: X1 b6 J( q0 s    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
- d+ ?/ ?. Z4 n+ N+ P   {for(k=0;k&lt;n;++k)
7 d/ {' w8 Z) R4 H& `/ r% _    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
    printf("%13.7e  ",p[i*n+k]);}
         
   printf("\n");}</P>
: [1 p2 \' s% U
3 [$ h$ M# \0 v* j7 z* {* ?<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
        for(v=0;v&lt;n;v++)
5 l- K8 E0 ]( Q/ M) a) k3 V+ V7 Z! K  {  c[i*n+v]=0.0;
   for(k=0;k&lt;n;k++)
0 x/ g0 ?8 w( w- _; C4 l    c[i*n+v]+=p[i*n+k]*a[k*n+v];
) f. j* L+ v& [! h) }+ b5 D  }</P>
<>
     for(i=0;i&lt;n;i++)
' `, j$ l0 A9 S; x! u- @3 |. \! N# J   for(v=0;v&lt;n;v++)
   {
  M% R; ]7 |' U( T. I    a[i*n+v]=0.0;
2 Y. g6 h% K8 U1 a    for(k=0;k&lt;n;k++)
3 F3 `  G+ z! }( `9 m  h     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
6 s/ y" v0 q0 O* M+ k! y2 z<>
}//最开始的for的结束的大括号
return;
/ M  X: p9 F! t3 d# L2 E. f' f  }</P>

" q. _3 M* ^% X# ?2 w/ a5 m4 L$ n9 `<>自己写的运行总是错误</P>

水木年华zzu

我已经把源代码发在计算数学板块了，你自己找下