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

aj6249

<># include "math.h"
  @4 E* k( j% w" W( o- x#include"stdio.h"</P>
2 D% \, z5 _$ G; y5 b0 s* t8 H# B  ^<>  void strp(a,p,c, n,u)
% X# I6 z2 G  ]! N, D$ T' q      int n;
   double a[],p[],c[],u[];
; O' N- S& P0 c, I{   
) _! L" y$ ]/ d9 b) f; c1 W   int i,j,k,v;
2 g3 }/ ^" E3 R   double sum,asum;
for(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
       for(i=0;i&lt;n;++i)//初始化u[]全为零
' ]: j0 p6 h7 l     u=0.0;</P>
<>
      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
' A3 {* F4 u) \5 S      {
       k=i*n+j;
    sum+=a[k]*a[k];
9 d& I" T3 w' a5 W' E' \   }
      asum=sqrt(sum);</P>
/ n* x  U* v2 }( C' O<>      for(i=j+1;i&lt;n;++i)
( T9 _9 {# t  m  P. a' }   {
   
    if(i==j+1)
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
" U# D' J1 M3 N- N$ G    else
     u=a[i*n+j];
   }
      </P>
<>
' C  ~* E2 ^- O) @- V   sum=0.0;  //实现P
4 v1 v4 F7 O# N- p) o   for(i=0;i&lt;n;++i)
    sum+=(u*u);</P>
+ h9 @. P1 ^6 z3 t0 B* C4 x<>   for(i=0;i&lt;n;++i)
# }4 V; K8 Z  ~; r# \   {for(k=0;k&lt;n;++k)
    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
& }# w( J* A( G& B, M- k1 M    printf("%13.7e  ",p[i*n+k]);}
         
. R/ \  E/ D6 l* D6 ^; N   printf("\n");}</P>

<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
0 d; V. V8 ]/ c& p# ?* V9 I$ M        for(v=0;v&lt;n;v++)
  {  c[i*n+v]=0.0;
! D4 G3 A/ w) {+ t' x   for(k=0;k&lt;n;k++)
7 g. A- H# X6 z* N. M+ `5 {1 s    c[i*n+v]+=p[i*n+k]*a[k*n+v];
  F  b2 W4 P0 L3 E) S; S/ E, |1 d  }</P>
0 |* q. V- F$ f$ _9 H& D; ?<>
& \* g1 L1 L$ _' f0 u- b. s7 L' _     for(i=0;i&lt;n;i++)
   for(v=0;v&lt;n;v++)
( R; w" Z" A7 e2 W" W8 R) Z   {
6 z7 u( e1 O  q3 Z6 W9 d8 H+ }    a[i*n+v]=0.0;
" \6 `( \$ M1 R" r! |/ [9 x2 ?    for(k=0;k&lt;n;k++)
     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
<>
}//最开始的for的结束的大括号
return;
  }</P>
) H- }6 w& ^  b5 [  B) b5 |* G
: P. t& _+ f& Y8 x% i, L3 N<>自己写的运行总是错误</P>

水木年华zzu

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