QQ登录

只需要一步,快速开始

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

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

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

2

主题

0

听众

19

积分

升级  14.74%

该用户从未签到

新人进步奖

跳转到指定楼层
1#
发表于 2005-4-24 22:33 |只看该作者 |倒序浏览
|招呼Ta 关注Ta
<># include "math.h"
- f- P4 K2 A7 R* g( ~& K7 [#include"stdio.h"</P>
3 u; {" U7 |4 `# k7 m( g+ X<>  void strp(a,p,c, n,u)
& O% ?3 f+ F  x      int n;
) I+ H: e6 @2 K% z0 p   double a[],p[],c[],u[];2 I; x: q$ \  D. ^0 j
{   " \% c% s; u* l6 V- E8 Z2 s- K
   int i,j,k,v;
3 _& m9 K- f  [& \/ d; W* S   double sum,asum;" E6 V5 N; D8 _& w5 L, M/ }8 D# G& \, y
for(j=0;j&lt;n-2;++j)! ]* u1 i( h: N' ^
{// 最开始的for 循环
: x! U' `5 K  K/ T  k       for(i=0;i&lt;n;++i)//初始化u[]全为零
% s4 z0 Q, W& K. S+ K     u=0.0;</P>4 V) M" e! \. D/ d# X
<>3 R4 Z$ E* W+ p0 Z
      sum=0.0;' R& u* p( K0 \
   for(i=j+1;i&lt;n;++i)//实现a
9 w3 V$ z: ]3 I# F1 d      {) L; D, E) y- b0 `
       k=i*n+j;
0 t( g: q% m7 {. k" U    sum+=a[k]*a[k];
) Z3 u3 b/ g# ?5 Y. B1 A/ G; p   }
4 T; Z( |7 T( V% [      asum=sqrt(sum);</P>
8 k0 e; k' Z1 P4 x/ ^: P2 A  B<>      for(i=j+1;i&lt;n;++i)2 K4 m7 G( m0 L" l3 k3 M4 L& Q: t
   {! ]$ ~/ F+ U( A5 `' `
    6 b' ?! m5 E7 \8 @- E  i
    if(i==j+1)
8 p% E9 ]  |$ h& s# U# A     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
+ B& h  F8 a) X1 ?& }$ _5 y    else
! X0 ]5 l9 V/ j3 ^( V& w8 V+ K8 O- r. r     u=a[i*n+j];6 v; T9 w$ N; p/ s2 z% Z
   }
( w. u1 m! B, j! u1 m0 f# Q      </P>
+ O/ o1 |7 }$ Q4 ?% V3 E' @8 s6 I<>
+ u+ s+ N6 Y) U( G   sum=0.0;  //实现P
* \0 a7 k. E( j9 ]# O   for(i=0;i&lt;n;++i)
* x0 A  z* R: X# }3 T0 Y% M% x    sum+=(u*u);</P>' A$ o  E# |& E1 E
<>   for(i=0;i&lt;n;++i): k: |7 H  w. d2 y
   {for(k=0;k&lt;n;++k)
3 p1 v' u8 s- x/ G  [' U" c7 k5 l    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;2 _( s- S' c5 G
    printf("%13.7e  ",p[i*n+k]);}3 v& G7 X& J  N9 [8 \
         5 P9 H3 k7 f" K
   printf("\n");}</P>
  @; e" }$ q- ?- u/ t; M1 }( _, K% W$ J, W" z
<>
# |6 X6 r" r1 Z7 F, D  V5 q+ p3 S  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘7 t0 N4 y3 I, Y2 q4 n
        for(v=0;v&lt;n;v++)
9 W# }1 p, k7 K' Y5 w  {  c[i*n+v]=0.0;4 |9 {& |" `8 C. d
   for(k=0;k&lt;n;k++); K) b, O) N5 k. ]
    c[i*n+v]+=p[i*n+k]*a[k*n+v];
  E% t& N& C" j/ ^  }</P>  h6 {' i" f$ [) w0 I" {7 D8 e
<>9 H8 L+ B" M) ]6 P7 y3 X
     for(i=0;i&lt;n;i++)/ v7 m& t( S' B7 h$ G
   for(v=0;v&lt;n;v++)
5 U! v# ^8 m% L7 g   {* W, k  o7 u$ T  ?
    a[i*n+v]=0.0;
' ?  z# D  [$ t. e; D& a    for(k=0;k&lt;n;k++)
5 R. ~4 L) J( a5 [; I: T     a[i*n+v]+=c[i*n+k]*p[k*n+v];, v8 M9 s0 `% Q3 L
   }</P>* F( j! c+ z  P6 f1 d- w) W
<>
2 j1 A8 L% y; Y6 k# n8 A}//最开始的for的结束的大括号
3 N) V3 q# g, ^; L  s return;   X, ]5 @0 Q3 h$ r
  }</P>
9 z* G* I5 _+ D8 B! l# b* g3 b; X, k2 V, }/ z+ F9 F
<>自己写的运行总是错误</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-7-19 16:00 , Processed in 0.439142 second(s), 63 queries .

    回顶部