QQ登录

只需要一步,快速开始

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

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

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

2

主题

0

听众

19

积分

升级  14.74%

该用户从未签到

新人进步奖

跳转到指定楼层
1#
发表于 2005-4-24 22:33 |只看该作者 |倒序浏览
|招呼Ta 关注Ta
<># include "math.h"
9 L9 W5 w4 ~: C% O& E. J#include"stdio.h"</P>
! j% y. U! A9 `* N<>  void strp(a,p,c, n,u)
. Q( y$ t4 @- L: r      int n;8 s. h3 A; E/ V7 A0 T' e+ R4 |8 u
   double a[],p[],c[],u[];
& I  x$ i( W6 ]{   
! p+ z" a$ u& s" p   int i,j,k,v;
% u' N& x9 z) k, ~   double sum,asum;
4 l/ X* @. `6 W0 v( E0 Ifor(j=0;j&lt;n-2;++j)$ t( N4 G* G3 Y. `4 ^) j
{// 最开始的for 循环
2 q2 q/ Y& k( ^5 x. o       for(i=0;i&lt;n;++i)//初始化u[]全为零
* n" k. n& J5 Q7 u& w     u=0.0;</P>' C5 W0 v. V/ p3 _
<>
) K6 A1 ?) N& J" A: L% n      sum=0.0;/ Y9 U# _2 G9 W- }# f! ]
   for(i=j+1;i&lt;n;++i)//实现a3 O7 B5 A, |" A: }% n+ l/ v
      {
. W+ B1 P1 H* F6 ~/ @% h       k=i*n+j;+ |7 J) D/ T, \9 a2 x
    sum+=a[k]*a[k];
" G- z1 G/ B2 |( y   }
# s! f3 T' G/ X$ k7 a      asum=sqrt(sum);</P>
( D* V1 g4 D$ t<>      for(i=j+1;i&lt;n;++i)8 t/ D& V; p2 [! W, |) M2 F  |6 `
   {
, }# z+ S5 a5 {" r+ O, |# O# s    0 w8 p# l) G0 n- @+ y5 t
    if(i==j+1)
) w8 l. k9 _6 {4 _4 U( q' |5 G* R8 J     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
# O3 s0 b  M. ?! r1 S* b    else
' v  I9 F  J/ \' g1 K: M* A     u=a[i*n+j];3 Y9 m2 f6 L$ }% H* I1 O/ W: u% e
   }
# e# x! H2 y+ z      </P>1 D. e& @* S/ @/ U! d( Z
<>; B, Z5 u2 C: o( E1 o
   sum=0.0;  //实现P
4 A) e# B1 w8 A# s! r   for(i=0;i&lt;n;++i)
3 w0 h  F# g' O    sum+=(u*u);</P>8 Q+ E2 F' n! `' M9 B
<>   for(i=0;i&lt;n;++i)
4 H; A  @) c2 M   {for(k=0;k&lt;n;++k), [7 ~! X1 P# }# D+ T
    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;- F9 F5 z+ a6 K. n+ I
    printf("%13.7e  ",p[i*n+k]);}
8 l' W, C5 h- R  l         
% h$ N7 S! T. Z8 d- ^& \   printf("\n");}</P>
' U! K3 @3 u% j
( ], V9 U& N& t# @<>9 H- T' f8 Z4 i+ s5 T
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘
$ q0 p) T6 f2 ]        for(v=0;v&lt;n;v++)
" _5 f" v2 P* O# v5 {  {  c[i*n+v]=0.0;+ M) F/ h9 H2 @
   for(k=0;k&lt;n;k++)1 ?$ j5 o3 J% O3 K* d/ d
    c[i*n+v]+=p[i*n+k]*a[k*n+v];
  w! X) G3 W% G. ^  }</P>4 G9 f8 s' n" ?- H) ^
<>& m) Y8 k; [4 t; y, p" t. r$ i
     for(i=0;i&lt;n;i++)  G, ?/ z$ o" Q1 |- `/ i
   for(v=0;v&lt;n;v++)2 I, {5 c' u7 G$ E7 p, N
   {( V  S5 N% Y) [1 }* }
    a[i*n+v]=0.0;+ }1 L$ y& k' p0 T3 c* Z
    for(k=0;k&lt;n;k++)3 z9 x+ ^! C4 z9 x# p+ H4 z' A
     a[i*n+v]+=c[i*n+k]*p[k*n+v];
0 y. G  g! y: ~( _; i   }</P>
9 y- \) |: _4 {! S6 a<>
/ X! V, l$ |: u, m5 @7 I- K$ m2 s}//最开始的for的结束的大括号
7 T. j9 t! ~1 }/ N- O return;
1 d4 _! Z  R# D7 E+ f% g7 a  }</P>
- j" T8 u  u: r# a
; A' [# ^' E% E/ O<>自己写的运行总是错误</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-18 08:39 , Processed in 0.555814 second(s), 63 queries .

    回顶部