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

aj6249

<># include "math.h"
: S& u. L0 ]2 |+ Z; ?# ^% Y- N#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
      int n;
   double a[],p[],c[],u[];
: u. A  N4 u, A" T{   
   int i,j,k,v;
   double sum,asum;
# Y. W5 I3 Z6 {: p' i+ Gfor(j=0;j&lt;n-2;++j)
2 i' N2 i8 p6 H! S6 x{// 最开始的for 循环
       for(i=0;i&lt;n;++i)//初始化u[]全为零
" s0 t, r- d' o3 ^7 ]$ D) z) s     u=0.0;</P>
<>
      sum=0.0;
. Z  ?4 j% C# k+ m1 V   for(i=j+1;i&lt;n;++i)//实现a
      {
3 P& M8 G! t5 E. l3 j7 Z/ ?3 d       k=i*n+j;
! t, b5 J. D$ E5 r6 m, T+ T' F    sum+=a[k]*a[k];
2 f' }( z$ L+ G   }
      asum=sqrt(sum);</P>
<>      for(i=j+1;i&lt;n;++i)
: H! L6 ^# a+ p- J   {
   
$ j% ^0 Z3 z( ?: @5 ?$ u( Z, h    if(i==j+1)
& Q7 @2 C9 P6 |* B* G1 c% x     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
    else
$ {* r; p/ N% g. N/ y+ J     u=a[i*n+j];
   }
      </P>
<>
   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;
    printf("%13.7e  ",p[i*n+k]);}
, S6 \/ Q. A$ O$ a' e( E% X         
) h( N" ~$ z& {   printf("\n");}</P>

<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
  A6 P4 b2 S5 o5 T' X# ~6 n0 _. O        for(v=0;v&lt;n;v++)
* \4 r: v9 p! J0 W$ Q. h# \+ ^& ^* S1 P  {  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# `/ P5 V5 v<>
; ?, `7 z* ?) J; D( C* u: c     for(i=0;i&lt;n;i++)
& b. s+ s5 r% _   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];
7 E  ~5 c% j$ {/ c. E   }</P>
: w7 w' ~1 o* t+ s' l, ^<>
}//最开始的for的结束的大括号
return;
  }</P>
3 r0 q* ]7 S$ F- h* ?( N( B' v
6 N+ s; H( }5 F- o<>自己写的运行总是错误</P>

水木年华zzu

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