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

aj6249

<># include "math.h"
+ ^9 L. H: }: e#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
      int n;
- W6 s4 b. b5 y2 H   double a[],p[],c[],u[];
{   
* Y5 E  E  Z- J9 e   int i,j,k,v;
5 Q/ L+ F6 u% E' }9 B$ @   double sum,asum;
for(j=0;j&lt;n-2;++j)
3 y* K0 }+ X- H8 ^2 t3 Q2 i5 f/ A% m{// 最开始的for 循环
       for(i=0;i&lt;n;++i)//初始化u[]全为零
' W1 L( [- x1 b7 K6 ^     u=0.0;</P>
<>
      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
      {
       k=i*n+j;
    sum+=a[k]*a[k];
- `* e- L5 \4 \, t   }
  {# C$ e  E1 C) p      asum=sqrt(sum);</P>
7 P/ {( k& y- e3 D<>      for(i=j+1;i&lt;n;++i)
   {
   
( K0 r& |7 z# v, U    if(i==j+1)
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
    else
& @% B( x4 F8 r5 D/ [8 K     u=a[i*n+j];
. P. }6 x% W, r8 ^9 H& W   }
# }! [  I$ u6 g" C, `# |      </P>
<>
( o8 M' g( A* V2 k, b   sum=0.0;  //实现P
/ M9 o$ ]* m- K/ v, t1 W7 u   for(i=0;i&lt;n;++i)
    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
( D% i: \5 o; {8 l" |" q   {for(k=0;k&lt;n;++k)
1 F" i) O0 X0 c0 e- \    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
    printf("%13.7e  ",p[i*n+k]);}
+ i  `' H4 j5 I* n  P         
   printf("\n");}</P>

: n$ Q$ p* ?2 I0 z; S* ?<>
0 v0 ]) |6 |% E  h  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
  ~: F; [0 b1 R        for(v=0;v&lt;n;v++)
! z  H' K7 z1 k+ 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];
  }</P>
<>
" r2 k* T# q# G7 f2 y5 i+ u, Q     for(i=0;i&lt;n;i++)
   for(v=0;v&lt;n;v++)
   {
    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];
) b' d6 k& R% I   }</P>
<>
}//最开始的for的结束的大括号
1 a( M" e2 w& }# }! h5 N- J return;
  }</P>
2 K4 U% r; b2 Q2 k4 K+ g
' j1 x9 ^0 ]6 X# T. `<>自己写的运行总是错误</P>

水木年华zzu

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