数学建模社区-数学中国

标题: (求助C语言实现Householder变换一般实矩阵为上Hessenberg矩阵的算法) [打印本页]
作者: aj6249    时间: 2005-4-24 22:33
标题: (求助C语言实现Householder变换一般实矩阵为上Hessenberg矩阵的算法)
<># include "math.h"
) _" j& n5 C/ D) V3 J* J5 F8 U#include"stdio.h"</P>
# ~3 t, H3 w* ?<>  void strp(a,p,c, n,u)" s6 U. z) P2 \8 A* I( S
      int n;
6 w1 P; t: ^' _  n" B   double a[],p[],c[],u[];
- C9 O; ]7 v( Q- |& k0 ^+ W{   ! C$ w" x& |& x8 F+ _
   int i,j,k,v;
3 Q' Q) k: W1 c   double sum,asum;
6 ~1 x: B6 J+ N3 O- Dfor(j=0;j&lt;n-2;++j). S/ [0 |5 ]) v8 d: J
{// 最开始的for 循环2 X7 J0 X" J' p; `3 }
       for(i=0;i&lt;n;++i)//初始化u[]全为零
) G: Z. P4 `. u& X1 i( U) \     u=0.0;</P>9 {7 P+ \6 T4 i! R" j- H6 d
<>8 T7 k/ B5 U4 O# N1 K
      sum=0.0;
+ U% t  ]& f; i1 I6 A/ x) F   for(i=j+1;i&lt;n;++i)//实现a' g, D7 _2 X) R) d# D9 I$ F
      {9 C9 U8 ?9 e6 c$ H
       k=i*n+j;# O' f5 A; i5 ~, z, B- J) g
    sum+=a[k]*a[k];- w  ]/ M& V7 o# P  E
   }
) Z! s/ C* T) O: S) ?3 r      asum=sqrt(sum);</P>
! z: M8 D2 ?- w+ d9 L<>      for(i=j+1;i&lt;n;++i)+ j2 N  r5 D6 o) ?
   {; c+ |6 E2 {" ?: @* K
   
" a1 s4 h/ M8 d+ a% {1 s( s    if(i==j+1)) M% C/ b& l& I9 p
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
' n& a7 \' N9 |    else0 z7 J! ]" O7 @. }/ J) n1 y
     u=a[i*n+j];8 b' m% R: M7 s) T
   }
3 j) K& }3 {+ r; t/ E7 t      </P>
5 A. I- q- a) u% r2 h# g<>
- i  H4 _: p8 `. V   sum=0.0;  //实现P. E- _+ T5 g; V$ c3 p9 U. F3 R
   for(i=0;i&lt;n;++i)
9 f1 l: V# {& h2 ~+ r$ f) _5 @    sum+=(u*u);</P>" [+ j7 p: }/ }2 a
<>   for(i=0;i&lt;n;++i)% j2 c, `1 X( U7 z3 V
   {for(k=0;k&lt;n;++k)$ Z* }9 t9 O. X& T% I7 q, `/ b& x; `- @
    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;+ [' o  p( ~7 F7 ~  b
    printf("%13.7e  ",p[i*n+k]);}9 H# F% ]4 F& P; E, S. g1 p
         
' m1 C8 d/ G- v3 M. w4 ?: {2 b; I   printf("\n");}</P>8 ?: J# X9 o' D4 v/ P9 o

# x  b* e1 g! Y$ j4 o- {<>- |1 D, |' c0 V4 o$ ~
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
. P, C6 Y8 c  \        for(v=0;v&lt;n;v++)
- d5 K& x" X9 [) K  {  c[i*n+v]=0.0;8 T. Y! z9 F( Z. k/ Q1 {# [+ W# C# o" P1 |
   for(k=0;k&lt;n;k++)
, ?* ]. Z( a$ k3 [0 c    c[i*n+v]+=p[i*n+k]*a[k*n+v];
' k. X% j# S3 l9 t9 f6 P  }</P>0 b* \! o2 s1 G( E
<>
7 z1 [& J8 {& f2 _     for(i=0;i&lt;n;i++)4 }2 A  S: p, y- x* j& u
   for(v=0;v&lt;n;v++)* c$ }# Q) G6 d+ m
   {
- s& N  Y& {1 Y6 j7 O    a[i*n+v]=0.0;8 \( I8 r+ g6 m) T: _/ g4 B; ]
    for(k=0;k&lt;n;k++)
; S' P0 L, G: r+ F  f/ b     a[i*n+v]+=c[i*n+k]*p[k*n+v];3 N) X" }1 U2 P; F
   }</P>
& p" }- k0 ?6 ]4 ]<>, X8 L' ^/ z& K4 s5 w1 K, W
}//最开始的for的结束的大括号
& A6 l% V! y7 x6 Q- w2 J  {$ n return;
( c" `' C3 R( `/ l4 N  A' s1 i) k7 O  }</P>
& F& d, g% I; p9 [& M" I
: X  B' Q& {2 M* F' @; G<>自己写的运行总是错误</P>
作者: 水木年华zzu    时间: 2009-1-20 23:24
我已经把源代码发在计算数学板块了，你自己找下



欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) Powered by Discuz! X2.5