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

aj6249

<># include "math.h"
#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
      int n;
   double a[],p[],c[],u[];
& R( \- ?- ]% E: R1 Y0 P{   
   int i,j,k,v;
   double sum,asum;
for(j=0;j&lt;n-2;++j)
7 o& F' u9 X1 A  I0 n: B) L{// 最开始的for 循环
' f7 \, z5 h+ |, Z# X       for(i=0;i&lt;n;++i)//初始化u[]全为零
     u=0.0;</P>
9 ^/ {6 h  b+ D. r, I<>
      sum=0.0;
6 s( C- X/ ^7 T   for(i=j+1;i&lt;n;++i)//实现a
5 U& @6 |( j6 Z      {
  ~$ {$ j: r! r; ]0 v       k=i*n+j;
! S3 ]$ h  l+ E2 J    sum+=a[k]*a[k];
   }
      asum=sqrt(sum);</P>
<>      for(i=j+1;i&lt;n;++i)
- H0 n% G: V. l   {
" N: s5 U) s; ?9 X  u3 z5 D' P   
    if(i==j+1)
# z- X, }2 {- p% D* f) P* N     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
8 G% j  I4 s2 ^. g* x+ R    else
3 D+ W3 I* t# j; _1 ?9 N     u=a[i*n+j];
/ X/ E, T+ R* h9 v  |5 X   }
. N: R" [' n& M      </P>
2 t  b6 q0 Q4 x0 A2 N4 p$ y8 \<>
2 O: K  W1 g3 a4 N5 h" B- f   sum=0.0;  //实现P
+ ~2 A3 ?9 G$ [   for(i=0;i&lt;n;++i)
    sum+=(u*u);</P>
4 C$ o' U/ b1 Y! X; s, a' H" P) u<>   for(i=0;i&lt;n;++i)
* {; M9 }; M9 m   {for(k=0;k&lt;n;++k)
7 s8 d6 k4 o2 d# J    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
    printf("%13.7e  ",p[i*n+k]);}
+ t" n6 O( f: c- a2 ?         
  _$ o5 l% v( R; d% F$ V: A   printf("\n");}</P>
, |- Q# j  o! ^% ^. p4 M
0 Z0 S/ b% d4 u4 b/ I<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
9 }1 T0 b( W* E( F3 p        for(v=0;v&lt;n;v++)
# d2 A0 G8 M/ k) m  {  c[i*n+v]=0.0;
' j/ P( M: X6 q+ p& q   for(k=0;k&lt;n;k++)
# I) r# P- s$ S  N9 g" \    c[i*n+v]+=p[i*n+k]*a[k*n+v];
% _: h0 Z7 l  h7 O  }</P>
! t4 y& W( [! t5 d! x! W7 ?: \& L+ V<>
     for(i=0;i&lt;n;i++)
   for(v=0;v&lt;n;v++)
   {
    a[i*n+v]=0.0;
- u! y( G: F3 B7 ~" L    for(k=0;k&lt;n;k++)
( x/ }, O  U3 R5 c  W% F     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
. F1 l) B4 W# w+ U0 B$ B# k0 r0 u<>
}//最开始的for的结束的大括号
" f$ q% X& j1 H& h return;
0 P) ?  R9 M$ m% ^2 V! @5 p  }</P>
% D6 Z+ c7 v4 u+ I/ j) d9 c# P) n( g
  w+ G5 l+ ~  {# j<>自己写的运行总是错误</P>

水木年华zzu

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