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

aj6249

<># include "math.h"
2 o  r# ~+ O. \# B#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
: I. A1 T0 R- V      int n;
5 ~3 R) K! U3 B+ d   double a[],p[],c[],u[];
{   
   int i,j,k,v;
: l$ R' [, E" ^  p( c. w   double sum,asum;
" U% ]. B& w6 t6 W0 \) B3 p7 Wfor(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
( X2 U, P, W7 T       for(i=0;i&lt;n;++i)//初始化u[]全为零
     u=0.0;</P>
& y( B8 z: j! X+ F% w<>
      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
      {
       k=i*n+j;
    sum+=a[k]*a[k];
   }
      asum=sqrt(sum);</P>
<>      for(i=j+1;i&lt;n;++i)
3 S0 I- X9 Y  d  `' r" ]   {
% Y0 d0 Z* N; o4 P$ ?   
    if(i==j+1)
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
    else
     u=a[i*n+j];
   }
# N+ U+ @* v( I% \; X  |7 Q      </P>
; l* d, m$ ~7 h& W4 E  s+ M2 \<>
- c" s" [( o* @" K   sum=0.0;  //实现P
* t* d9 M8 B6 j; S) p   for(i=0;i&lt;n;++i)
5 g4 {4 P& Q1 i1 v/ R0 r    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
3 O9 ^! t0 S5 i4 E% z/ Z% @! W, c   {for(k=0;k&lt;n;++k)
* g/ b0 O2 l+ d& d: J4 ~    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
    printf("%13.7e  ",p[i*n+k]);}
/ e6 l. o* w+ f, S/ a         
   printf("\n");}</P>

<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
        for(v=0;v&lt;n;v++)
, M+ t# x4 B, n6 Y/ I2 N6 G  {  c[i*n+v]=0.0;
   for(k=0;k&lt;n;k++)
! V3 Z* w" H3 C1 D' z" H    c[i*n+v]+=p[i*n+k]*a[k*n+v];
  }</P>
# |# I4 E( y! ~<>
     for(i=0;i&lt;n;i++)
& D* F) o7 H9 x& O   for(v=0;v&lt;n;v++)
4 N7 A+ K  L3 L3 A5 [   {
    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>
<>
}//最开始的for的结束的大括号
$ {% |6 q' S1 ~8 O3 E  v# i, n3 H) F return;
2 a) H- v) V0 V- K1 Q9 V! [  }</P>

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

水木年华zzu

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