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

aj6249

<># include "math.h"
2 ~# L6 B. u# d% Q! q#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
9 O( l0 t2 _* g+ Q2 d6 S      int n;
+ b7 u7 C3 N6 o   double a[],p[],c[],u[];
{   
& I7 H5 s' i7 T# F   int i,j,k,v;
3 I0 f- P6 m5 N0 {) Z   double sum,asum;
for(j=0;j&lt;n-2;++j)
( m) V% C5 r5 x% N{// 最开始的for 循环
  D) ?/ u8 U# P/ A/ Q9 ~! c5 ?7 q1 f       for(i=0;i&lt;n;++i)//初始化u[]全为零
     u=0.0;</P>
  A8 L: r2 h6 q5 R" p" p<>
- a, U7 y1 }  t; L% b3 c8 p! e# Y      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
      {
       k=i*n+j;
; d! D: R4 b( i2 R5 j, W    sum+=a[k]*a[k];
8 c+ w: V! i) Z" a0 [   }
      asum=sqrt(sum);</P>
+ _9 C# l7 X3 q1 w<>      for(i=j+1;i&lt;n;++i)
   {
6 ?) w. n0 I+ M/ E& E% I2 q   
0 F0 o" F4 d# x, `! o. ^    if(i==j+1)
, H" o  Y+ |  T1 B     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
+ }+ n+ i% }" b; W/ }2 B$ i    else
     u=a[i*n+j];
5 G4 V  t! [5 s   }
      </P>
<>
! r2 C9 Z: `( l; O, N   sum=0.0;  //实现P
5 |$ J$ J0 T9 R& m2 a: z   for(i=0;i&lt;n;++i)
    sum+=(u*u);</P>
6 w8 h' I8 ]/ |! K$ a<>   for(i=0;i&lt;n;++i)
6 s/ D, A% X& c: i* L& S$ Q+ Z' d  b   {for(k=0;k&lt;n;++k)
& u1 j" H4 e; N4 L    {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>

<>
  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++)
+ O5 v5 u- |3 E6 t  G    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++)
# m& s! y  r8 S, }6 \; @  h, g; _   {
    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>
/ `$ B- @, {+ H, C5 d3 T9 O; W<>
}//最开始的for的结束的大括号
return;
8 S3 r8 Q" `9 Q* B; g: q% X  }</P>

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

水木年华zzu

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