QQ登录

只需要一步,快速开始

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

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

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

2

主题

0

听众

19

积分

升级  14.74%

该用户从未签到

新人进步奖

跳转到指定楼层
1#
发表于 2005-4-24 22:33 |只看该作者 |倒序浏览
|招呼Ta 关注Ta
<># include "math.h"$ u. b5 w* G+ n5 k) [, _% Z
#include"stdio.h"</P>
6 D$ ?' Y# N3 ]1 P<>  void strp(a,p,c, n,u)
' G! A- [; |+ o/ O6 ~      int n;
" B) c9 U& K; u' p0 g) ?   double a[],p[],c[],u[];
. K6 {0 t) f# A# \6 [) G2 `! s{   
1 H( u8 W9 K* p1 e  w   int i,j,k,v;
3 \6 n. F% J, B" f3 [   double sum,asum;! l4 N3 a/ ?; a1 q0 Y0 B5 U
for(j=0;j&lt;n-2;++j)
  V7 _1 f0 E! l9 A& O{// 最开始的for 循环
, Y2 r8 @8 \' d% P% X       for(i=0;i&lt;n;++i)//初始化u[]全为零) U; j4 J& ^1 I' \5 R) E3 o
     u=0.0;</P>
/ |2 ~/ R1 p  |. P, t" S<>! C( s$ I  z3 H
      sum=0.0;' P7 h# H* m8 C+ ?% q, ^
   for(i=j+1;i&lt;n;++i)//实现a, M- ]2 f* |/ X/ Z
      {
# y: d5 l% ~7 C' s, E$ G& i2 M       k=i*n+j;) w" z5 F" m1 X$ K
    sum+=a[k]*a[k];
" ?" m: b, q; W( t   }
# {1 f' X& Z0 c& k7 L      asum=sqrt(sum);</P>$ T: u$ y4 a0 \/ _9 B
<>      for(i=j+1;i&lt;n;++i)
. E" l/ r3 k- X" _3 x7 H& Q   {
# l; D! O9 i' Q" F1 F5 P+ a0 j    - p; Y! t5 `5 D: u) c
    if(i==j+1)
: O% k. ~* f1 x6 B     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
! u. ^0 D, ^; f" l8 A" J" T0 N    else
$ z4 F. D2 D+ C# S; ]) V* K' q9 K     u=a[i*n+j];% a2 L( z! L" V/ G
   }
3 z$ E% X* X! f, q0 p6 g      </P>
: s  k; h0 [. K, A2 q<>
% r+ T* M$ I0 C   sum=0.0;  //实现P) Z3 F- j1 Y/ t0 r* V8 C+ m. r
   for(i=0;i&lt;n;++i)
( x* j) N0 t4 W6 e: a+ j/ V: ]. n    sum+=(u*u);</P>: P; R% M' a6 i, u: X
<>   for(i=0;i&lt;n;++i)* ~; ]3 \+ K/ t1 g  L. {
   {for(k=0;k&lt;n;++k), m" l0 G0 w* E+ j, U7 t
    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
7 ^4 T0 z& V3 T6 K2 e: \+ K  ^    printf("%13.7e  ",p[i*n+k]);}( p$ T5 K* B5 n. \9 _) }
         3 f0 d# |  d# H! w, K: K# j8 ?6 ~
   printf("\n");}</P>
, Z0 n6 e# r9 n" }" P3 N/ s1 j
/ J6 B% n2 x# T. t1 O4 ]( G$ n1 s, o<>; w6 q- g9 K/ p% v' A
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘/ Q5 o8 _; W) o4 o4 w( w
        for(v=0;v&lt;n;v++)" c3 W# G- V- ^7 g% L) L# P9 V
  {  c[i*n+v]=0.0;
, C8 A8 z$ N9 X) j5 B   for(k=0;k&lt;n;k++)
& T4 n2 c) d. ~( u: b  y. I    c[i*n+v]+=p[i*n+k]*a[k*n+v];
. ^8 Q" @/ p8 k  v  C! o$ V2 _  }</P>. A  ?0 Q( N2 A2 F0 w
<>
0 n$ ]0 y" a7 _, r     for(i=0;i&lt;n;i++)- d6 `/ x4 }: F! c( c9 m% b8 P+ w
   for(v=0;v&lt;n;v++): S; l; Y5 n( R* o( p9 }) G
   {
1 U- v& x# R" n; L9 z/ v. _% R    a[i*n+v]=0.0;
7 D* S+ i% M9 z& l    for(k=0;k&lt;n;k++)  w: C3 w% D" S
     a[i*n+v]+=c[i*n+k]*p[k*n+v];; {3 F5 u( q: p$ [* Q
   }</P>
1 w4 D( ^( f+ u. ^% Y<>- p: m0 v, L$ t' Y3 I1 N
}//最开始的for的结束的大括号
" s4 `1 H+ q4 ^ return;
% b/ H- n: i" p7 v; U  c& S" c5 g3 d  }</P>8 B. `+ o% o' u; J: s. F; \' ?

  j" m6 ?( Y* ~' `. `$ t+ z<>自己写的运行总是错误</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-4-17 07:17 , Processed in 0.392065 second(s), 62 queries .

    回顶部