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

aj6249

<># include "math.h"
5 d: Q! B! s: S" H#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
      int n;
   double a[],p[],c[],u[];
{   
   int i,j,k,v;
   double sum,asum;
for(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
       for(i=0;i&lt;n;++i)//初始化u[]全为零
( {8 G2 R: I# h4 |     u=0.0;</P>
<>
      sum=0.0;
/ r& Z& ^. m7 H' w0 s' V   for(i=j+1;i&lt;n;++i)//实现a
      {
       k=i*n+j;
: q$ }' ^" w4 }    sum+=a[k]*a[k];
   }
9 Q# J* w' \/ c5 ~; j4 I' t      asum=sqrt(sum);</P>
- O$ }; C7 u8 H3 K<>      for(i=j+1;i&lt;n;++i)
+ r$ l5 p& `: u# [( u2 }7 e   {
   
- w+ B0 b* Z1 m* n5 R% ~! p    if(i==j+1)
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
, Z; R- o$ F, |4 I    else
     u=a[i*n+j];
   }
0 W* p/ Y  d: o      </P>
<>
/ f5 h& R3 X6 R( r4 v/ R   sum=0.0;  //实现P
   for(i=0;i&lt;n;++i)
    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
! N- x% t# ?* B5 D/ I: k. \* c$ o   {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]);}
. P  q" e1 b# G5 g         
$ ]0 W$ j9 @) J   printf("\n");}</P>
6 d; V5 F7 x' |; I  ~  d6 k
<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
6 P5 b7 a/ m7 e' D( z- G        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];
) `8 c: s' ~& Z6 Q: Y5 h7 e  }</P>
& H; f& J9 n7 @<>
     for(i=0;i&lt;n;i++)
   for(v=0;v&lt;n;v++)
   {
& ?  d6 W5 I) l$ g4 W5 V# v1 M1 {    a[i*n+v]=0.0;
! [: a& z7 y7 y' u* y% ~    for(k=0;k&lt;n;k++)
; b5 `4 X  U7 D  \9 h     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
<>
8 `4 X0 |2 W) t0 R: L7 |}//最开始的for的结束的大括号
return;
  }</P>
$ E# @  \1 H+ x" V5 S
4 p, N3 q- U1 A7 ?5 Q. m5 A& s<>自己写的运行总是错误</P>

水木年华zzu

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