数学建模社区-数学中国

标题: (求助C语言实现Householder变换一般实矩阵为上Hessenberg矩阵的算法) [打印本页]
作者: aj6249    时间: 2005-4-24 22:33
标题: (求助C语言实现Householder变换一般实矩阵为上Hessenberg矩阵的算法)
<># include "math.h"
- W  g: _2 E' M( I1 ^3 u( x  O+ J#include"stdio.h"</P>
/ ^) ?- o& q( ~6 c% ]. p# o<>  void strp(a,p,c, n,u)
+ J) T, Z4 J& l6 t2 Z1 \- p      int n;
2 L. l7 }) x$ V8 H6 Z* [- I" C7 T5 f   double a[],p[],c[],u[];
; E/ J& s& C: N5 |4 m{   ; ?/ N, B7 O5 R" y. @
   int i,j,k,v;' C, N# o+ I) x- q# S7 j0 j
   double sum,asum;6 H3 t# P! c) j0 s2 @  i
for(j=0;j&lt;n-2;++j)
" s  ~) ^0 w5 ~" I  E{// 最开始的for 循环
& M* q; I8 e% ?3 S0 \* Y       for(i=0;i&lt;n;++i)//初始化u[]全为零
+ w0 s2 e$ a! D8 F     u=0.0;</P>
- z5 i) P$ v5 |. Q; A5 z6 i6 c<>/ T# w# i1 V; k2 I3 i
      sum=0.0;- h3 N% a1 u" [) }
   for(i=j+1;i&lt;n;++i)//实现a
9 d& K, \0 m4 X' z      {
) y6 O) c# M: B- ?- c       k=i*n+j;& p! M' v. w$ a. ]6 T
    sum+=a[k]*a[k];
# B& S% l( J0 t$ p' f   }1 }8 r7 C6 q+ D1 s
      asum=sqrt(sum);</P>
' _/ l9 f; w; }2 Q$ ~5 u<>      for(i=j+1;i&lt;n;++i)# D, d2 f1 i1 d6 i, b
   {
+ H* l. e: `. A2 I   
1 w: h6 x3 S" A; O; o    if(i==j+1)) F$ s9 u& T. H) m$ y1 M6 N0 g' `
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;8 p; V1 q, [9 k- ]
    else
6 q; }7 j( |! Y  j" D     u=a[i*n+j];6 [. Q" i  r+ m5 V! i
   }% v6 _. m% X$ [+ I
      </P>; t9 w6 l4 h: ]/ p6 d
<>% c7 @. F: R4 o& D% A! h- e0 M8 H1 G
   sum=0.0;  //实现P% R& J/ j# [. K
   for(i=0;i&lt;n;++i)6 V% o5 g7 R/ v2 Z3 h
    sum+=(u*u);</P>' x7 a& n8 w( s, c) C
<>   for(i=0;i&lt;n;++i)
% i/ A) s1 U" o   {for(k=0;k&lt;n;++k)
. T! C" m9 u9 V6 {+ R: j) _    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
1 }0 N7 t7 w) R6 X" A% G    printf("%13.7e  ",p[i*n+k]);}; o5 S$ [; {+ m6 U. ~% d3 m: i  R: g
         6 i) B+ _( B8 m
   printf("\n");}</P>/ K( V: r6 b. |0 W
* }9 f2 Y* R9 I2 w" `, a, r
<>
6 `* f9 M, O) d, l  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
  q! M2 E' _- f" Z" K7 g( B4 Y' t6 O* {        for(v=0;v&lt;n;v++)
- C- I$ ~3 c5 n. E% y  {  c[i*n+v]=0.0;
, _+ E! W5 f- W1 r# ~   for(k=0;k&lt;n;k++). N& c: v- _  E: ^
    c[i*n+v]+=p[i*n+k]*a[k*n+v];
" w( j# t( Q5 k4 l. l  }</P>2 G' l1 G! \' D8 \9 H6 u( ?
<>' T: }$ ^' ?) j. h* D8 Y4 ?- k( D
     for(i=0;i&lt;n;i++)
4 q% E& C* A, e" o" p3 |   for(v=0;v&lt;n;v++)$ ^8 O3 ~6 p' F- M6 p! d( }
   {. S$ W  O9 _$ ^2 X- y
    a[i*n+v]=0.0;
: C' V4 q! D; d    for(k=0;k&lt;n;k++)! p) X( Y5 d- X, S% o: `- \
     a[i*n+v]+=c[i*n+k]*p[k*n+v];
, F! A' ~- \* z& \   }</P>
: u. K0 }  \6 E. E$ R<>. ?: _+ b% y8 F$ P# ]8 t8 d6 o
}//最开始的for的结束的大括号  @( k/ a3 _  [8 Z; D- ]
return;
! @/ f4 k( g5 g9 l9 N  }</P>: T$ z! I9 Y- \; R/ |# L

/ }& \5 _8 a, a: c# f+ l<>自己写的运行总是错误</P>
作者: 水木年华zzu    时间: 2009-1-20 23:24
我已经把源代码发在计算数学板块了，你自己找下



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