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

aj6249

<># include "math.h"
- f- P4 K2 A7 R* g( ~& K7 [#include"stdio.h"</P>
3 u; {" U7 |4 `# k7 m( g+ X<>  void strp(a,p,c, n,u)
& O% ?3 f+ F  x      int n;
) I+ H: e6 @2 K% z0 p   double a[],p[],c[],u[];
{   
   int i,j,k,v;
3 _& m9 K- f  [& \/ d; W* S   double sum,asum;
for(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
: x! U' `5 K  K/ T  k       for(i=0;i&lt;n;++i)//初始化u[]全为零
% s4 z0 Q, W& K. S+ K     u=0.0;</P>
<>
      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
9 w3 V$ z: ]3 I# F1 d      {
       k=i*n+j;
0 t( g: q% m7 {. k" U    sum+=a[k]*a[k];
) Z3 u3 b/ g# ?5 Y. B1 A/ G; p   }
4 T; Z( |7 T( V% [      asum=sqrt(sum);</P>
8 k0 e; k' Z1 P4 x/ ^: P2 A  B<>      for(i=j+1;i&lt;n;++i)
   {
   
    if(i==j+1)
8 p% E9 ]  |$ h& s# U# A     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
+ B& h  F8 a) X1 ?& }$ _5 y    else
! X0 ]5 l9 V/ j3 ^( V& w8 V+ K8 O- r. r     u=a[i*n+j];
   }
( w. u1 m! B, j! u1 m0 f# Q      </P>
+ O/ o1 |7 }$ Q4 ?% V3 E' @8 s6 I<>
+ u+ s+ N6 Y) U( G   sum=0.0;  //实现P
* \0 a7 k. E( j9 ]# O   for(i=0;i&lt;n;++i)
* x0 A  z* R: X# }3 T0 Y% M% x    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
   {for(k=0;k&lt;n;++k)
3 p1 v' u8 s- x/ G  [' U" c7 k5 l    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
    printf("%13.7e  ",p[i*n+k]);}
         
   printf("\n");}</P>
  @; e" }$ q- ?- u
<>
# |6 X6 r" r1 Z7 F, D  V5 q+ p3 S  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
        for(v=0;v&lt;n;v++)
9 W# }1 p, k7 K' Y5 w  {  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];
  E% t& N& C" j/ ^  }</P>
<>
     for(i=0;i&lt;n;i++)
   for(v=0;v&lt;n;v++)
5 U! v# ^8 m% L7 g   {
    a[i*n+v]=0.0;
' ?  z# D  [$ t. e; D& a    for(k=0;k&lt;n;k++)
5 R. ~4 L) J( a5 [; I: T     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
<>
2 j1 A8 L% y; Y6 k# n8 A}//最开始的for的结束的大括号
3 N) V3 q# g, ^; L  s return;
  }</P>
9 z* G* I5 _+ D8 B! l# b
<>自己写的运行总是错误</P>

水木年华zzu

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