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

aj6249

<># include "math.h"
#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
* S3 o8 n% f7 ]3 q' \  Q) W0 ?1 K      int n;
# o0 U/ r* ]: y   double a[],p[],c[],u[];
{   
   int i,j,k,v;
+ K: U2 }  @& y( i/ V) m$ Y   double sum,asum;
- @# Q2 T9 ~9 l$ A; B8 |for(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
2 p% \. `3 L% v* y5 x! }8 {6 A3 W       for(i=0;i&lt;n;++i)//初始化u[]全为零
8 D0 F; N& S5 C; X% `; E, V     u=0.0;</P>
" d1 j; c" d7 m$ g/ ?<>
      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
0 Z1 B' _+ l5 [  L      {
) R: S2 v- c. g/ L       k=i*n+j;
    sum+=a[k]*a[k];
   }
. }- C5 W) H" C, t' P      asum=sqrt(sum);</P>
$ D! y/ d3 d) A/ @* s& S; g  c2 j3 R) |7 ?# n<>      for(i=j+1;i&lt;n;++i)
! m3 v" b8 \2 `. x) d+ L   {
   
6 V. ^: c1 ^! }% _& ?- Q1 o    if(i==j+1)
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
! l9 ^/ b" X. B/ u8 g4 |/ B    else
     u=a[i*n+j];
   }
1 c1 S. C- Z( a2 I9 W/ E      </P>
' O5 j( y6 O% [4 M% B. ^<>
3 d$ J( _2 V5 f. ?   sum=0.0;  //实现P
& x. G2 S' @/ s" y6 w) K& U   for(i=0;i&lt;n;++i)
# q$ S5 D" S2 }2 ?) u    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
6 h9 Q  t- r, b, w! l) a2 N0 k" H" c   {for(k=0;k&lt;n;++k)
    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
8 |9 K7 U- l# m# Z- T( _3 w    printf("%13.7e  ",p[i*n+k]);}
- h) l# `$ P  z+ h5 }& o         
   printf("\n");}</P>

1 ]8 D5 x- s5 l  q<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
        for(v=0;v&lt;n;v++)
# q$ X$ c3 G. u$ c' M" J- E  {  c[i*n+v]=0.0;
# M6 I  E3 O. {, U   for(k=0;k&lt;n;k++)
+ R. x7 g7 f( c) @( \    c[i*n+v]+=p[i*n+k]*a[k*n+v];
  }</P>
<>
) U) C# \. L5 Q7 ~) J# K     for(i=0;i&lt;n;i++)
   for(v=0;v&lt;n;v++)
4 x: h7 i) v, n: }- c   {
    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];
! \0 K3 e- t. A   }</P>
<>
7 w- a& D2 j; s+ R" s  v}//最开始的for的结束的大括号
3 M! o0 X! h0 R5 M1 ]. ]' ? return;
" t2 ]& x4 f/ v1 }9 |  }</P>
  t$ h# |: X& ~$ T1 h2 m/ X
<>自己写的运行总是错误</P>

水木年华zzu

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