QQ登录

只需要一步,快速开始

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

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

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

2

主题

0

听众

19

积分

升级  14.74%

该用户从未签到

新人进步奖

跳转到指定楼层
1#
发表于 2005-4-24 22:33 |只看该作者 |倒序浏览
|招呼Ta 关注Ta
<># include "math.h"( m' g$ D1 s/ r  I% E8 j4 b
#include"stdio.h"</P>8 o/ C9 g; N6 V# {
<>  void strp(a,p,c, n,u)
* S3 o8 n% f7 ]3 q' \  Q) W0 ?1 K      int n;
# o0 U/ r* ]: y   double a[],p[],c[],u[];+ u4 v1 e) u" H; |8 F* D' Z
{   ( T3 A3 _3 g! k$ A% }
   int i,j,k,v;
+ K: U2 }  @& y( i/ V) m$ Y   double sum,asum;
- @# Q2 T9 ~9 l$ A; B8 |for(j=0;j&lt;n-2;++j)4 I. V5 i. V; i* u6 G8 j. d0 b& Z
{// 最开始的for 循环
2 p% \. `3 L% v* y5 x! }8 {6 A3 W       for(i=0;i&lt;n;++i)//初始化u[]全为零
8 D0 F; N& S5 C; X% `; E, V     u=0.0;</P>
" d1 j; c" d7 m$ g/ ?<>- ~0 s) ~5 m: _* s* p9 A
      sum=0.0;/ ]; K8 L# P# U2 f% l) h
   for(i=j+1;i&lt;n;++i)//实现a
0 Z1 B' _+ l5 [  L      {
) R: S2 v- c. g/ L       k=i*n+j;$ N2 h5 g  \2 N6 ?9 L7 O( d
    sum+=a[k]*a[k];2 @, i% s9 K% {8 X' S
   }
. }- C5 W) H" C, t' P      asum=sqrt(sum);</P>
$ D! y/ d3 d) A/ @* s& S; g  c2 j3 R) |7 ?# n<>      for(i=j+1;i&lt;n;++i)
! m3 v" b8 \2 `. x) d+ L   {& j9 U( \3 P3 _- S
   
6 V. ^: c1 ^! }% _& ?- Q1 o    if(i==j+1)+ ?% f+ A) H; d3 U
     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;
! l9 ^/ b" X. B/ u8 g4 |/ B    else: O3 V2 q( H5 R9 r' Q
     u=a[i*n+j];& c+ u8 ?7 ?/ q6 s, C
   }
1 c1 S. C- Z( a2 I9 W/ E      </P>
' O5 j( y6 O% [4 M% B. ^<>
3 d$ J( _2 V5 f. ?   sum=0.0;  //实现P
& x. G2 S' @/ s" y6 w) K& U   for(i=0;i&lt;n;++i)
# q$ S5 D" S2 }2 ?) u    sum+=(u*u);</P>, T3 D0 E2 d- d" l) u. G2 d5 h) n
<>   for(i=0;i&lt;n;++i)
6 h9 Q  t- r, b, w! l) a2 N0 k" H" c   {for(k=0;k&lt;n;++k)2 j. j0 L2 N2 W7 U) {
    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
8 |9 K7 U- l# m# Z- T( _3 w    printf("%13.7e  ",p[i*n+k]);}
- h) l# `$ P  z+ h5 }& o         $ T, S* H  R1 y, e9 ?
   printf("\n");}</P># c- U) o2 F9 k  U! l

1 ]8 D5 x- s5 l  q<>" j5 s! H, K5 N1 h  R
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘2 V, t* _1 k- X
        for(v=0;v&lt;n;v++)
# q$ X$ c3 G. u$ c' M" J- E  {  c[i*n+v]=0.0;
# M6 I  E3 O. {, U   for(k=0;k&lt;n;k++)
+ R. x7 g7 f( c) @( \    c[i*n+v]+=p[i*n+k]*a[k*n+v];) A8 }( {) F, Q; x
  }</P>  q8 T9 V% H5 }! t
<>
) U) C# \. L5 Q7 ~) J# K     for(i=0;i&lt;n;i++)$ N- E  [. x3 `5 Q
   for(v=0;v&lt;n;v++)
4 x: h7 i) v, n: }- c   {7 W' H% l% K% w
    a[i*n+v]=0.0;! Y" d0 {, x1 M5 G$ o0 M4 x
    for(k=0;k&lt;n;k++)" }6 ]4 l0 `) L; T3 f; O. p0 i
     a[i*n+v]+=c[i*n+k]*p[k*n+v];
! \0 K3 e- t. A   }</P>6 }3 u0 @& F/ \4 P! v# |
<>
7 w- a& D2 j; s+ R" s  v}//最开始的for的结束的大括号
3 M! o0 X! h0 R5 M1 ]. ]' ? return;
" t2 ]& x4 f/ v1 }9 |  }</P>
  t$ h# |: X& ~$ T1 h2 m/ X* d, ?3 q( J* y- n" H
<>自己写的运行总是错误</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 07:32 , Processed in 0.310731 second(s), 63 queries .

    回顶部