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

aj6249

<># include "math.h"
7 G" G" j4 E' U0 x- W  U7 f& y#include"stdio.h"</P>
, w( R. G4 d( A<>  void strp(a,p,c, n,u)
( T/ V) G- d6 U( P& z6 w      int n;
   double a[],p[],c[],u[];
{   
   int i,j,k,v;
   double sum,asum;
for(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
! o: h' g2 {+ D. _* A4 U       for(i=0;i&lt;n;++i)//初始化u[]全为零
' |: ~/ V9 c; m4 ?  b     u=0.0;</P>
* m; z- y2 W' U  A2 g: s<>
      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
- O/ ]. l" d0 A' h0 z      {
       k=i*n+j;
    sum+=a[k]*a[k];
6 l* A9 H, {& l" c& o+ u" V: }" s   }
  v6 E0 a, {  K2 k' ~6 `+ L      asum=sqrt(sum);</P>
( P3 v) G- G: a/ j* q0 W/ |, P/ q8 ]<>      for(i=j+1;i&lt;n;++i)
1 f: s; O; o  ?& [   {
   
2 y( a* @  c9 l# ^# t    if(i==j+1)
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
; Z) Z7 A$ p) F: |# w/ j' d    else
/ C, i" L& y, i2 J8 R     u=a[i*n+j];
   }
      </P>
<>
   sum=0.0;  //实现P
   for(i=0;i&lt;n;++i)
" _* c- _' o' ?0 {7 d' R0 E6 `2 m! f* Y    sum+=(u*u);</P>
) C- y# m7 Z' v" T4 B<>   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]);}
         
   printf("\n");}</P>

* t: A, L8 b0 U7 R4 v$ s& U/ o' H<>
% V( }! L  c: X% u  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
+ \; S$ O( j7 k        for(v=0;v&lt;n;v++)
  {  c[i*n+v]=0.0;
   for(k=0;k&lt;n;k++)
    c[i*n+v]+=p[i*n+k]*a[k*n+v];
  }</P>
<>
     for(i=0;i&lt;n;i++)
   for(v=0;v&lt;n;v++)
   {
    a[i*n+v]=0.0;
    for(k=0;k&lt;n;k++)
     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
. _& ~3 H" W6 M- q+ j<>
}//最开始的for的结束的大括号
* y' ~$ b" V/ P9 M  H% Z return;
  }</P>
! A* g$ \( k1 O, i$ }) C
<>自己写的运行总是错误</P>

水木年华zzu

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