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

aj6249

<># include "math.h"
! f& ^+ P: `9 q4 X#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
      int n;
   double a[],p[],c[],u[];
9 y; P  {# w% f: e' K{   
   int i,j,k,v;
" @! X5 L; P/ W! P" y   double sum,asum;
3 [2 ~: H+ h4 i2 t; S  `: ]for(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
2 ?7 J7 V3 K& J6 N       for(i=0;i&lt;n;++i)//初始化u[]全为零
     u=0.0;</P>
! q) K# Z! ^3 _8 D/ a$ {<>
      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
) F6 P1 B9 ?. q2 z; t7 f      {
       k=i*n+j;
! d  R( ^, b2 U& \% r) V    sum+=a[k]*a[k];
* @! J1 U* G) A7 {# f& V   }
      asum=sqrt(sum);</P>
/ ]' j) h$ z: ^, Z! \; Q5 z<>      for(i=j+1;i&lt;n;++i)
   {
- @9 M% p$ I8 _) }# f. K   
    if(i==j+1)
: m9 c0 t8 I0 M. s9 U     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
    else
. Q3 g; g* k! N" k     u=a[i*n+j];
# u4 b. O$ G7 E% Q" t+ Z   }
      </P>
<>
   sum=0.0;  //实现P
' N3 n, R5 H$ q. w: {2 C   for(i=0;i&lt;n;++i)
3 S% f% D7 m# {* ~' y0 d8 J    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
$ {, q0 R- |/ S& N   {for(k=0;k&lt;n;++k)
    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
6 z% D4 `6 [# ?    printf("%13.7e  ",p[i*n+k]);}
! `  |" l3 X$ Y; E$ G4 X: I8 A         
   printf("\n");}</P>
  N, d, E5 D/ g: @+ r& m1 R- q% V
5 z4 }" Q( C" ^2 p2 Y. g<>
# J& l6 I" D9 k  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
' E* m8 \: {. D* l        for(v=0;v&lt;n;v++)
  {  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>
& l8 F6 Y3 r, D, c) g<>
     for(i=0;i&lt;n;i++)
4 c- J4 S, I  R# Y3 Y   for(v=0;v&lt;n;v++)
  |; ]9 e' s1 }. [! w1 H   {
2 Q5 {/ s+ \- A/ N$ @* U    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>
4 t, s! c% Y8 j7 n, Y2 i/ A9 ?2 m! D<>
}//最开始的for的结束的大括号
6 ~; {( Z5 O7 q% [9 ? return;
  }</P>

2 j8 z* D3 v; m! Q$ M: W7 _<>自己写的运行总是错误</P>

水木年华zzu

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