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

aj6249

<># include "math.h"
; M1 [- {* Q* y. P6 J1 y#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
7 i; ?9 C' [" A      int n;
   double a[],p[],c[],u[];
{   
" z; S0 k$ w7 r% K$ r4 R0 w/ B   int i,j,k,v;
   double sum,asum;
5 t7 u+ G1 o* \! z7 X3 h& }for(j=0;j&lt;n-2;++j)
6 H$ _2 v7 L3 d( {+ }+ y{// 最开始的for 循环
       for(i=0;i&lt;n;++i)//初始化u[]全为零
, e5 i% r) R& O3 v     u=0.0;</P>
) _' b+ C- F! ~1 t5 W" d5 b<>
      sum=0.0;
  _  {+ ?" k0 I; r& p- P7 M& M   for(i=j+1;i&lt;n;++i)//实现a
" O- v6 b# g: I# H! E/ U      {
       k=i*n+j;
  `/ E# C' C* a: Z8 y. [    sum+=a[k]*a[k];
   }
      asum=sqrt(sum);</P>
, {" i; s$ G% X* j" S- I  e<>      for(i=j+1;i&lt;n;++i)
   {
' _6 L7 d* P' \: k; }+ {. a8 w   
    if(i==j+1)
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
    else
5 r1 ^+ ?6 {2 m* M     u=a[i*n+j];
   }
      </P>
<>
   sum=0.0;  //实现P
, d9 X" S) F6 D5 N3 G8 k   for(i=0;i&lt;n;++i)
    sum+=(u*u);</P>
: ]& X0 i) k& Q* c! `/ t% O  C/ F2 z<>   for(i=0;i&lt;n;++i)
) ~# Z1 d/ h" N2 R: c* n: K   {for(k=0;k&lt;n;++k)
) y' J! x, A+ z: e, z, r$ Y5 _    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
! d8 e/ L) \. ]3 x7 \0 O    printf("%13.7e  ",p[i*n+k]);}
         
   printf("\n");}</P>
1 M: f" M+ P0 J# t
<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
8 ~9 i- F/ Q. X6 p        for(v=0;v&lt;n;v++)
  {  c[i*n+v]=0.0;
3 E+ J4 y  w: X; Y( y0 n' }" H   for(k=0;k&lt;n;k++)
+ P0 |; H" F9 `    c[i*n+v]+=p[i*n+k]*a[k*n+v];
  }</P>
<>
     for(i=0;i&lt;n;i++)
- F5 s$ D  ]) E6 t5 s0 `   for(v=0;v&lt;n;v++)
   {
    a[i*n+v]=0.0;
3 d% f, v5 ]7 u2 t    for(k=0;k&lt;n;k++)
     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
<>
0 `3 m, o- H! f0 x  P1 v6 m0 @}//最开始的for的结束的大括号
return;
$ N, q. ]: C& d2 O  U" a5 \1 c  }</P>
5 Y+ b( p4 b  L+ J5 G+ f5 y- X
, `9 H8 W) \, t+ {" V4 m, E<>自己写的运行总是错误</P>

水木年华zzu

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