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

aj6249

<># include "math.h"
#include"stdio.h"</P>
) \4 m- o% s0 E' z% g" g<>  void strp(a,p,c, n,u)
! [& y* Y% i6 p# i  J7 a( f" Q      int n;
   double a[],p[],c[],u[];
! z8 G6 I. n7 G1 a$ c; s. d9 {1 Q4 s" N{   
   int i,j,k,v;
   double sum,asum;
for(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
) b' @2 t  p. J3 t  @       for(i=0;i&lt;n;++i)//初始化u[]全为零
! H- i% J! b! s  f+ |1 N     u=0.0;</P>
% w6 |, b7 Z& i+ @% J! `- @# M<>
# @) w6 Y4 [8 F: e      sum=0.0;
) L4 t$ q! I* K3 h2 h6 C1 ~9 v2 b   for(i=j+1;i&lt;n;++i)//实现a
      {
9 j# N5 {: j1 ^# s       k=i*n+j;
    sum+=a[k]*a[k];
   }
" w, t2 G  H, e' V1 p0 g      asum=sqrt(sum);</P>
6 z6 j' Y6 ~/ E! I<>      for(i=j+1;i&lt;n;++i)
   {
   
. I' S( D$ d9 ^* F    if(i==j+1)
2 i2 W" w- y( d% ?     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
- ?) @5 l: f$ Q  X1 H    else
     u=a[i*n+j];
4 o  L& q8 ^1 L9 B1 ^1 U   }
/ S* k2 I* f0 {  i/ W      </P>
* o, L2 z- z( z: G- R/ X, X# w<>
- R" L) M* B3 t0 P   sum=0.0;  //实现P
* M6 h4 N& N: M% H# n   for(i=0;i&lt;n;++i)
    sum+=(u*u);</P>
" E/ a9 @; a0 b% F; ~; w<>   for(i=0;i&lt;n;++i)
   {for(k=0;k&lt;n;++k)
. {/ r- z4 M' A# x  t8 O( q( i    {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>
8 p( F8 ]8 C% u
<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
3 Q5 J! W& K2 _/ H5 {" `        for(v=0;v&lt;n;v++)
+ c) l5 R$ g( A: f2 M& n  {  c[i*n+v]=0.0;
+ O+ L  s  L" c   for(k=0;k&lt;n;k++)
% a) E  P: u& v" Z    c[i*n+v]+=p[i*n+k]*a[k*n+v];
+ ?* K$ B- v9 [7 K% R9 a: h! f  }</P>
<>
$ ^3 d. I7 H% \( W( F. I     for(i=0;i&lt;n;i++)
   for(v=0;v&lt;n;v++)
   {
    a[i*n+v]=0.0;
    for(k=0;k&lt;n;k++)
, G  V9 x/ _) X. A3 R9 I" ?' v     a[i*n+v]+=c[i*n+k]*p[k*n+v];
. v0 L9 X/ x$ s   }</P>
* P: y! L% ]! ]6 {! @$ p<>
0 B/ }3 K4 W# Y( X# D}//最开始的for的结束的大括号
6 S+ c7 [: {8 ^ return;
2 X4 }) ^& p- n8 E6 Q; y1 G  }</P>
5 O0 e( A/ E5 m" B
- Q9 w8 Z/ G9 w) j/ y% c<>自己写的运行总是错误</P>

水木年华zzu

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