QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 6308|回复: 1
打印 上一主题 下一主题

(求助C语言实现Householder变换一般实矩阵为上Hessenberg矩阵的算法)

[复制链接]
字体大小: 正常 放大
aj6249        

2

主题

0

听众

19

积分

升级  14.74%

该用户从未签到

新人进步奖

跳转到指定楼层
1#
发表于 2005-4-24 22:33 |只看该作者 |倒序浏览
|招呼Ta 关注Ta
<># include "math.h"
1 l, l5 \7 _8 g  G# m#include"stdio.h"</P>* O; ?. Q9 o" P
<>  void strp(a,p,c, n,u)1 \- M0 E1 O8 @: l, g7 J8 H
      int n;6 L1 h9 h; z) [0 O
   double a[],p[],c[],u[];! ^' S& y" s) }# a- d
{   ; I# D( T" n' O1 M
   int i,j,k,v;1 S5 C0 H4 `2 t1 V! x" a
   double sum,asum;" \1 [+ L) j' F* s' G" Q$ S% W: P9 _! C
for(j=0;j&lt;n-2;++j)
( d1 w9 @1 ^! }5 X- o# \5 v: ^{// 最开始的for 循环
: N: C# U! u" L' C  z6 Z       for(i=0;i&lt;n;++i)//初始化u[]全为零! x8 [" K6 i8 F) k/ M
     u=0.0;</P>
. ^! a& h6 e1 Q<>) j5 A# m( t% n3 l5 h3 ]4 }
      sum=0.0;
8 {* f8 \/ v4 V) _   for(i=j+1;i&lt;n;++i)//实现a
6 T" A* [8 Y1 T      {
% z0 L: t) ^7 u" t$ \, W) E       k=i*n+j;
7 H+ v) C1 R' Y5 O, \    sum+=a[k]*a[k];
. T0 m& v% H; G9 ^- c' ?; S; S   }
/ R# K  M" D' H( f      asum=sqrt(sum);</P>
2 j' r4 ^) m5 }2 i7 k3 C6 l<>      for(i=j+1;i&lt;n;++i)
. B, v! O6 \6 S) x( J% X   {& A" K( G& ^6 t0 c
   
. m$ {1 O( v, m5 z    if(i==j+1)% `( a7 s2 M' I% j, u1 C+ I, f" L2 i
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
+ \" t: H/ Z! j' X2 {7 Q; l    else, ~" l1 O/ }  D, t7 D
     u=a[i*n+j];
6 d# Z1 ]8 [+ I/ G   }
2 K& k4 W: S6 A1 A1 E$ ^$ g' X3 t      </P>; Z) }! ~' |" q4 s3 w1 T
<>
. W( E+ `% ?1 l& ^; p   sum=0.0;  //实现P! q. ^1 C; c, \( a$ i! @
   for(i=0;i&lt;n;++i)
5 Z/ ?* Q) E% C    sum+=(u*u);</P>7 M- F/ c* O  w3 I
<>   for(i=0;i&lt;n;++i). u  V5 A  J5 M6 x5 Y& |
   {for(k=0;k&lt;n;++k)
' y: s0 j/ i$ f    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;1 n% ~; ]7 J( Z$ C; N5 W% \: I* L
    printf("%13.7e  ",p[i*n+k]);}
8 \1 X; E* x- \4 M; n) ]6 \* t         + U+ j, Z) n6 k4 F4 h0 q* H  H
   printf("\n");}</P>  ?0 g4 k3 ~1 ]! b. s& ?2 s
( |! U: G$ o% Z6 a  {9 i
<>: `( G3 [- J7 C6 \2 i7 {" i- [
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘' y3 Z* `: b) w; u4 ]
        for(v=0;v&lt;n;v++)2 L+ }9 Q; P9 X0 K" i0 |8 Y
  {  c[i*n+v]=0.0;
" E% B1 j& S- {   for(k=0;k&lt;n;k++); y; U7 S/ K, M& L# {, u" c
    c[i*n+v]+=p[i*n+k]*a[k*n+v];
3 `  t0 P" v' [$ _  }</P>
% w, z" q+ R2 H$ k1 @<>
+ f+ U  d+ O" [( O' U     for(i=0;i&lt;n;i++)8 m% Q" |) N6 I0 I  n( r' @( X
   for(v=0;v&lt;n;v++)9 }, c6 H1 A" ~5 I
   {2 i& O( P# ]+ {
    a[i*n+v]=0.0;
- w# N1 H; L0 x6 J1 M2 L1 w  o    for(k=0;k&lt;n;k++)* ~: |* c) o( v) y6 u2 H4 h2 g
     a[i*n+v]+=c[i*n+k]*p[k*n+v];
1 c. R) y7 d& F. N   }</P>$ J4 J5 U4 h# x
<>6 R* V9 P! J# P# {
}//最开始的for的结束的大括号# n3 K/ B% Q. `' o) n5 B; ?
return;
. B# ?2 Q2 x$ H' t  }</P>
% U8 q+ o9 i- C; r2 u9 X
" R& L; j/ s* m" T2 z) r<>自己写的运行总是错误</P>
zan
转播转播0 分享淘帖0 分享分享0 收藏收藏0 支持支持0 反对反对0 微信微信

21

主题

7

听众

3435

积分

升级  47.83%

  • TA的每日心情

    2014-5-25 20:58
  • 签到天数: 20 天

    [LV.4]偶尔看看III

    新人进步奖 优秀斑竹奖

    群组Matlab讨论组

    群组小草的客厅

    群组数学趣味、游戏、IQ等

    群组C 语言讨论组

    群组我行我数

    回复

    使用道具 举报

    您需要登录后才可以回帖 登录 | 注册地址

    qq
    收缩
    • 电话咨询

    • 04714969085
    fastpost

    关于我们| 联系我们| 诚征英才| 对外合作| 产品服务| QQ

    手机版|Archiver| |繁體中文 手机客户端  

    蒙公网安备 15010502000194号

    Powered by Discuz! X2.5   © 2001-2013 数学建模网-数学中国 ( 蒙ICP备14002410号-3 蒙BBS备-0002号 )     论坛法律顾问:王兆丰

    GMT+8, 2026-9-3 06:37 , Processed in 0.280473 second(s), 63 queries .

    回顶部