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

aj6249

<># include "math.h"
) V# g" b8 D9 V#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
      int n;
   double a[],p[],c[],u[];
{   
   int i,j,k,v;
   double sum,asum;
for(j=0;j&lt;n-2;++j)
; x  P# ]: i/ w{// 最开始的for 循环
2 I9 q2 U4 T) j$ v       for(i=0;i&lt;n;++i)//初始化u[]全为零
     u=0.0;</P>
<>
2 n7 ]( S4 h% v" U8 Q  A      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
) B2 u6 C( v- y  t1 C2 g/ y      {
       k=i*n+j;
+ t6 p5 G" T7 I! D    sum+=a[k]*a[k];
' C3 x1 w, B# n" j+ C7 w   }
      asum=sqrt(sum);</P>
5 u  m1 w- f- I9 |<>      for(i=j+1;i&lt;n;++i)
' v- C* l" s9 E  s   {
   
    if(i==j+1)
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
    else
# r' _7 Z+ F9 @; h     u=a[i*n+j];
   }
% A: k7 ], q( F8 w' d2 F      </P>
1 n8 {, r- T1 s) J6 R, U4 Y<>
   sum=0.0;  //实现P
% X9 B* b4 G( x' ?   for(i=0;i&lt;n;++i)
% v7 R% z8 }7 l' q- D$ q1 u& r, i8 n    sum+=(u*u);</P>
+ n, Y7 _: q% c<>   for(i=0;i&lt;n;++i)
) Z9 x# H% d' m# x   {for(k=0;k&lt;n;++k)
* o0 l8 l4 ]8 s" W# q    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
( W& R, Y+ y+ V& o- h7 Q    printf("%13.7e  ",p[i*n+k]);}
         
2 c4 T$ l) q$ P0 K2 I$ t   printf("\n");}</P>

& t+ T0 ?3 |8 J9 P5 B0 s4 K<>
2 F8 ~- G+ K6 H* U& E2 s: O  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
        for(v=0;v&lt;n;v++)
  {  c[i*n+v]=0.0;
& S! M0 u* G0 D3 @/ E   for(k=0;k&lt;n;k++)
3 {1 H% g# B. X1 p, w    c[i*n+v]+=p[i*n+k]*a[k*n+v];
  }</P>
6 l4 V, D% Y1 `/ r: H; K* C<>
     for(i=0;i&lt;n;i++)
7 V0 X' e- h) ^: x! ^   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];
" D3 L8 r& D. t1 @- m3 D3 Y$ N" Y8 I   }</P>
<>
}//最开始的for的结束的大括号
/ d0 [; ?1 I: t return;
  }</P>

: C0 a7 o: \! K2 g9 u$ e( g9 \( h<>自己写的运行总是错误</P>

水木年华zzu

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