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

aj6249

<># include "math.h"
( {+ G  p" `( B; m( s8 X3 R4 i#include"stdio.h"</P>
5 O  q0 A: J3 G9 a# \. j<>  void strp(a,p,c, n,u)
& i8 b5 z1 E$ z3 X% S      int n;
   double a[],p[],c[],u[];
4 x' }( n# y- O{   
   int i,j,k,v;
   double sum,asum;
3 j  {/ n$ Q: Mfor(j=0;j&lt;n-2;++j)
( B. K9 J/ h) w  @. n{// 最开始的for 循环
       for(i=0;i&lt;n;++i)//初始化u[]全为零
% k( U2 s; b$ C/ M9 W# ~, |     u=0.0;</P>
7 h/ ?3 m! G! ~<>
* l$ A* |% G1 @" p      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
      {
/ o" `/ A/ Z6 C+ X* o       k=i*n+j;
    sum+=a[k]*a[k];
9 ^) |0 u* Y, Q  }+ T   }
      asum=sqrt(sum);</P>
<>      for(i=j+1;i&lt;n;++i)
   {
+ [+ L( _. O8 ]7 I9 y, Q   
    if(i==j+1)
5 c; @1 ~2 R; j! W  u* d     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
; b4 D0 x$ p9 o7 {6 O    else
' P% K; w% ~! s     u=a[i*n+j];
   }
      </P>
- i% K$ w+ ~1 s<>
7 L' h) h% Z. b* N/ P; k+ q6 V   sum=0.0;  //实现P
   for(i=0;i&lt;n;++i)
0 o. n2 T) k% c5 P. l    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
: u1 _2 I& A$ E/ ~$ p) b" ]) q   {for(k=0;k&lt;n;++k)
7 R! E6 }% ]: F1 i( A# x4 A    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
1 X. F) i/ j$ `( I' f    printf("%13.7e  ",p[i*n+k]);}
         
   printf("\n");}</P>
5 W! s7 E6 t4 `+ ^
) p' N& L+ ]! p# A<>
) [" V+ P( P+ P3 g  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
        for(v=0;v&lt;n;v++)
1 V# {0 |4 v7 ?; d4 H7 d+ B0 v$ F9 |" H  {  c[i*n+v]=0.0;
   for(k=0;k&lt;n;k++)
    c[i*n+v]+=p[i*n+k]*a[k*n+v];
$ j/ ~$ G, x% Q* y4 ~" q  }</P>
<>
0 F* w* U6 N& F& I$ X) C. T     for(i=0;i&lt;n;i++)
. ~! _# Z' d' E4 _, X: p   for(v=0;v&lt;n;v++)
3 F' v/ Y: t# r1 R   {
4 z3 ~+ M6 u6 v+ T" T& H    a[i*n+v]=0.0;
+ N: S6 B* ?* x6 {* ~/ h    for(k=0;k&lt;n;k++)
( z; d( E0 J) r     a[i*n+v]+=c[i*n+k]*p[k*n+v];
9 @: o6 s( w% A$ U% R. X+ ]   }</P>
1 y5 ]0 e4 v7 p& J9 m$ x! Y<>
) g! X; Z# H- R& X" L/ l}//最开始的for的结束的大括号
return;
) M6 t' N/ G3 N: s" z. B  }</P>

  j2 d; S+ P9 I- k  q<>自己写的运行总是错误</P>

水木年华zzu

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