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

aj6249

<># include "math.h"
#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
      int n;
& ]) s; C( `- s* P7 P3 T! r   double a[],p[],c[],u[];
4 L& {4 N( v; A' L/ y% o9 O{   
5 y8 p$ J7 V+ z* a$ n" B* ]: G5 V   int i,j,k,v;
6 h& A' M' B3 A" W& u8 s* L" D' z   double sum,asum;
for(j=0;j&lt;n-2;++j)
' E0 T5 q1 ~- H! Z{// 最开始的for 循环
       for(i=0;i&lt;n;++i)//初始化u[]全为零
7 A. i* d. t/ a3 p8 Z- F     u=0.0;</P>
<>
' L( Y# \' v+ ~; S& T7 i5 |      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
      {
       k=i*n+j;
& G& j* ~* ~9 j$ E/ O& S7 L1 U2 @    sum+=a[k]*a[k];
   }
      asum=sqrt(sum);</P>
' y5 |" J* T2 c: `( {' N<>      for(i=j+1;i&lt;n;++i)
   {
   
    if(i==j+1)
+ R+ B: A  s. F3 F5 N8 g+ ?     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
4 X& H" O& r5 q! J5 @$ p& i: N( _    else
     u=a[i*n+j];
  o$ `* K) [, ^, u" u1 z   }
      </P>
<>
6 \( G( E+ T" k   sum=0.0;  //实现P
   for(i=0;i&lt;n;++i)
    sum+=(u*u);</P>
<>   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;
5 `: R; J% B- r1 i; T; R    printf("%13.7e  ",p[i*n+k]);}
         
   printf("\n");}</P>

<>
! ~6 R9 u( D" s  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
: f# o0 e6 X+ s" `& F9 T+ N: G        for(v=0;v&lt;n;v++)
  {  c[i*n+v]=0.0;
+ \; Q# [2 n1 I2 I5 d  V) n   for(k=0;k&lt;n;k++)
2 S6 O: [; l) I9 x3 o7 |4 g    c[i*n+v]+=p[i*n+k]*a[k*n+v];
  }</P>
& X8 d( D2 U# J7 X3 Q- A9 s7 A: ~<>
# n2 f* m+ o3 f# Z' v9 k$ d+ E     for(i=0;i&lt;n;i++)
8 I# n2 y1 _8 F$ V* T1 X! n   for(v=0;v&lt;n;v++)
! T/ s) C9 I# L+ L   {
8 C" m3 Z- @7 J2 c0 ^% G    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];
5 Q  ~* K( s! T) K9 Z   }</P>
' N) q1 l* T0 s3 W+ q, @3 g<>
}//最开始的for的结束的大括号
+ ~' w: }! [9 T9 U5 d return;
  }</P>
0 \' T# V3 |, Q. w
$ p8 T4 i' f5 |6 K8 @# I<>自己写的运行总是错误</P>

水木年华zzu

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