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

aj6249

<># include "math.h"
#include"stdio.h"</P>
1 L0 M  v" L; `" C3 S0 w8 E<>  void strp(a,p,c, n,u)
' y" Y! v+ }- {( k0 W9 i2 i      int n;
   double a[],p[],c[],u[];
{   
3 _9 D- ]3 Z- B# k: g: O   int i,j,k,v;
   double sum,asum;
for(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
7 l: n! v! F6 y: n" Q0 b6 y% F- w1 S       for(i=0;i&lt;n;++i)//初始化u[]全为零
5 u! ^0 Y& X( r     u=0.0;</P>
' f/ e7 N8 A. E5 O& X$ p( H<>
! o" g9 E7 o* j, c( N" n  b      sum=0.0;
0 `1 Z1 [  m8 M9 y4 J   for(i=j+1;i&lt;n;++i)//实现a
8 U- n* I5 e5 ^% b2 P. i      {
       k=i*n+j;
    sum+=a[k]*a[k];
   }
      asum=sqrt(sum);</P>
<>      for(i=j+1;i&lt;n;++i)
; {$ {5 t8 j7 J+ h) _- {   {
' s5 T( l+ ?& p7 I9 I% g9 W- y   
1 v0 I. }# {! Z, N7 n) b6 E    if(i==j+1)
! B+ S  u0 w; ]" C6 n     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
. y4 N2 \6 X& i# x    else
4 \8 d7 a' X! M9 E$ e) v     u=a[i*n+j];
   }
      </P>
<>
% N4 R2 q. W- v" T! a8 d4 O" \   sum=0.0;  //实现P
   for(i=0;i&lt;n;++i)
6 ^; w7 Y/ m( s3 @0 ^3 {    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
1 k1 a- T% ?+ e) t6 u# |) k   {for(k=0;k&lt;n;++k)
    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
& k: Q" V% n6 K3 z    printf("%13.7e  ",p[i*n+k]);}
! \3 T6 s6 A' f# P# R$ D0 |         
7 g( Z; d" e) K! ?+ ^   printf("\n");}</P>
( ~0 K' _$ y* J9 H) V6 b4 O, L
2 z" n" F, c; A2 H* r<>
- N, n  u$ N0 I  ]7 F+ s# n4 |" D  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
/ C5 ]% _2 j7 Z6 o- p        for(v=0;v&lt;n;v++)
  {  c[i*n+v]=0.0;
   for(k=0;k&lt;n;k++)
7 T: q1 ~: P4 }" l' @    c[i*n+v]+=p[i*n+k]*a[k*n+v];
9 h) E/ V6 b: ~  }</P>
$ Z' t* D% S  I+ A0 Z* g+ N3 @; G4 D<>
     for(i=0;i&lt;n;i++)
   for(v=0;v&lt;n;v++)
- O6 F  H* {) i' v. J" h+ ?   {
6 g1 l5 n" n5 x; ?( n0 ?    a[i*n+v]=0.0;
8 Q  K$ K( R* i& @    for(k=0;k&lt;n;k++)
     a[i*n+v]+=c[i*n+k]*p[k*n+v];
7 T3 C  m6 ~% M: }" c) M   }</P>
<>
0 A) |1 _. [$ `4 d9 t8 ~}//最开始的for的结束的大括号
return;
1 ~8 s; ]4 E. f' g& M  }</P>

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

水木年华zzu

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