数学建模社区-数学中国

标题: (求助C语言实现Householder变换一般实矩阵为上Hessenberg矩阵的算法) [打印本页]
作者: aj6249    时间: 2005-4-24 22:33
标题: (求助C语言实现Householder变换一般实矩阵为上Hessenberg矩阵的算法)
<># include "math.h"2 T* |, w/ }) l# a
#include"stdio.h"</P>) K% p! r3 o! q. D' C
<>  void strp(a,p,c, n,u)
* C/ ^2 P/ B( }; K! Z      int n;+ O- V/ A7 f( ?( c% u# _! x
   double a[],p[],c[],u[];" G$ c$ M1 C' r2 d. |
{   : u/ O+ A8 J; h- O4 z! }
   int i,j,k,v;
, i) L: D( `0 i- E, m   double sum,asum;: G1 b3 L% u3 z  n0 C
for(j=0;j&lt;n-2;++j)& ^: o& Z" k1 Y: }/ _5 g9 f1 v
{// 最开始的for 循环
+ m$ ^; v! z4 K! t/ N1 H       for(i=0;i&lt;n;++i)//初始化u[]全为零) D& a: G3 t1 I7 P$ m
     u=0.0;</P>& ~' ]( X, {+ e7 t* i& @
<>
9 A  \- m& i9 b7 @7 m) d' V& S6 I      sum=0.0;
: P8 }( N$ t2 j4 b   for(i=j+1;i&lt;n;++i)//实现a
: |2 O$ J, A; x, E      {; i7 H: C* ~- I) I- J& R, k
       k=i*n+j;
# L5 w0 F. C! X; t    sum+=a[k]*a[k];6 m* |$ W! _- e& |& y
   }
+ W& |5 t" u! L5 \6 A' Z      asum=sqrt(sum);</P>
* R7 P; u( @- v( B<>      for(i=j+1;i&lt;n;++i)) K4 t9 K$ E, j" l7 S/ q9 j0 N- t, i
   {
0 e/ l$ }0 _/ S% g4 Z   
2 |* j' u$ o' a# _$ R    if(i==j+1)" o0 R! z0 y6 k6 _$ w3 r) e
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
" g5 ^/ Y+ [0 a1 M# F4 g6 S: E, E    else
" N8 J2 i4 b7 M/ ]! D     u=a[i*n+j];" A8 @3 D# _6 Z- d* n" a$ G: u& c
   }
6 V# s' n+ z  [5 X. |9 U" \      </P>
0 l$ m& ?* e; [6 p1 O<>
- I. m  k( C% B+ A4 y) f   sum=0.0;  //实现P, ]2 a$ g# E7 T5 D6 F* _! j
   for(i=0;i&lt;n;++i)( y2 Q5 o+ K- K* b" |( l
    sum+=(u*u);</P>9 K; @% _; r7 S+ n4 y1 L% f
<>   for(i=0;i&lt;n;++i)
! ^0 u7 p) Z& C% K) q3 j" R   {for(k=0;k&lt;n;++k)
3 b. m- ]% J, k/ R( u7 G3 D& d" [2 ?    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;: t: T7 B: R$ g6 D9 z% ?. N
    printf("%13.7e  ",p[i*n+k]);}8 ]4 a: e) |& ^6 C- o
         5 m& B5 ?  N8 T" ~( }' K" `
   printf("\n");}</P>
; l' a' R" Y  Z3 [" x7 r0 u$ p
  R6 M* Q/ i3 a9 `$ O0 T$ \3 H<>
, w. y) y0 T0 p3 @  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
/ b: i/ k' S4 b6 f        for(v=0;v&lt;n;v++): x  y& u# W; q- [$ F% \1 b3 g, ?
  {  c[i*n+v]=0.0;
) N8 P7 m* d! Y$ o6 D. r  C   for(k=0;k&lt;n;k++)
" N8 t2 V) U% B( s5 a    c[i*n+v]+=p[i*n+k]*a[k*n+v];8 w, o; R; ]6 M% u3 a6 a( m' D  q
  }</P>! ^+ v2 d3 x6 B1 g6 ]4 a6 h
<>4 p' i# x& K0 j" C' u/ U
     for(i=0;i&lt;n;i++)& {$ N/ [" n: F* \# _0 e0 q" W
   for(v=0;v&lt;n;v++)7 d- n0 ?) H! r+ l( ?5 m6 x
   {6 l( U" q3 X: |: }, k
    a[i*n+v]=0.0;5 A: q6 ]! V  \$ h$ t- _. F
    for(k=0;k&lt;n;k++)
* q. T6 n, ^: e     a[i*n+v]+=c[i*n+k]*p[k*n+v];
0 i$ [* Z# }% h- [' v   }</P>
( ~2 k3 _  u% ~+ h<>
5 I9 ~* e" H; i$ |  v* \}//最开始的for的结束的大括号
3 ~: B( @' D2 l; P& a return; 2 H8 ?+ x- s, S% W0 S9 n" @: g- o
  }</P>2 w  a2 X! m5 w& x; K6 s; M

3 n4 Y6 r4 X1 F! U<>自己写的运行总是错误</P>
作者: 水木年华zzu    时间: 2009-1-20 23:24
我已经把源代码发在计算数学板块了，你自己找下



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