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

aj6249

<># include "math.h"
#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
9 n3 l9 @! \9 W, w      int n;
   double a[],p[],c[],u[];
0 u, V5 c( V( |{   
   int i,j,k,v;
   double sum,asum;
for(j=0;j&lt;n-2;++j)
7 d/ ]  y# |8 k+ B& \- b{// 最开始的for 循环
. c8 [) c4 A6 R' [6 \# n# b       for(i=0;i&lt;n;++i)//初始化u[]全为零
$ x$ a" l+ g' v4 O5 O     u=0.0;</P>
<>
" r* u6 R# y+ D. Z8 }6 W      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
* u' D8 i' @! y; `      {
       k=i*n+j;
% c) N  [! ?. }    sum+=a[k]*a[k];
   }
$ F$ T+ J; \# l" F$ J) q$ R* w7 O      asum=sqrt(sum);</P>
<>      for(i=j+1;i&lt;n;++i)
   {
6 v2 e2 s9 T7 _3 X5 V0 R- r   
! C. r6 g% h( Q5 m    if(i==j+1)
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
    else
0 G2 @4 w& E2 L# B: A, s; r     u=a[i*n+j];
   }
2 F$ U5 C5 \) @      </P>
: h9 ?/ @6 Z! W<>
   sum=0.0;  //实现P
/ ?. X) F- `# \5 R  H   for(i=0;i&lt;n;++i)
& I( T5 X4 S# ?1 x+ o. i. [) j    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
   {for(k=0;k&lt;n;++k)
  n: l0 ^% G" y# {  C2 i    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
3 H, W0 p7 C, U( s" y* o4 G5 g    printf("%13.7e  ",p[i*n+k]);}
, Q8 i$ j1 G* w( L6 L8 M8 V         
   printf("\n");}</P>

<>
, k, Q2 g. l7 O+ q" _  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
; B. B% Q% C4 B; u! p+ B        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>
. h: O* f% e9 L  c! a. p2 Z<>
     for(i=0;i&lt;n;i++)
$ A/ ]) }7 K" i$ T   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];
0 L+ k( c2 a+ B   }</P>
6 _# R: U7 m) G- [) c- `<>
: B6 P+ @  ^0 q4 ^6 j}//最开始的for的结束的大括号
return;
  }</P>

<>自己写的运行总是错误</P>

水木年华zzu

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