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

aj6249

<># include "math.h"
#include"stdio.h"</P>
# S4 h4 N  @5 C<>  void strp(a,p,c, n,u)
      int n;
. V  ]9 X% j# ?; y2 c0 c. [3 T   double a[],p[],c[],u[];
{   
* o* @, E- Z# l. U; x0 U" E% I   int i,j,k,v;
: _+ A  V/ A3 t! ]0 n. w  p   double sum,asum;
for(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
       for(i=0;i&lt;n;++i)//初始化u[]全为零
( w+ V( G8 e7 \/ t     u=0.0;</P>
<>
; s  Z$ D9 G/ y' y8 `) q# A/ }      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
      {
& l* P7 c) l& q( i. u( ^       k=i*n+j;
    sum+=a[k]*a[k];
   }
      asum=sqrt(sum);</P>
<>      for(i=j+1;i&lt;n;++i)
7 t' c6 z; x$ D1 I& f   {
   
) {. f/ Y# K/ z; J- P2 E    if(i==j+1)
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
* z2 B1 _4 H3 G% \$ h    else
# G6 K. C* d/ x! H+ b     u=a[i*n+j];
   }
* S/ X* f- y. |3 @4 b      </P>
<>
   sum=0.0;  //实现P
   for(i=0;i&lt;n;++i)
$ `9 b' A; k) \9 s: \+ s    sum+=(u*u);</P>
. w5 m" v, N1 A- o<>   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;
8 v  ^: C7 h+ a6 D" o4 Y9 k    printf("%13.7e  ",p[i*n+k]);}
! a7 l0 h( p6 ?+ N' H         
   printf("\n");}</P>
$ ]! k9 I9 @  V6 X2 d5 s6 C( C
<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
: U+ |7 `9 w* I4 H        for(v=0;v&lt;n;v++)
0 o" o, [3 K9 d9 v, B+ H  {  c[i*n+v]=0.0;
7 Y, o; A4 E: w( x. Q   for(k=0;k&lt;n;k++)
& O1 H  }5 N" e8 }8 T- F: o) H1 X    c[i*n+v]+=p[i*n+k]*a[k*n+v];
: R# |2 U: {. q2 j7 [* z  }</P>
<>
     for(i=0;i&lt;n;i++)
   for(v=0;v&lt;n;v++)
' _9 b7 S0 T* G+ y7 X2 S   {
+ ?8 I4 m4 b$ M5 D; Y$ t    a[i*n+v]=0.0;
    for(k=0;k&lt;n;k++)
, k, e( \" ]. P  H. u     a[i*n+v]+=c[i*n+k]*p[k*n+v];
/ {+ a  V" G7 X, M8 p& f) A   }</P>
; B& K! F0 D% `1 @0 J<>
}//最开始的for的结束的大括号
return;
  }</P>
7 S1 G2 a0 R" f' Q0 K' }
6 U) Q0 O4 }: G( e' I& p<>自己写的运行总是错误</P>

水木年华zzu

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