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

aj6249

<># include "math.h"
! i  Q) s# L* ?( J#include"stdio.h"</P>
, q9 B! {1 ?& ]) w) y5 x! m. h/ D<>  void strp(a,p,c, n,u)
      int n;
   double a[],p[],c[],u[];
{   
   int i,j,k,v;
& D1 E: z4 k' _4 ~   double sum,asum;
for(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
- U0 n$ ^& u% Y$ z- l       for(i=0;i&lt;n;++i)//初始化u[]全为零
. a# B3 Y, p( M- t     u=0.0;</P>
<>
2 H* `% _; j8 r/ [  H" b' D      sum=0.0;
% V) `$ W& e" W# V; Y   for(i=j+1;i&lt;n;++i)//实现a
      {
  W# h8 a7 m: X       k=i*n+j;
    sum+=a[k]*a[k];
6 O0 I8 R! v9 `0 `: }: f# k, h   }
  l. h4 Z* W" ^) e; _# K      asum=sqrt(sum);</P>
( E( H4 ?9 ~# b2 `9 U  {3 {8 Z4 `4 E<>      for(i=j+1;i&lt;n;++i)
   {
   
. `# @& o; ?  H( m6 {    if(i==j+1)
7 O: _/ l6 E" y& t# Q1 b     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
    else
. Q; X! |' l5 o7 ^" }1 i- T& t     u=a[i*n+j];
! u+ N9 ?; k# Q# K" b+ h   }
, l0 v0 z! V& P      </P>
<>
/ [7 {: z' _1 l. u6 l' _+ f# a   sum=0.0;  //实现P
' N# \2 V- }; T* b   for(i=0;i&lt;n;++i)
3 z3 k* w; v) h  p4 ?5 o/ k    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
   {for(k=0;k&lt;n;++k)
    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
    printf("%13.7e  ",p[i*n+k]);}
0 w+ L0 n2 ]- [# i# F         
( [( L( T  [+ _( T/ S* H4 E   printf("\n");}</P>

" a& X5 _: O& t5 ^! V0 ]( |( ^<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
        for(v=0;v&lt;n;v++)
  {  c[i*n+v]=0.0;
/ N! T. R# w: D- K+ `   for(k=0;k&lt;n;k++)
    c[i*n+v]+=p[i*n+k]*a[k*n+v];
9 l- h' g& N. v' z; d; }  }</P>
, f3 b- l3 `2 q<>
     for(i=0;i&lt;n;i++)
$ r4 C2 O" M& T   for(v=0;v&lt;n;v++)
   {
    a[i*n+v]=0.0;
0 y! p5 d2 o  i1 S5 K5 X! l    for(k=0;k&lt;n;k++)
- D! }) ~' C. a     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
" A- p1 o+ X2 z+ L/ _<>
  S- }  `! ~, y0 @% i}//最开始的for的结束的大括号
return;
  }</P>

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

水木年华zzu

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