数学建模社区-数学中国

标题: (求助C语言实现Householder变换一般实矩阵为上Hessenberg矩阵的算法) [打印本页]
作者: aj6249    时间: 2005-4-24 22:33
标题: (求助C语言实现Householder变换一般实矩阵为上Hessenberg矩阵的算法)
<># include "math.h"/ @1 U4 K: S- d+ G) z0 O& b, N
#include"stdio.h"</P>2 w7 |) w+ N# E7 u9 Q  E( X
<>  void strp(a,p,c, n,u)
9 s+ k- L: p6 D7 l6 H  T7 D      int n;; k9 J3 _/ v- L! G
   double a[],p[],c[],u[];
. X4 d- K+ _  A" x{   5 K% a$ d/ D* u- j7 Y
   int i,j,k,v;* \; w" @( l, p6 P/ F/ R
   double sum,asum;
& K3 N: k5 K6 Q! ?' ffor(j=0;j&lt;n-2;++j)+ Q) g9 {. o& \% z$ G* n# O
{// 最开始的for 循环
/ ^0 V$ P0 g4 h9 g  g. L$ d6 g9 D* C- T       for(i=0;i&lt;n;++i)//初始化u[]全为零
* Q8 L- F' r; m$ M1 T) Y     u=0.0;</P>$ k$ M5 _; |8 L" v! M
<>
% |2 X/ }" Z! i, q- P9 Q      sum=0.0;
% l. `9 ?# p- b# A5 W5 e$ \# d   for(i=j+1;i&lt;n;++i)//实现a0 G. |. j2 u" s, u9 l
      {
; Y6 Q4 G! _3 o" Q( {       k=i*n+j;+ H# C7 b+ y* v( T/ j9 l( e
    sum+=a[k]*a[k];
6 ]6 }3 ?2 ]# y  v% M   }! I/ S% f5 h# s2 {; S/ B2 i# f
      asum=sqrt(sum);</P>0 K/ K1 ^/ ?3 J
<>      for(i=j+1;i&lt;n;++i)
6 Q4 y1 b; L6 d# x& W   {) y0 D* g' E& D' x! X1 T! V- x
    , Y$ S. B! D2 F
    if(i==j+1)/ k( {7 m4 {; X8 t" ?9 @
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;1 `4 }0 F* W2 J$ W8 F: F7 `
    else& F- `! d* m- b3 o  _) Q& o% D( f
     u=a[i*n+j];4 R2 c' K' Z' l8 q; i% @! ~
   }
4 i& W1 D" ?- }1 i8 r- Y5 H; H      </P>
4 n$ u. V' W/ G( _<>
8 S5 K* ]0 e* E8 d   sum=0.0;  //实现P
/ |, Y; o( I, e- ^. k" @# m   for(i=0;i&lt;n;++i)
1 J! ^" \1 c8 G: w    sum+=(u*u);</P># V1 u6 ?' v# R& Y8 e
<>   for(i=0;i&lt;n;++i)- F# W3 q3 B0 h0 y' q
   {for(k=0;k&lt;n;++k)' @  ^7 B( f/ m" h1 }( W7 P
    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;) A5 x" C5 G% f1 z- e1 R# y. h
    printf("%13.7e  ",p[i*n+k]);}- Y( R2 C) E% m
         
; N$ ~- ^* _0 E0 K2 g   printf("\n");}</P>
. K- P6 H5 H) f% O7 L; G% U9 u( b5 b; a* b( r
<>
" Y( b- l0 M; ~5 U1 O  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
6 J6 Z5 W6 S' i# U% U! k        for(v=0;v&lt;n;v++)' Z, v6 F3 {- o
  {  c[i*n+v]=0.0;. f/ m" A7 X' s) O8 |5 K
   for(k=0;k&lt;n;k++)
+ S5 [. M: q3 K% {" P    c[i*n+v]+=p[i*n+k]*a[k*n+v];" E% ]7 D# \" m2 Y) Q- ?/ n
  }</P>
# O, G8 d  J" L! }( p<>
! B/ i- y3 l! T  ^8 z4 p, W     for(i=0;i&lt;n;i++)
+ M  ~6 z5 t/ z8 n9 }   for(v=0;v&lt;n;v++)# Z7 [6 O$ j9 S4 R
   {
; n' V0 X9 U6 N$ n    a[i*n+v]=0.0;
4 A/ t- E, X( |: G    for(k=0;k&lt;n;k++)7 Q; W0 N$ J) r
     a[i*n+v]+=c[i*n+k]*p[k*n+v];/ J# M% j$ E/ n3 [6 _# _
   }</P>
, g* h1 e9 j0 a<>, ]7 v2 i4 i) r2 o
}//最开始的for的结束的大括号2 w7 f9 G/ n3 V; }$ Z" |2 P
return; " N2 T9 r; {3 M9 g# }. L
  }</P>
; q# t! x) j$ b0 a) `) t: u
$ `; |0 A5 u+ C8 Y( t1 K* p<>自己写的运行总是错误</P>
作者: 水木年华zzu    时间: 2009-1-20 23:24
我已经把源代码发在计算数学板块了，你自己找下



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