QQ登录

只需要一步,快速开始

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

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

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

2

主题

0

听众

19

积分

升级  14.74%

该用户从未签到

新人进步奖

跳转到指定楼层
1#
发表于 2005-4-24 22:33 |只看该作者 |倒序浏览
|招呼Ta 关注Ta
<># include "math.h"2 D# w- C4 R, E% E; ]7 q
#include"stdio.h"</P>
  O: p  l  E7 l- j. J<>  void strp(a,p,c, n,u)9 C* }. E, D& S& n" A
      int n;, l" y0 R& h) ]  e+ V' g7 L
   double a[],p[],c[],u[];
7 ?$ P/ F3 v2 |$ T{   : {$ h# c; Y5 M( P: r& U2 Z1 G
   int i,j,k,v;2 E  s, @* O0 L% R* @$ c; N7 o
   double sum,asum;
( m) r7 U* r" e6 D  Zfor(j=0;j&lt;n-2;++j)
( R6 l5 u# [  O' A8 H: I{// 最开始的for 循环
. Y  W2 H# l8 Q       for(i=0;i&lt;n;++i)//初始化u[]全为零) O' J( |0 Q3 ?
     u=0.0;</P>8 M1 |3 A) A! P6 q) o
<>; P5 M7 o2 z) ~* I
      sum=0.0;, ]+ V5 W- ^" x- P: s: _, m- x
   for(i=j+1;i&lt;n;++i)//实现a; `& t0 M) k4 }/ v
      {
# d* S4 e! k* R, a       k=i*n+j;
# m- P+ Q# X/ R    sum+=a[k]*a[k];
3 Z# c/ ]  R. P1 m* t% @( N   }
7 i) f8 ]! [8 [" V" q% E      asum=sqrt(sum);</P>
: g! `3 R) T1 G8 L. F3 J<>      for(i=j+1;i&lt;n;++i)
8 m* ]3 V# s9 n' y" q/ D. W3 E0 R* y   {
9 n5 c% x/ O/ X1 }2 W, H- ?   
: B; {- y- V( F% y6 {$ _    if(i==j+1)
8 O1 m3 E) J$ \) [; m% j4 S8 ?. V     u=a[i*n+j]+(a[i*n+j]&gt;0 ? 1.0:-1.0)*asum;4 Z% B4 L7 z7 U( o6 l1 E) y
    else
4 ~! }3 ?  P) `# m  V9 I8 l     u=a[i*n+j];
! t3 a+ ]( x+ H: I$ ?3 ?   }, r5 n# o1 o+ w6 O
      </P>) T$ X4 \) w* V4 c# c
<>; [# |0 V% r! D% X, v1 o
   sum=0.0;  //实现P
0 q, I4 N% S3 y9 t0 w$ r  D   for(i=0;i&lt;n;++i)
8 x4 ~* W# o: D7 Z  S9 p: E  L    sum+=(u*u);</P>
  f0 b9 g* o# g# z5 T; R# L<>   for(i=0;i&lt;n;++i)
( X" o. p8 ~4 ~   {for(k=0;k&lt;n;++k)8 h7 k+ b9 Q8 x4 U! B# }+ [8 X
    {p[i*n+k]=(i==k?1.0:0.0)+((-2.0)*u*u[k])/sum;
( p& t4 c/ d: s6 u! h4 a    printf("%13.7e  ",p[i*n+k]);}
% j4 W* Q+ ?: e( ?  P" m         4 J  Z& |1 k# P. E$ O
   printf("\n");}</P>: e% v0 q! Q, c

8 ^- ~  {7 h* n5 d' i<>. Z1 N3 g( P8 }  o
  for(i=0;i&lt;n;i++) //实现最后的矩阵相乘/ ^8 u5 _" K+ j) L& u# e
        for(v=0;v&lt;n;v++)
# U* G9 Y0 E) \  V9 Z  {  c[i*n+v]=0.0;
" L- ]2 n( y' N/ }. }4 ?   for(k=0;k&lt;n;k++)
1 T# N4 X( ?9 R1 s    c[i*n+v]+=p[i*n+k]*a[k*n+v];
. Q+ E4 {% @& F0 }2 C  }</P>" ~3 d9 n8 E7 S1 q
<>3 s' n4 S- }) R0 M! ?' }
     for(i=0;i&lt;n;i++)2 i) `9 L+ `7 }6 e" [' [
   for(v=0;v&lt;n;v++)
" S; o5 ~6 o" Y. H. |# X6 Q   {* y3 Z2 H2 H# }% v, z+ z; F
    a[i*n+v]=0.0;
, `8 V7 }6 ]: e1 G) P+ ~' N; n    for(k=0;k&lt;n;k++)
! p$ W6 W5 c6 i# l2 e6 S& k6 I     a[i*n+v]+=c[i*n+k]*p[k*n+v];
; f) ^! U3 P9 k" q$ ]7 D   }</P>
6 Z9 ^! y- U& Z5 i% r( }  W  f<>) c" L& h  l* f) N, `, }5 \9 c3 }' o
}//最开始的for的结束的大括号) j! P- q5 P' A& y( \
return;
. M% ~  ~2 @# B0 L6 F. N3 R  }</P>, ^) n+ Z0 H( `- T0 {$ p& }
* r: S* A% i6 Z) w2 t9 ?- U& e
<>自己写的运行总是错误</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, 2025-8-1 03:33 , Processed in 0.493767 second(s), 63 queries .

    回顶部