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

aj6249

<># include "math.h"
#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
! B; H5 l5 I% Y0 P- p  @' c8 {2 V      int n;
( W) J- r5 N1 }* M+ p   double a[],p[],c[],u[];
. _! l5 k0 A" j* b( V! p{   
% M: }0 y9 e0 B' g* N% {! c: ~   int i,j,k,v;
- F  d$ h% g  m8 g   double sum,asum;
for(j=0;j&lt;n-2;++j)
0 ^7 L* P4 i* m2 \6 S" g% T: K) b{// 最开始的for 循环
       for(i=0;i&lt;n;++i)//初始化u[]全为零
) [, t# Y1 k3 f6 r% T' E( c     u=0.0;</P>
<>
      sum=0.0;
- {+ r5 V; W7 O   for(i=j+1;i&lt;n;++i)//实现a
      {
       k=i*n+j;
' m- [1 t8 u4 E; g5 l    sum+=a[k]*a[k];
: e: a/ ?0 l: ~+ q: v/ x5 v   }
5 y: r5 |- N$ i4 P      asum=sqrt(sum);</P>
3 }1 w$ u9 t9 U! n<>      for(i=j+1;i&lt;n;++i)
3 m* @7 [* X1 h9 l8 w   {
( D9 o# [# O8 v( B   
% v9 _4 K' e+ n& U4 M- }    if(i==j+1)
( q1 z# Y$ r5 x- d     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
    else
6 r0 B2 K' i) J- w  `     u=a[i*n+j];
7 ~( Y& G1 j# C+ F5 r   }
      </P>
<>
   sum=0.0;  //实现P
) W0 S- \  J# R. h7 b   for(i=0;i&lt;n;++i)
& r3 ~3 J- Y$ F9 S0 g: N5 M9 L    sum+=(u*u);</P>
" f* G) [5 K3 s4 {<>   for(i=0;i&lt;n;++i)
   {for(k=0;k&lt;n;++k)
    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
8 n$ p7 V" J% x9 G- m: w    printf("%13.7e  ",p[i*n+k]);}
         
   printf("\n");}</P>

8 f2 ~  L0 c( o; q9 q/ P<>
& J  T9 |! b0 a  d( o7 G( _- I  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
: O. b, r4 y6 X  l4 a' T6 w        for(v=0;v&lt;n;v++)
" ]5 e5 o3 a, Y; F+ p  {  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];
  }</P>
4 a# D4 ?$ C" ?0 o# A4 V* A+ i* L<>
7 F$ _) q8 r4 O% C     for(i=0;i&lt;n;i++)
7 B; C  L$ ^2 U6 p( Z7 N+ [1 p+ T   for(v=0;v&lt;n;v++)
   {
7 X, \0 b' x6 m    a[i*n+v]=0.0;
    for(k=0;k&lt;n;k++)
     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
<>
}//最开始的for的结束的大括号
return;
  }</P>
1 |2 {0 Q; l  b7 l4 _; N( N
<>自己写的运行总是错误</P>

水木年华zzu

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