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

aj6249

<># include "math.h"
+ W6 t* R, V  ?6 H: o+ B#include"stdio.h"</P>
. r3 g+ q1 `0 g6 g( u! k<>  void strp(a,p,c, n,u)
      int n;
5 B" [0 L4 \8 \+ s$ N, n8 d   double a[],p[],c[],u[];
7 o0 x( L" p  a! J4 c& `3 f{   
' V, ]9 a4 {/ C$ s   int i,j,k,v;
# y# B" j' q, h% u- g   double sum,asum;
% W$ E5 J) C2 Efor(j=0;j&lt;n-2;++j)
5 P/ x* G  @$ d- Y) A{// 最开始的for 循环
' K; P) Y1 y# {/ B/ ^       for(i=0;i&lt;n;++i)//初始化u[]全为零
     u=0.0;</P>
<>
' P, q, ~- S- p: N0 ]3 T, ]/ T      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
      {
) S9 q0 M6 e/ z# ?8 B4 q' K( w& P       k=i*n+j;
    sum+=a[k]*a[k];
9 V' K) u% ?, \; ~# v/ Y   }
      asum=sqrt(sum);</P>
  f  c# g7 }, k) H+ q' s& a9 z<>      for(i=j+1;i&lt;n;++i)
9 b- f7 F; Q3 l   {
; Q  R* j2 @! T. T. i/ J   
    if(i==j+1)
- S* h+ x( F' a3 e3 ~$ u     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
2 {% w5 }- o* ?/ Z1 s    else
     u=a[i*n+j];
   }
      </P>
& e$ p  S3 b# ^' h<>
! x1 ^' G8 `! r; A   sum=0.0;  //实现P
   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)
# a6 T4 ^: V: @5 }  I* R7 e    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
    printf("%13.7e  ",p[i*n+k]);}
5 z- ]$ Q! f5 g# r0 G8 I         
7 @7 b& O4 f8 a8 T& T   printf("\n");}</P>
' c/ p$ Y+ j" @6 K
- B" W) E1 \8 c: P; }. p<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
        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>
0 m- H2 ~$ f7 f. |/ f<>
     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% Y- J# R3 _- s! p- n/ j     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
<>
}//最开始的for的结束的大括号
return;
  }</P>
/ @, H, V5 {8 B$ t; m( y( F
  {6 z' a* _* k# `, m9 Z- J" B<>自己写的运行总是错误</P>

水木年华zzu

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