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

aj6249

<># include "math.h"
* L# _0 E. G' W$ a$ X#include"stdio.h"</P>
* g5 ^- ?! d5 g1 |<>  void strp(a,p,c, n,u)
      int n;
   double a[],p[],c[],u[];
{   
8 K5 k9 g% R9 T  o" i4 ?   int i,j,k,v;
   double sum,asum;
for(j=0;j&lt;n-2;++j)
4 S( M, Y( q& R0 _+ n8 G& h{// 最开始的for 循环
4 o" Z: V: G# D. F' q       for(i=0;i&lt;n;++i)//初始化u[]全为零
     u=0.0;</P>
7 W6 ]' l& n7 |$ q) f<>
      sum=0.0;
6 a- D$ w0 A3 Y   for(i=j+1;i&lt;n;++i)//实现a
      {
) o- j+ C* i7 X' Y6 `. S       k=i*n+j;
, `8 E9 M7 m% k    sum+=a[k]*a[k];
1 A3 b' d9 ]# B% G3 c' f   }
      asum=sqrt(sum);</P>
<>      for(i=j+1;i&lt;n;++i)
   {
   
; Y9 N0 X# j5 B4 G# S0 {9 z. z4 e4 E    if(i==j+1)
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
    else
5 r: u# E" ^' F7 w+ X2 g     u=a[i*n+j];
   }
- }9 u6 M: C$ _4 N' I      </P>
' d: g( s* k% s, i6 H& ~' _<>
   sum=0.0;  //实现P
" Z3 b2 m! @5 ?   for(i=0;i&lt;n;++i)
    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
   {for(k=0;k&lt;n;++k)
8 r4 \2 s0 d0 V- ]  Q    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
    printf("%13.7e  ",p[i*n+k]);}
         
6 F! [0 ]5 _5 g0 @8 }   printf("\n");}</P>
% ]0 J8 u7 _/ Q0 p* p* }; Q3 Z
9 \3 Y/ W* X: I<>
- y/ {; T, i4 U$ k  m  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
        for(v=0;v&lt;n;v++)
  {  c[i*n+v]=0.0;
: J: t; A; S$ B3 v9 `; v0 {   for(k=0;k&lt;n;k++)
8 ~8 u% @: H2 [* A8 b' b    c[i*n+v]+=p[i*n+k]*a[k*n+v];
  }</P>
% w/ Q& _2 Q; M9 v- X/ Q) i<>
" F' K- z) U+ I     for(i=0;i&lt;n;i++)
" M6 W3 y) g! B* ~: ^+ @3 |1 f   for(v=0;v&lt;n;v++)
   {
    a[i*n+v]=0.0;
    for(k=0;k&lt;n;k++)
: G' n. p( A3 [: n- H3 g" x* s! s     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
<>
, o; D) Z2 y8 g}//最开始的for的结束的大括号
* d" C- Y: u+ M) A, a$ I return;
# e) A% N' l4 f6 S3 u; \) W6 q) W, X2 h  }</P>

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

水木年华zzu

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