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

aj6249

<># include "math.h"
$ u  f# Y9 l' t, \/ M#include"stdio.h"</P>
7 M' _0 S! |  o. f5 J% J! B<>  void strp(a,p,c, n,u)
; ?9 I. N9 K) D7 _1 U      int n;
( ~  S0 ~; u, b3 X# `   double a[],p[],c[],u[];
{   
   int i,j,k,v;
   double sum,asum;
" [0 {" e+ ~! C; }" _. A6 Afor(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
0 i6 R/ q0 v/ E8 s5 u1 U       for(i=0;i&lt;n;++i)//初始化u[]全为零
# Y" R" I- g$ z3 a     u=0.0;</P>
! S& X0 q4 c& m& q<>
% J9 j* ]) ]( Y" j) E      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
% P% n7 e4 r& h' N1 G      {
3 R1 T7 q9 w3 U% b       k=i*n+j;
    sum+=a[k]*a[k];
7 L: S0 N% B3 \4 C* Z2 G  H2 `- U   }
/ {$ |* b) h/ @; o: p8 X      asum=sqrt(sum);</P>
7 K( K+ J% A) @9 I; f5 C- J* Z/ `  ~<>      for(i=j+1;i&lt;n;++i)
1 A1 y2 G5 y1 ?$ v$ C/ V   {
+ ~& |) }6 C* m   
    if(i==j+1)
- `0 X$ m8 E- \9 r     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
    else
     u=a[i*n+j];
* V( v+ y4 G# S5 y+ m   }
      </P>
* D1 r, J  q4 K8 h0 J' [<>
   sum=0.0;  //实现P
   for(i=0;i&lt;n;++i)
  d) s" }" R* h    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;
9 C- F: @7 m0 Y    printf("%13.7e  ",p[i*n+k]);}
         
   printf("\n");}</P>

<>
$ F5 U+ q! Q5 H9 T( W  B  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
( U9 R# v7 R- W6 N# O8 d        for(v=0;v&lt;n;v++)
1 z- F% U+ t5 W" L) Q  {  c[i*n+v]=0.0;
+ q6 m% a/ G& L   for(k=0;k&lt;n;k++)
    c[i*n+v]+=p[i*n+k]*a[k*n+v];
* x9 P) x$ V& R  }</P>
' e2 w: f! ^0 c* K' N; ^<>
* v! ~( t0 ~! l( A     for(i=0;i&lt;n;i++)
   for(v=0;v&lt;n;v++)
   {
    a[i*n+v]=0.0;
$ N% V+ {: V) O* ^2 X4 Z1 }    for(k=0;k&lt;n;k++)
     a[i*n+v]+=c[i*n+k]*p[k*n+v];
# b8 N$ Q' Q/ m6 x   }</P>
4 h  V7 |' r/ L! T<>
}//最开始的for的结束的大括号
9 K3 O0 ?" G7 X1 E, }) f; M return;
  }</P>
7 |2 U$ k; q; @
, f8 @% a5 n: m<>自己写的运行总是错误</P>

水木年华zzu

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