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

aj6249

<># include "math.h"
#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
( U3 a" p# M" Z* `/ K  i      int n;
   double a[],p[],c[],u[];
{   
   int i,j,k,v;
- B3 a7 X' q, x4 q/ q8 |9 t8 D2 _   double sum,asum;
5 _% \# }2 i8 q# B: W  Xfor(j=0;j&lt;n-2;++j)
( |/ g% j4 }  e2 L{// 最开始的for 循环
6 x3 f# d9 u5 |+ I! t# a       for(i=0;i&lt;n;++i)//初始化u[]全为零
     u=0.0;</P>
<>
      sum=0.0;
   for(i=j+1;i&lt;n;++i)//实现a
" x4 ^& |5 n8 _* X1 D( l# _& }      {
       k=i*n+j;
    sum+=a[k]*a[k];
( T5 t7 e* E* Z, d7 u3 K   }
- m+ b- A/ J* S; O      asum=sqrt(sum);</P>
1 r! E' q) R; S" e( z8 q7 Z<>      for(i=j+1;i&lt;n;++i)
   {
   
+ k! X3 x' U/ C$ F0 }+ y/ a* e# C    if(i==j+1)
: m' l; i3 i1 q, b     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
( e, s$ U! y0 `1 D& \6 L. i    else
     u=a[i*n+j];
   }
$ n% G- I) [6 Z! ~# v4 h+ N      </P>
<>
   sum=0.0;  //实现P
" g% [0 R0 \( y, q* Z$ T0 ?   for(i=0;i&lt;n;++i)
* e, T' g* X! K" P% V6 C( d4 S5 }# r  N: r    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
  @2 t* g' c, O  y" v7 R   {for(k=0;k&lt;n;++k)
1 }- u. W* B  I8 ^- F7 `    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
    printf("%13.7e  ",p[i*n+k]);}
: l9 P% \( b0 j* g  @         
   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++)
6 \3 N  J2 T& J" P& k0 {4 \8 I    c[i*n+v]+=p[i*n+k]*a[k*n+v];
8 T. V3 D0 ^4 T" a% s  }</P>
<>
( e- U9 A& Y1 T9 B9 v     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++)
+ y/ {! y  h3 x1 D  Z     a[i*n+v]+=c[i*n+k]*p[k*n+v];
   }</P>
2 _% K# l* e2 e<>
}//最开始的for的结束的大括号
return;
! y/ ~" }5 X3 n6 E, k- U  }</P>
2 X, n! _6 f1 J$ c
5 I7 `: i5 c: v! @<>自己写的运行总是错误</P>

水木年华zzu

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