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

aj6249

<># include "math.h"
#include"stdio.h"</P>
# B  ^; b7 `/ K5 n, l& v- u, q: _1 G<>  void strp(a,p,c, n,u)
      int n;
   double a[],p[],c[],u[];
! U7 t! h1 N* k# N: e7 `+ a{   
   int i,j,k,v;
7 T* N5 r* Q$ M8 b! F* r% d   double sum,asum;
0 q6 l" S, B/ X1 Rfor(j=0;j&lt;n-2;++j)
* N& N1 S- I4 v% F  Q! m{// 最开始的for 循环
       for(i=0;i&lt;n;++i)//初始化u[]全为零
" v7 A3 L0 K& c* S8 O# d' q; _     u=0.0;</P>
! S3 f+ t% e5 {<>
' b9 ?7 R4 R( }" R' q  [      sum=0.0;
* E8 R" O  M* ?4 @: l/ u& d$ }  O   for(i=j+1;i&lt;n;++i)//实现a
      {
       k=i*n+j;
    sum+=a[k]*a[k];
* V: `: P! i  o5 [* b   }
      asum=sqrt(sum);</P>
<>      for(i=j+1;i&lt;n;++i)
   {
   
    if(i==j+1)
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
& ~# M7 Y- L* R- n  b1 ~* I    else
     u=a[i*n+j];
" P0 i5 V7 z0 B/ O$ N   }
      </P>
; `3 R8 ?% q+ c/ v% c<>
* W$ d* O0 |9 W3 |  {/ {4 I   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)
- p. F3 L" D& M9 @    {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>

5 u$ Z1 j  ~3 b3 B<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
5 r0 o, M1 j. j        for(v=0;v&lt;n;v++)
( d" `2 t3 r3 ~) u, n8 k+ F, k  {  c[i*n+v]=0.0;
: T% T/ f( J' w& F. q8 g& R   for(k=0;k&lt;n;k++)
6 H* C* a" W  ^4 c' {    c[i*n+v]+=p[i*n+k]*a[k*n+v];
  }</P>
<>
     for(i=0;i&lt;n;i++)
3 K: i/ b& Q/ z; @# |   for(v=0;v&lt;n;v++)
0 O: k2 e' K1 z5 n, f4 U   {
    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];
2 H3 e7 z% h. t, d( }& C1 S   }</P>
<>
}//最开始的for的结束的大括号
0 ]% R0 r0 H. ~6 E; Q0 m return;
  }</P>

- v) I8 u8 q& ^& ?9 x<>自己写的运行总是错误</P>

水木年华zzu

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