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

aj6249

<># include "math.h"
1 `* d8 S; G1 X3 ~" p#include"stdio.h"</P>
<>  void strp(a,p,c, n,u)
& ~/ U; F* W8 {9 X& N$ ^, ^      int n;
* \; {9 i% l$ U( z: |! {   double a[],p[],c[],u[];
{   
   int i,j,k,v;
% e0 b/ V  R$ b; k; j7 K   double sum,asum;
for(j=0;j&lt;n-2;++j)
{// 最开始的for 循环
       for(i=0;i&lt;n;++i)//初始化u[]全为零
0 @' g9 J" r) j0 w     u=0.0;</P>
% J- k$ F8 E; x: c<>
      sum=0.0;
- n  P& O& u- ?( k   for(i=j+1;i&lt;n;++i)//实现a
! `3 S" m% K: p      {
8 e8 g$ j* {* {+ H! Z       k=i*n+j;
8 c' j9 h7 n, j: L0 L$ }    sum+=a[k]*a[k];
   }
) c, R$ w: h$ ?6 @: A; S6 J6 f, {      asum=sqrt(sum);</P>
<>      for(i=j+1;i&lt;n;++i)
   {
' u0 K: t4 n1 H' l   
, N/ F: Y3 U# z0 I5 j. y5 ?    if(i==j+1)
: L, j: R7 L$ y7 R     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
    else
     u=a[i*n+j];
   }
! H( O/ c1 t7 {1 }      </P>
<>
   sum=0.0;  //实现P
   for(i=0;i&lt;n;++i)
    sum+=(u*u);</P>
<>   for(i=0;i&lt;n;++i)
* H; j7 |+ _) U3 `4 b: k   {for(k=0;k&lt;n;++k)
( X3 \  \$ p3 P" n' N% S" N8 X$ }    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
0 S4 [8 t$ y: A    printf("%13.7e  ",p[i*n+k]);}
         
   printf("\n");}</P>

<>
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
2 h3 T/ X, O3 N, O% N5 W        for(v=0;v&lt;n;v++)
: G/ t# K' s6 i8 l  {  c[i*n+v]=0.0;
   for(k=0;k&lt;n;k++)
# E3 Y$ f8 z. T1 h* `. |) ~2 D# I# J7 N    c[i*n+v]+=p[i*n+k]*a[k*n+v];
' _: `5 C/ v- g( E: p6 r# x  }</P>
$ \* X6 K1 a3 I% S% s2 c<>
! t% c' d, B# b( M/ Y0 M0 t" p     for(i=0;i&lt;n;i++)
   for(v=0;v&lt;n;v++)
   {
# h2 k1 U; {3 f$ C, ]8 o* l    a[i*n+v]=0.0;
0 |  N" F( c/ C- E    for(k=0;k&lt;n;k++)
     a[i*n+v]+=c[i*n+k]*p[k*n+v];
! L* ~* K3 _: p8 I- T* w   }</P>
5 u. V  Q' `) d1 s<>
}//最开始的for的结束的大括号
return;
  }</P>

6 F1 \% H. a- r. @! W4 v# [% T/ A( v<>自己写的运行总是错误</P>

水木年华zzu

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