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

aj6249

<># include "math.h"
, F$ O2 h: o/ V# ~2 G; J0 C  H4 N#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
      int n;
   double a[],p[],c[],u[];
{   
5 ]' Y; ^; S" D) ]* U# }9 `   int i,j,k,v;
- `2 p: ~1 `; l: W   double sum,asum;
6 `3 [" z, ~3 xfor(j=0;j&lt;n-2;++j)
5 B: J% O1 `" Q{// 最开始的for 循环
       for(i=0;i&lt;n;++i)//初始化u[]全为零
     u=0.0;</P>
0 I( u1 [9 {1 X' J/ ]7 h0 [4 b9 e<>
' B- ^4 F- _) x+ }/ |8 ^      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
      {
       k=i*n+j;
% o7 j1 O( t1 i& O' b: K4 u$ p; Y    sum+=a[k]*a[k];
   }
      asum=sqrt(sum);</P>
! G# b7 R" J0 l# W<>      for(i=j+1;i&lt;n;++i)
% B& S1 F0 X0 h3 J8 F   {
% S; U" U( ~% ]   
    if(i==j+1)
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
    else
1 t* J0 x1 `& E- X* ]" Y     u=a[i*n+j];
4 N. @9 i6 `7 a3 Q3 [8 q   }
      </P>
# `$ j3 X0 x0 E, W8 y, `! W; l2 X8 ^<>
   sum=0.0;  //实现P
8 D9 t9 \" d# }   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)
7 c2 T6 F& P+ P; Q4 C3 t! r$ p    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
    printf("%13.7e  ",p[i*n+k]);}
         
   printf("\n");}</P>
$ q" r2 G3 i8 t, A* Q0 x& _
* a3 j& {/ I& H<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
* |* i6 ]4 P. \' m) Z& J! }        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>
! _( W. [$ |& K5 l+ \2 J- H<>
     for(i=0;i&lt;n;i++)
   for(v=0;v&lt;n;v++)
   {
    a[i*n+v]=0.0;
" O- w' R2 I0 W+ w  o) Q( X    for(k=0;k&lt;n;k++)
& F% X. [& q# r( T! V4 F5 X6 J' \     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
# C5 M. R; O4 V7 F) U- p7 j# Q<>
}//最开始的for的结束的大括号
3 [# |6 v8 P7 M$ h8 x# z return;
  }</P>

7 f+ f4 T; \0 ?& q; k<>自己写的运行总是错误</P>

水木年华zzu

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