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

aj6249

<># include "math.h"
+ ?+ q7 n1 G& E9 ?2 Z. h* b#include"stdio.h"</P>
4 p# p8 y3 h8 o9 J! b0 U# z<>  void strp(a,p,c, n,u)
6 `; U( C) R: V$ x) Y" |5 u      int n;
4 Z# o* p3 i9 G2 n! h$ Y   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[]全为零
  W. {7 p. U$ s1 s) n) h     u=0.0;</P>
) {+ d6 ]# E7 d8 {<>
4 h* n( Z! E3 L6 Q      sum=0.0;
. ~8 d7 _3 H7 d4 i6 i" t* I2 d   for(i=j+1;i&lt;n;++i)//实现a
      {
       k=i*n+j;
    sum+=a[k]*a[k];
   }
" d( t' I  T) T      asum=sqrt(sum);</P>
<>      for(i=j+1;i&lt;n;++i)
   {
6 `, z0 H. [4 v! K# r5 k* x$ k. q   
    if(i==j+1)
- F2 B" z, x1 i7 f' t. m  O: [     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
$ l* @1 X" E& u, y5 r7 K0 h4 k: X    else
3 A' k0 h2 Z2 G+ K2 z# w     u=a[i*n+j];
   }
' F/ \4 B- Z) @% n/ }' R% u. t  {% ~9 V, G      </P>
<>
   sum=0.0;  //实现P
( y; p' a* ]% o3 z# U   for(i=0;i&lt;n;++i)
9 Y0 L0 M9 W1 s# @    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
0 {( i& x( X) P   {for(k=0;k&lt;n;++k)
) w7 O3 H, E  O0 W0 x    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
    printf("%13.7e  ",p[i*n+k]);}
" E' Y. ], s" N7 r, b) \$ h         
   printf("\n");}</P>
& o9 F, m8 x! ~2 U
<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
$ j# c. F1 n% v, V8 N        for(v=0;v&lt;n;v++)
4 W- [' Q) \+ k1 P/ |2 l5 _  {  c[i*n+v]=0.0;
   for(k=0;k&lt;n;k++)
6 S8 p; ]  d! ?) i5 C    c[i*n+v]+=p[i*n+k]*a[k*n+v];
  }</P>
$ i: G6 a8 h- E7 q8 \- E6 g<>
     for(i=0;i&lt;n;i++)
   for(v=0;v&lt;n;v++)
   {
    a[i*n+v]=0.0;
/ E4 H* X# n, R+ E( k( Q. v    for(k=0;k&lt;n;k++)
! V. a1 e9 S9 [6 X3 |     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
3 u6 {, S) `! H6 l+ h3 G<>
}//最开始的for的结束的大括号
return;
( s6 n' @- Y+ o' M$ k5 _& ?  }</P>

: P, Y: O* X: A<>自己写的运行总是错误</P>

水木年华zzu

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