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

aj6249

<># include "math.h"
! a2 @6 Q2 B% Z6 a- I6 w#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
7 r7 M3 g4 O5 x8 ~7 S      int n;
0 w: C2 a1 a! z7 P7 l0 _1 B/ F) t2 n   double a[],p[],c[],u[];
* e* B3 v5 S6 r8 P! t{   
   int i,j,k,v;
" v5 [7 q) w# H   double sum,asum;
3 r* F2 [& s' l; d3 W' y/ Z+ dfor(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
" \( _- o; F$ g8 z6 j1 S' Y, j* c1 p       for(i=0;i&lt;n;++i)//初始化u[]全为零
     u=0.0;</P>
<>
3 a& k$ n! N* @" Z$ A      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
      {
       k=i*n+j;
    sum+=a[k]*a[k];
   }
, ?7 Q) X- @2 r6 V- u: g      asum=sqrt(sum);</P>
* f2 m. d* m+ _- u4 N; h<>      for(i=j+1;i&lt;n;++i)
   {
3 D( _  c& u" W6 l4 p   
/ n3 {8 p2 h; M3 s' ^- D3 K    if(i==j+1)
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
8 {- J7 c1 p! m    else
     u=a[i*n+j];
  w+ G, d* u. ~. e5 o. J   }
      </P>
5 A: v' t" ^) E4 |& i) p* R<>
   sum=0.0;  //实现P
$ U* v2 r( q' Q7 ~( @# t   for(i=0;i&lt;n;++i)
! b  l, l5 j0 l% Z9 X. }    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
   {for(k=0;k&lt;n;++k)
% O! h1 m$ r- f' J) C    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
+ ?1 v% G7 Q* W! n: Y4 _( W    printf("%13.7e  ",p[i*n+k]);}
         
   printf("\n");}</P>
" k8 m( d# [" A5 t* G. A
<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
        for(v=0;v&lt;n;v++)
+ {0 r: [' C% x9 g  {  c[i*n+v]=0.0;
   for(k=0;k&lt;n;k++)
& e: C; g  m$ I5 `, T    c[i*n+v]+=p[i*n+k]*a[k*n+v];
2 W' _* i1 p$ C: O  V  }</P>
<>
     for(i=0;i&lt;n;i++)
8 E( w( _" K; B: ]   for(v=0;v&lt;n;v++)
   {
    a[i*n+v]=0.0;
( h! i0 ]( c! c0 y! N2 X  H    for(k=0;k&lt;n;k++)
     a[i*n+v]+=c[i*n+k]*p[k*n+v];
$ j6 w3 C# n: q   }</P>
/ v6 F8 n4 N3 g<>
) q: q3 N& c/ j2 ?: u3 w; i% o! B}//最开始的for的结束的大括号
return;
  }</P>

: M* v! D3 b) Q- T: r, a<>自己写的运行总是错误</P>

水木年华zzu

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