应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » Fortran 求解复数矩阵SVD时，用到cgesvd，结果与matlab不同，为什么？
161/2返回列表12下一页
查看: 2856  |  回复: 15
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

长路漫漫

新虫 (初入文坛)

[求助] Fortran 求解复数矩阵SVD时，用到cgesvd，结果与matlab不同，为什么？

RT，很小的一个测试矩阵，代码如下
      integer M,N,MN,Mmax, Nmax,p, nu, nv
         parameter  (M=6,N=4,MN=min(M,N))
          complex :: A(M,N),AA(M,N)
          character , parameter :: job = 'A'
          real                    S(MN)
        complex  :: U(M,M), VT(N,N), V(N,N)
        complex  :: work(1000)
        real  :: rwork(1000)
          integer    info,rw,cl,lwork,kk,i,j,k
         do rw=1,M
             do cl=1,N
                 A(rw,cl)=cmplx(rw,2*cl)
             enddo
          enddo
          lwork=MAX(1,2*MIN(M,N)+MAX(M,N))
       call CGESVD( JOB, JOB, M, N, A, M, S, U, M, VT, N,   &
      &                  WORK, LWORK, RWORK, INFO )
        V=transpose(conjg(VT))  
        aa=u*s*v'
得到的u,v只有前2 列与matlab相同，而后面几列的差距很大，S的结果两边差异非常小，基本一致。
而矩阵aa与原始的一样。有大神可告知原因么？

[ Last edited by 长路漫漫 on 2013-6-26 at 11:09 ]
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2013-06-26 11:02:21
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
回帖支持 ( 显示支持度最高的前 50 名 )

somomo91

专家顾问 (职业作家)

【答案】应助回帖

★ ★ ★ ★
感谢参与，应助指数 +1
jjdg: 金币+1, 感谢参与 2013-06-27 00:31:32
长路漫漫: 金币+3, ★★★很有帮助 2013-06-27 10:22:44
这个问题几乎是公认的
MATLAB 虽然调用 CGESVD，但是它对输入参数进行了预处理，
并且因为编译器的不同，数据精度处理也不大相同，
因此导致两个结果不大相同，
赞一下(1人) 回复此楼
2楼2013-06-26 21:02:44
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通回帖

somomo91

专家顾问 (职业作家)

【答案】应助回帖

★
jjdg: 金币+1, 感谢参与 2013-06-27 00:31:39
具体的讨论，见LAPACK官方论坛：
http://icl.cs.utk.edu/lapack-for ... ?f=2&t=1689

当然还有其他地方，都说过类似的问题的
赞一下 回复此楼
3楼2013-06-26 21:03:55
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

长路漫漫

新虫 (初入文坛)
引用回帖:
3楼: Originally posted by somomo91 at 2013-06-26 21:03:55
具体的讨论，见LAPACK官方论坛：
http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=1689

当然还有其他地方，都说过类似的问题的

谢谢，非常有帮助，只是，我的矩阵很小，而且对精度要求也没有那么高，给你看看一个4×3的复数矩阵得到的两种结果
矩阵A
M=4;
N=3;
for ii=1:M
for jj=1:N
a(ii,jj)=ii+i*2*jj;
end
end
matlab：

u =

  Columns 1 through 4

  -0.2028 - 0.2666i   0.5235 + 0.3708i   0.4803 + 0.0812i   0.3898 - 0.1821i
  -0.2559 - 0.2442i   0.3188 + 0.2616i  -0.4234 - 0.1552i  -0.3399 + 0.2662i
  -0.3091 - 0.2217i   0.1142 + 0.1524i  -0.3466 + 0.1235i   0.0841 - 0.0690i
  -0.3622 - 0.1993i  -0.0905 + 0.0432i   0.3628 - 0.0013i  -0.5831 + 0.1887i
  -0.4153 - 0.1768i  -0.2952 - 0.0660i  -0.3939 - 0.1528i   0.3245 - 0.3249i
  -0.4684 - 0.1544i  -0.4998 - 0.1753i   0.3208 + 0.1047i   0.1246 + 0.1211i

  Columns 5 through 6

   0.0100 - 0.1874i  -0.0737 + 0.1148i
   0.0538 - 0.3784i   0.2691 - 0.3108i
  -0.1977 + 0.7171i   0.2468 + 0.2400i
   0.0048 + 0.2312i  -0.5080 - 0.0163i
   0.3183 - 0.1634i  -0.4324 - 0.0181i
  -0.1892 - 0.2192i   0.4982 - 0.0096i

