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

aj6249

<># include "math.h"
8 |1 ^: F1 F- J- D: E0 j#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
      int n;
- w) P3 E+ J& s# x+ E   double a[],p[],c[],u[];
- J1 T4 t7 S% B# l7 j/ l{   
   int i,j,k,v;
. \" J+ F" F) Q/ @9 `   double sum,asum;
. H7 B$ g+ a- F0 Ffor(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
' a9 O" v& l* T( a       for(i=0;i&lt;n;++i)//初始化u[]全为零
* A% \2 Y* i# {( @8 x, I  d     u=0.0;</P>
; x. J/ a4 i  G; W<>
6 F$ k: h% ^' U/ O      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
- i" U  P" W" Q- y      {
3 z# l& e# d' c2 l6 t4 E( ]% v       k=i*n+j;
% X; m: g! x' z+ J! `. |    sum+=a[k]*a[k];
; A/ y$ T6 z' r5 C0 ^   }
! U# Z5 `; ~0 `/ {: @0 ]9 a& h# v      asum=sqrt(sum);</P>
<>      for(i=j+1;i&lt;n;++i)
3 T3 O* M+ X) o8 Y1 S" J' k   {
   
9 h5 @4 f, d0 A$ b) H2 O    if(i==j+1)
2 [7 f( S2 d: z/ }5 M9 |7 g; S     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
8 d9 @* \, t$ V    else
     u=a[i*n+j];
1 X) K: w; G4 n+ b1 P/ O7 w/ W   }
      </P>
<>
   sum=0.0;  //实现P
   for(i=0;i&lt;n;++i)
! N: D6 y" a1 K$ c9 s$ T    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
& {5 k8 G2 d: r( ?+ j1 @  z   {for(k=0;k&lt;n;++k)
    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
# P/ t4 ]) S0 }    printf("%13.7e  ",p[i*n+k]);}
         
   printf("\n");}</P>
8 S9 [$ I6 Y: Y7 N
<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
        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];
" F" Y$ U/ L( d/ @  }</P>
<>
     for(i=0;i&lt;n;i++)
  N5 d7 E' G( P& s+ m* ^   for(v=0;v&lt;n;v++)
   {
    a[i*n+v]=0.0;
5 [  P) ~, s5 K7 Z2 @8 ]+ J    for(k=0;k&lt;n;k++)
# X1 k- D; [6 P+ q' U     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
<>
}//最开始的for的结束的大括号
- u$ X2 m& b8 E+ T5 \, v return;
3 Z7 g1 K$ ?; `) h7 F7 h4 U  }</P>

+ ^% Q( Z* S! V) ]- k<>自己写的运行总是错误</P>

水木年华zzu

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