博文
fortran lapack的dgesvd和matlab SVD结果不一致，请大神指点！
||
对于A矩阵，matlab SVD分解结果如下：
A =
2.2700 -1.5400 1.1500 -1.9400
0.2800 -1.6700 0.9400 -0.7800
-0.4800 -3.0900 0.9900 -0.2100
1.0700 1.2200 0.7900 0.6300
-2.3500 2.9300 -1.4500 2.3000
0.6200 -7.3900 1.0300 -2.5700
>> [U,S,V] = svd(A,0)
U =
0.2774 -0.6003 -0.1277 0.1323
0.2020 -0.0301 0.2805 0.7034
0.2918 0.3348 0.6453 0.1906
-0.0938 -0.3699 0.6781 -0.5399
-0.4213 0.5266 0.0413 -0.0575
0.7816 0.3353 -0.1645 -0.3957
S =
9.9966 0 0 0
0 3.6831 0 0
0 0 1.3569 0
0 0 0 0.5000
V =
0.1921 -0.8030 0.0041 -0.5642
-0.8794 -0.3926 -0.0752 0.2587
0.2140 -0.2980 0.7827 0.5027
-0.3795 0.3351 0.6178 -0.6017
但在linux用fortran lapack得到的结果如下：
主函数：
* DGESVD Example Program Text
* NAG Copyright 2005.
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER (NIN=5,NOUT=6)
INTEGER MMAX, NB, NMAX, IERR
PARAMETER (MMAX=10,NB=64,NMAX=8)
INTEGER LDA, LDVT, LWORK
PARAMETER (LDA=MMAX,LDVT=NMAX,
+ LWORK=MMAX+4*NMAX+NB*(MMAX+NMAX))
* .. Local Scalars ..
DOUBLE PRECISION EPS, SERRBD
INTEGER I, IFAIL, INFO, J, LWKOPT, M, N
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), DUMMY(1,1), RCONDU(NMAX),
+ RCONDV(NMAX), S(NMAX), UERRBD(NMAX),
+ VERRBD(NMAX), VT(LDVT,NMAX), WORK(LWORK)
* .. External Functions ..
DOUBLE PRECISION DLAMCH
EXTERNAL DLAMCH
* .. External Subroutines ..
EXTERNAL DDISNA, DGESVD, X04CAF
* .. Executable Statements ..
OPEN(NIN,FILE='dgesvd-ex.d',STATUS='OLD',IOSTAT=IERR)
OPEN(NOUT,FILE='dgesvd-ex.r',STATUS='REPLACE',IOSTAT=IERR)
WRITE (NOUT,*) 'DGESVD Example Program Results'
WRITE (NOUT,*)
* Skip heading in data file
READ (NIN,*)
READ (NIN,*) M, N
IF (M.LE.MMAX .AND. N.LE.NMAX) THEN
*
* Read the m by n matrix A from data file
*
READ (NIN,*) ((A(I,J),J=1,N),I=1,M)
*
* Compute the singular values and left and right singular vectors
* of A (A = U*S*(V**T), m.ge.n)
*
CALL DGESVD('Overwrite A by U','Singular vectors (V)',M,N,A,
+ LDA,S,DUMMY,1,VT,LDVT,WORK,LWORK,INFO)
LWKOPT = WORK(1)
*
IF (INFO.EQ.0) THEN
*
* Print solution
*
WRITE (NOUT,*) 'Singular values'
WRITE (NOUT,99999) (S(J),J=1,N)
*
IFAIL = 0
CALL X04CAF('General',' ',M,N,A,LDA,
+ 'Left singular vectors (first n columns of U)',
+ IFAIL)
WRITE (NOUT,*)
CALL X04CAF('General',' ',N,N,VT,LDVT,
+ 'Right singular vectors by row (V**T)',IFAIL)
*
* Get the machine precision, EPS and compute the approximate
* error bound for the computed singular values. Note that for
* the 2-norm, S(1) = norm(A)
*
EPS = DLAMCH('Eps')
SERRBD = EPS*S(1)
*
* Call DDISNA (F08FLF) to estimate reciprocal condition
* numbers for the singular vectors
*
CALL DDISNA('Left',M,N,S,RCONDU,INFO)
CALL DDISNA('Right',M,N,S,RCONDV,INFO)
*
* Compute the error estimates for the singular vectors
*
DO 20 I = 1, N
UERRBD(I) = SERRBD/RCONDU(I)
VERRBD(I) = SERRBD/RCONDV(I)
20 CONTINUE
*
* Print the approximate error bounds for the singular values
* and vectors
*
WRITE (NOUT,*)
WRITE (NOUT,*) 'Error estimate for the singular values'
WRITE (NOUT,99998) SERRBD
WRITE (NOUT,*)
WRITE (NOUT,*)
+ 'Error estimates for the left singular vectors'
WRITE (NOUT,99998) (UERRBD(I),I=1,N)
WRITE (NOUT,*)
WRITE (NOUT,*)
+ 'Error estimates for the right singular vectors'
WRITE (NOUT,99998) (VERRBD(I),I=1,N)
ELSE
WRITE (NOUT,99997) 'Failure in DGESVD. INFO =', INFO
END IF
*
* Print workspace information
*
IF (LWORK.LT.LWKOPT) THEN
WRITE (NOUT,*)
WRITE (NOUT,99996) 'Optimum workspace required = ', LWKOPT,
+ 'Workspace provided = ', LWORK
END IF
ELSE
WRITE (NOUT,*) 'MMAX and/or NMAX too small'
END IF
STOP
*
99999 FORMAT (3X,(8F8.4))
99998 FORMAT (4X,1P,6E11.1)
99997 FORMAT (1X,A,I4)
99996 FORMAT (1X,A,I5,/1X,A,I5)
END
输入文件dgesvd-ex.d为：
DGESVD Example Program Data
6 4 :Values of M and N
2.27 -1.54 1.15 -1.94
0.28 -1.67 0.94 -0.78
-0.48 -3.09 0.99 -0.21
1.07 1.22 0.79 0.63
-2.35 2.93 -1.45 2.30
0.62 -7.39 1.03 -2.57 :End of matrix A
输出结果文件dgesvd-ex.r为：
DGESVD Example Program Results
Singular values
9.9966 3.6831 1.3569 0.5000
Left singular vectors (first n columns of U)
1 2 3 4
1 -0.2774 -0.6003 -0.1277 0.1323
2 -0.2020 -0.0301 0.2805 0.7034
3 -0.2918 0.3348 0.6453 0.1906
4 0.0938 -0.3699 0.6781 -0.5399
5 0.4213 0.5266 0.0413 -0.0575
6 -0.7816 0.3353 -0.1645 -0.3957
Right singular vectors by row (V**T)
1 2 3 4
1 -0.1921 0.8794 -0.2140 0.3795
2 -0.8030 -0.3926 -0.2980 0.3351
3 0.0041 -0.0752 0.7827 0.6178
4 -0.5642 0.2587 0.5027 -0.6017
Error estimate for the singular values
1.1E-15
Error estimates for the left singular vectors
1.8E-16 4.8E-16 1.3E-15 2.2E-15
Error estimates for the right singular vectors
1.8E-16 4.8E-16 1.3E-15 1.3E-15
这里输出的是V的转置，为什么U、V都是第一列符号相反呢？请大神指教！非常感谢
https://blog.sciencenet.cn/blog-350278-490543.html
上一篇：ubuntu GMT4.5.7配置文件.gmtdefaults4(本人用的)
下一篇：IGS台站在中国（GMT画图）
A =
2.2700 -1.5400 1.1500 -1.9400
0.2800 -1.6700 0.9400 -0.7800
-0.4800 -3.0900 0.9900 -0.2100
1.0700 1.2200 0.7900 0.6300
-2.3500 2.9300 -1.4500 2.3000
0.6200 -7.3900 1.0300 -2.5700
>> [U,S,V] = svd(A,0)
U =
0.2774 -0.6003 -0.1277 0.1323
0.2020 -0.0301 0.2805 0.7034
0.2918 0.3348 0.6453 0.1906
-0.0938 -0.3699 0.6781 -0.5399
-0.4213 0.5266 0.0413 -0.0575
0.7816 0.3353 -0.1645 -0.3957
S =
9.9966 0 0 0
0 3.6831 0 0
0 0 1.3569 0
0 0 0 0.5000
V =
0.1921 -0.8030 0.0041 -0.5642
-0.8794 -0.3926 -0.0752 0.2587
0.2140 -0.2980 0.7827 0.5027
-0.3795 0.3351 0.6178 -0.6017
但在linux用fortran lapack得到的结果如下：
主函数：
* DGESVD Example Program Text
* NAG Copyright 2005.
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER (NIN=5,NOUT=6)
INTEGER MMAX, NB, NMAX, IERR
PARAMETER (MMAX=10,NB=64,NMAX=8)
INTEGER LDA, LDVT, LWORK
PARAMETER (LDA=MMAX,LDVT=NMAX,
+ LWORK=MMAX+4*NMAX+NB*(MMAX+NMAX))
* .. Local Scalars ..
DOUBLE PRECISION EPS, SERRBD
INTEGER I, IFAIL, INFO, J, LWKOPT, M, N
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), DUMMY(1,1), RCONDU(NMAX),
+ RCONDV(NMAX), S(NMAX), UERRBD(NMAX),
+ VERRBD(NMAX), VT(LDVT,NMAX), WORK(LWORK)
* .. External Functions ..
DOUBLE PRECISION DLAMCH
EXTERNAL DLAMCH
* .. External Subroutines ..
EXTERNAL DDISNA, DGESVD, X04CAF
* .. Executable Statements ..
OPEN(NIN,FILE='dgesvd-ex.d',STATUS='OLD',IOSTAT=IERR)
OPEN(NOUT,FILE='dgesvd-ex.r',STATUS='REPLACE',IOSTAT=IERR)
WRITE (NOUT,*) 'DGESVD Example Program Results'
WRITE (NOUT,*)
* Skip heading in data file
READ (NIN,*)
READ (NIN,*) M, N
IF (M.LE.MMAX .AND. N.LE.NMAX) THEN
*
* Read the m by n matrix A from data file
*
READ (NIN,*) ((A(I,J),J=1,N),I=1,M)
*
* Compute the singular values and left and right singular vectors
* of A (A = U*S*(V**T), m.ge.n)
*
CALL DGESVD('Overwrite A by U','Singular vectors (V)',M,N,A,
+ LDA,S,DUMMY,1,VT,LDVT,WORK,LWORK,INFO)
LWKOPT = WORK(1)
*
IF (INFO.EQ.0) THEN
*
* Print solution
*
WRITE (NOUT,*) 'Singular values'
WRITE (NOUT,99999) (S(J),J=1,N)
*
IFAIL = 0
CALL X04CAF('General',' ',M,N,A,LDA,
+ 'Left singular vectors (first n columns of U)',
+ IFAIL)
WRITE (NOUT,*)
CALL X04CAF('General',' ',N,N,VT,LDVT,
+ 'Right singular vectors by row (V**T)',IFAIL)
*
* Get the machine precision, EPS and compute the approximate
* error bound for the computed singular values. Note that for
* the 2-norm, S(1) = norm(A)
*
EPS = DLAMCH('Eps')
SERRBD = EPS*S(1)
*
* Call DDISNA (F08FLF) to estimate reciprocal condition
* numbers for the singular vectors
*
CALL DDISNA('Left',M,N,S,RCONDU,INFO)
CALL DDISNA('Right',M,N,S,RCONDV,INFO)
*
* Compute the error estimates for the singular vectors
*
DO 20 I = 1, N
UERRBD(I) = SERRBD/RCONDU(I)
VERRBD(I) = SERRBD/RCONDV(I)
20 CONTINUE
*
* Print the approximate error bounds for the singular values
* and vectors
*
WRITE (NOUT,*)
WRITE (NOUT,*) 'Error estimate for the singular values'
WRITE (NOUT,99998) SERRBD
WRITE (NOUT,*)
WRITE (NOUT,*)
+ 'Error estimates for the left singular vectors'
WRITE (NOUT,99998) (UERRBD(I),I=1,N)
WRITE (NOUT,*)
WRITE (NOUT,*)
+ 'Error estimates for the right singular vectors'
WRITE (NOUT,99998) (VERRBD(I),I=1,N)
ELSE
WRITE (NOUT,99997) 'Failure in DGESVD. INFO =', INFO
END IF
*
* Print workspace information
*
IF (LWORK.LT.LWKOPT) THEN
WRITE (NOUT,*)
WRITE (NOUT,99996) 'Optimum workspace required = ', LWKOPT,
+ 'Workspace provided = ', LWORK
END IF
ELSE
WRITE (NOUT,*) 'MMAX and/or NMAX too small'
END IF
STOP
*
99999 FORMAT (3X,(8F8.4))
99998 FORMAT (4X,1P,6E11.1)
99997 FORMAT (1X,A,I4)
99996 FORMAT (1X,A,I5,/1X,A,I5)
END
输入文件dgesvd-ex.d为：
DGESVD Example Program Data
6 4 :Values of M and N
2.27 -1.54 1.15 -1.94
0.28 -1.67 0.94 -0.78
-0.48 -3.09 0.99 -0.21
1.07 1.22 0.79 0.63
-2.35 2.93 -1.45 2.30
0.62 -7.39 1.03 -2.57 :End of matrix A
输出结果文件dgesvd-ex.r为：
DGESVD Example Program Results
Singular values
9.9966 3.6831 1.3569 0.5000
Left singular vectors (first n columns of U)
1 2 3 4
1 -0.2774 -0.6003 -0.1277 0.1323
2 -0.2020 -0.0301 0.2805 0.7034
3 -0.2918 0.3348 0.6453 0.1906
4 0.0938 -0.3699 0.6781 -0.5399
5 0.4213 0.5266 0.0413 -0.0575
6 -0.7816 0.3353 -0.1645 -0.3957
Right singular vectors by row (V**T)
1 2 3 4
1 -0.1921 0.8794 -0.2140 0.3795
2 -0.8030 -0.3926 -0.2980 0.3351
3 0.0041 -0.0752 0.7827 0.6178
4 -0.5642 0.2587 0.5027 -0.6017
Error estimate for the singular values
1.1E-15
Error estimates for the left singular vectors
1.8E-16 4.8E-16 1.3E-15 2.2E-15
Error estimates for the right singular vectors
1.8E-16 4.8E-16 1.3E-15 1.3E-15
这里输出的是V的转置，为什么U、V都是第一列符号相反呢？请大神指教！非常感谢
https://blog.sciencenet.cn/blog-350278-490543.html
上一篇：ubuntu GMT4.5.7配置文件.gmtdefaults4(本人用的)
下一篇：IGS台站在中国（GMT画图）
扫一扫，分享此博文
全部作者的其他最新博文
全部精选博文导读
相关博文
- • Mutual Expressing Capabilities and Limits Among DIKWP （初学者版）
- • Expressing Capabilities Among Conscious Spaces （初学者版）
- • 大语言模型相亲择偶观：DIKWP实测（初学者版）
- • 12 philosophical problems by DIKWP Semantic Mathematic（初学者版）
- • Philosophical Technologiza to Artificial Consciousness（初学者版）
- • Discovering the Normal Distribution with New DIKWP （初学者版）