s =

   32.8054         0         0         0
         0    2.7938         0         0
         0         0    0.0000         0
         0         0         0    0.0000
         0         0         0         0
         0         0         0         0


v =

  -0.3202                      -0.7730             0.2225             0.5005         
  -0.3972 + 0.1228i  -0.3530 - 0.0509i  -0.1119 - 0.3241i  -0.7494 + 0.1441i
  -0.4742 + 0.2455i   0.0671 - 0.1017i  -0.4436 + 0.6482i  -0.0026 - 0.2882i
  -0.5512 + 0.3683i   0.4871 - 0.1526i   0.3331 - 0.3241i   0.2515 + 0.1441i

zegsvd结果：

u=
(-0.202824375016695,-0.266642178187997) (-0.255949364969120,-0.244187270428356)
(-0.309074354921545,-0.221732362668715) (-0.362199344873970,-0.199277454909075)
(-0.415324334826394,-0.176822547149434) (-0.468449324778819,-0.154367639389794)

(0.523495728897005,0.370844937064488) (0.318833666239141,0.261621268331258)
(0.114171603581277,0.152397599598028)
(-9.049045907658722E-002,4.317393086479759E-002)
(-0.295152521734451,-6.604973786843274E-002)
(-0.499814584392315,-0.175273406601663)

(-0.340495458990843,0.541393498514851) (0.337828171542821,-0.540655966822066)
(0.160806639553680,-0.275424645436289)
(4.866929136520519E-002,7.466583397862278E-002)
(-0.228593892608515,0.132598643065634)
(2.178524913765310E-002,6.742263669924775E-002)

(-0.213595095465232,1.770559300328978E-002)
(-8.118515338632207E-002,-0.101568035512989)
(0.142978695055282,-2.696899503189933E-002)
(0.159236069089761,0.159954322471254)
(0.645307867526083,7.874251718869874E-002)
(-0.652742382819571,-0.127865402118354)

(0.101959655905054,7.515513816888129E-002)
(0.385377896491782,-6.292538458577177E-002)
(-0.389645765970497,-8.452303330092807E-002)
(-0.660124183763325,0.118679170354774) (0.437875799945744,-0.107863393308084)
(0.124556597391242,6.147750267112817E-002)

(-5.714988270508767E-002,4.358750066103587E-002)
(0.164737090936639,-0.319067431094656) (-0.215932231409642,0.700403583527349)
(9.841764078738405E-002,-0.548295272551258)
(7.776246243304019E-002,5.371201559391532E-002)
(-6.783508004233359E-002,6.965960386361419E-002)



s=
32.8054068874654     
   2.79379293216862     
  9.854419884063358E-016
  6.434838164926571E-016


v=
(-0.320194781107884,0.000000000000000E+000)
(-0.397196078916262,0.122773730945919) (-0.474197376724640,0.245547461891838)
(-0.551198674533019,0.368321192837758)

(-0.772965265811650,0.000000000000000E+000)
(-0.352952325959982,-5.085805228873791E-002)
(6.706061389168552E-002,-0.101716104577475)
(0.487073553743353,-0.152574156866213)

(-0.282895611404907,0.000000000000000E+000)
(0.402086157056524,-0.348691818157092)
(4.451452010167405E-002,0.697383636314183)
(-0.163705065753291,-0.348691818157091)

(-0.469009672659151,0.000000000000000E+000)
(0.610331955712300,0.210322709402046) (0.186365106552852,-0.420645418804091)
(-0.327687389606002,0.210322709402045)

