数学建模社区-数学中国

标题: (求助C语言实现Householder变换一般实矩阵为上Hessenberg矩阵的算法) [打印本页]
作者: aj6249    时间: 2005-4-24 22:33
标题: (求助C语言实现Householder变换一般实矩阵为上Hessenberg矩阵的算法)
<># include "math.h"
0 t  ]/ V7 Z6 H/ J! ^# h/ u- m#include"stdio.h"</P>
$ ^" x) f+ r* N, ~<>  void strp(a,p,c, n,u)
% u5 o4 S* g3 O% I      int n;
* n. x" S1 \& x& @   double a[],p[],c[],u[];; V, a) F8 r5 @9 D7 F) b
{   
8 U% ?6 j2 G# j! S5 U$ U   int i,j,k,v;
8 |1 e4 o/ c  f5 i4 l" `; z   double sum,asum;( D$ P* w6 J. L* Z" D- A0 o
for(j=0;j&lt;n-2;++j)( o- V* f' J0 ~+ ^) Y2 n3 n
{// 最开始的for 循环
5 o) i0 B$ I% [  ^3 h       for(i=0;i&lt;n;++i)//初始化u[]全为零
. f* g1 s' o4 C5 s3 z     u=0.0;</P>
) |2 a6 D  W2 ?' L( L1 N<>1 {' m: ^. w* E$ M/ P
      sum=0.0;6 G' Z: N, ?. V- M* d
   for(i=j+1;i&lt;n;++i)//实现a
: ]! q6 T2 o# T3 Y; H      {
7 Q( p9 D+ |. ]       k=i*n+j;
0 G- O; C' e9 Y. r- \! F    sum+=a[k]*a[k];! y! ~/ S( D/ n5 {! w% J4 ^
   }5 A+ E! |/ x. u- J3 N  w2 e4 f
      asum=sqrt(sum);</P>& [" U+ _1 X+ x# ~% ^: F) f" K
<>      for(i=j+1;i&lt;n;++i), f% S/ D) w' D# ~9 C
   {
+ |2 q8 f; S1 a) T9 H' t7 m   
; k( J6 u, q0 V    if(i==j+1)3 i8 W! \" V3 I6 D
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;( E* q" l. v* D* I* i  h, g( H
    else
) g  E* s- s% X; V; }     u=a[i*n+j];% x* Y  p/ M. V# w$ s. O) i. o
   }) U* }5 _1 r' N% C
      </P>
( |# h. U% [. i5 i# f<>% ^( C0 x$ z) G4 I( w
   sum=0.0;  //实现P1 j* x0 X" {: [# f, ?! q
   for(i=0;i&lt;n;++i)
0 d6 X4 N9 i2 ^6 u* o    sum+=(u*u);</P>
5 C3 x; s4 `- X2 g8 L) @/ @5 k<>   for(i=0;i&lt;n;++i)1 _. I9 Z* E- G. K  f
   {for(k=0;k&lt;n;++k)
4 s, f6 d% N. P: Q$ E    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
2 A5 D6 f6 U* d/ |6 z    printf("%13.7e  ",p[i*n+k]);}
) P# p7 s# [) h. }: g         7 K) N1 Q1 O8 W; w3 ]4 {
   printf("\n");}</P>9 c6 c& N2 ~" C& z
* x4 Y  f* h6 ?* C! f: j
<>
. b* w  {9 y6 ]  q  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘/ O5 k  u, q( `
        for(v=0;v&lt;n;v++)
; S. ?3 z7 |8 T1 Y5 d: g8 ?* x" V; P  {  c[i*n+v]=0.0;
2 k# `- Y. ?& F; M2 g   for(k=0;k&lt;n;k++). L4 ]7 P( \% D& Q2 }# B
    c[i*n+v]+=p[i*n+k]*a[k*n+v];
6 t6 v  v4 x" ]: d( F  }</P>6 f% M8 w% L$ }3 v4 Q7 W
<>
+ G: A# b+ g( j     for(i=0;i&lt;n;i++)) m2 x6 ]" N3 M
   for(v=0;v&lt;n;v++)% o/ O/ Y; f* C) R. d8 @( z2 J
   {
3 r  v! y8 H* ]    a[i*n+v]=0.0;
  J# Z" G, L( W+ g) z$ `& ^& V+ Q8 x    for(k=0;k&lt;n;k++)
4 n+ P/ X3 a; I& i/ G; O# A8 k     a[i*n+v]+=c[i*n+k]*p[k*n+v];$ _- L4 J  X: q/ W" X3 V
   }</P>
) u; i( K; \* g6 ]<>
+ H* @) h0 d4 I8 m. A" W}//最开始的for的结束的大括号
2 S% z7 \3 J  O4 r+ s return;
7 Y3 `! j5 x3 R/ o3 X" r  }</P>7 ~! z& a3 P" h7 o; L1 o  P9 g

( J6 c/ A2 ^& A& G/ X<>自己写的运行总是错误</P>
作者: 水木年华zzu    时间: 2009-1-20 23:24
我已经把源代码发在计算数学板块了，你自己找下



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