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

aj6249

<># include "math.h"
#include"stdio.h"</P>
7 v3 E1 _7 T$ E* M<>  void strp(a,p,c, n,u)
' v: X6 z6 m- S( L2 {  e" ?( S      int n;
8 u$ ~7 w( E, x  z! E$ @) D/ a( k   double a[],p[],c[],u[];
5 U0 Y* P# Q, O* ?8 [1 d5 S" A, H/ E{   
( @9 Y& w+ [4 ^2 E) x   int i,j,k,v;
$ a. X# W0 N; ^& ^   double sum,asum;
" p0 l: G) F& qfor(j=0;j&lt;n-2;++j)
1 j' h% d: z$ |3 Y( x# S! l1 S. Y{// 最开始的for 循环
8 V$ o  {# r0 \       for(i=0;i&lt;n;++i)//初始化u[]全为零
$ P0 l9 y( r1 ]' q2 s. T     u=0.0;</P>
1 N1 F! X; S1 ]<>
" k2 G/ V) c* m8 M  ~5 r( {) C9 n      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
      {
       k=i*n+j;
0 }8 k0 J8 ]( I4 o. x+ a    sum+=a[k]*a[k];
   }
- Y  W  P5 u& g0 {7 D# z      asum=sqrt(sum);</P>
<>      for(i=j+1;i&lt;n;++i)
; F  P% p% q2 W5 A. K7 Y. k! g   {
   
/ C3 b  l/ `; @6 l* l. L& H    if(i==j+1)
. A. O. z. \/ V+ P     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
    else
     u=a[i*n+j];
# F1 a& }' a5 x. s5 t   }
      </P>
) ~( h+ ]* \1 p; m+ n( z<>
7 d: l+ C' q: ~% m- C   sum=0.0;  //实现P
% u* [$ ?; T3 c5 V; u0 F( j   for(i=0;i&lt;n;++i)
9 O% K& Z3 Y( h+ v8 F+ V    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
1 N2 {% x9 g% m( Q' p   {for(k=0;k&lt;n;++k)
  `$ d# q/ C9 q0 v/ u+ |7 I    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
8 D0 F/ B! V' ^( J    printf("%13.7e  ",p[i*n+k]);}
         
+ j3 u  f2 U% P! d$ K0 h( ~: ?   printf("\n");}</P>
& @/ r' ~% ]* w$ u
  q1 \( D2 A1 w& Z# C<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
        for(v=0;v&lt;n;v++)
1 K1 x4 F2 X1 K* ^; u  {  c[i*n+v]=0.0;
3 R( ?- Y# _3 V   for(k=0;k&lt;n;k++)
3 M- V) B2 \5 s% W9 g% [    c[i*n+v]+=p[i*n+k]*a[k*n+v];
! I) D- F4 D- M6 j' G9 d  }</P>
<>
     for(i=0;i&lt;n;i++)
* N. Z9 b* _8 i   for(v=0;v&lt;n;v++)
   {
    a[i*n+v]=0.0;
    for(k=0;k&lt;n;k++)
# y1 \6 x$ S' _1 I: D     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
' W( ^% S) u: G+ c<>
}//最开始的for的结束的大括号
return;
3 `: y( D1 w6 ^$ ~9 d  }</P>

<>自己写的运行总是错误</P>

水木年华zzu

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