很明显可发现s基本一致，而u和s从第三列开始就区别非常大，跟你
具体的讨论，见LAPACK官方论坛：
http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=1689
讨论的内容还挺不一样的，我之前实验过DGESVD这个，处理实数矩阵时，对于一般精度的矩阵，U,S,V结果基本差不多，误差非常小了，而这个CGESVD和ZESVD处理复数矩阵，误差非常大，就LAPACK官方论坛上，也没有搜索到与CGESVD和ZESVD处理复数矩阵的类似帖子
为这个问题纠结了很久，麻烦大神进一步解答哈，谢谢～～
赞一下 回复此楼
4楼2013-06-27 11:20:53
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

长路漫漫

新虫 (初入文坛)
引用回帖:
4楼: Originally posted by 长路漫漫 at 2013-06-27 11:20:53
谢谢，非常有帮助，只是，我的矩阵很小，而且对精度要求也没有那么高，给你看看一个4×3的复数矩阵得到的两种结果
矩阵A
M=4;
N=3;
for ii=1:M
for jj=1:N
a(ii,jj)=ii+i*2*jj;
end
end
matlab：

u = ...

这里有些写错了，现在进行修正
做的矩阵式6×4的，
M=6;
N=4;
for ii=1:M
for jj=1:N
a(ii,jj)=ii+i*2*jj;
end
end

出了点bug～
赞一下 回复此楼
5楼2013-06-27 15:34:04
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

somomo91

专家顾问 (职业作家)

【答案】应助回帖

引用回帖:
5楼: Originally posted by 长路漫漫 at 2013-06-27 15:34:04
这里有些写错了，现在进行修正
做的矩阵式6×4的，
M=6;
N=4;
for ii=1:M
for jj=1:N
a(ii,jj)=ii+i*2*jj;
end
end

出了点bug～...

S 矩阵是一样的
Fortran下精度为 （E-016）
9.854419884063358E-016
6.434838164926571E-016
MATLAB 默认显示4位零，你可以在MATLAB下 用 format long 显示一下
赞一下 回复此楼
6楼2013-06-27 16:37:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

somomo91

专家顾问 (职业作家)
很明显可发现s基本一致，而u和s从第三列开始就区别非常大
你想说 u 和 v 区别大吧  ？
赞一下 回复此楼
7楼2013-06-27 16:38:36
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

长路漫漫

新虫 (初入文坛)
引用回帖:
7楼: Originally posted by somomo91 at 2013-06-27 16:38:36
很明显可发现s基本一致，而u和s从第三列开始就区别非常大
你想说 u 和 v 区别大吧  ？

嗯嗯～就是，因为，接下来的运算用到了U,V，所以需要较为精确点的
我网上找了很多代码，都无法匹配上matlab的，很抓狂，求指导！！
谢谢～
赞一下 回复此楼
8楼2013-06-27 16:59:43
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

somomo91

专家顾问 (职业作家)
引用回帖:
8楼: Originally posted by 长路漫漫 at 2013-06-27 16:59:43
嗯嗯～就是，因为，接下来的运算用到了U,V，所以需要较为精确点的
我网上找了很多代码，都无法匹配上matlab的，很抓狂，求指导！！
谢谢～...

那你希望是什么样的 ？
用FORTRAN的对，还是MATLAB的对 ？
赞一下 回复此楼
9楼2013-06-27 17:52:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

长路漫漫

新虫 (初入文坛)
引用回帖:
9楼: Originally posted by somomo91 at 2013-06-27 17:52:57
那你希望是什么样的 ？
用FORTRAN的对，还是MATLAB的对 ？...

啊？？有点懵。。。我希望？呵呵，没有希望，因为得到的SVD结果不和matlab近似，所以后来的运算，结果就很不对，所以，如果按这样说，我就是希望matlab对的
呵呵
赞一下 回复此楼
10楼2013-06-27 17:55:49
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 长路漫漫 的主题更新
161/2返回列表12下一页
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