REAL A(3,3), U(3,3), V(3,3), SIGMA(3), WORK(3),Y(3),C(3),Y0(3)
INTEGER I,IERR, J, M, N, NM
OPEN (6,FILE="SVD.OUT",STATUS="UNKNOWN",FORM="FORMATTED")
OPEN (5,FILE= "SVD.IN",STATUS="UNKNOWN",FORM="FORMATTED")
140 FORMAT(3I5)
150 FORMAT(4E15.7)
READ(5,140) NM,M,N
DO 131 I=1,M
READ(5,150) (A(I,J),J=1,N)
131 CONTINUE
READ (5,150) (Y(I),I=1,M)
CALL SVD(NM,M,N,A,SIGMA,.TRUE.,U,.TRUE.,V,IERR,WORK)
IF(IERR.NE.0) WRITE (6,2) IERR
2 FORMAT(15H TROUBLE.IERR=,I4)
WRITE(6,120)
120 FORMAT(/'МАТРИЦА А')
DO 121 I=1,M
WRITE(6,130) (A(I,J),J=1,N)
130 FORMAT(8E15.7)
121 CONTINUE
WRITE (6,160) (Y(I),I=1,N)
160 FORMAT(/'ПРАВЫЕ ЧАСТИ'/8E15.7)
210 FORMAT(/'СИНГУЛЯРНЫЕ ЧИСЛА')
WRITE(6,210)
DO 3 J=1,N
WRITE(6,6) SIGMA(J)
3 CONTINUE
SMA=SIGMA(1)
SMI=SIGMA(1)
DO 211 J=2,N
IF(SIGMA(J).GT.SMA) SMA=SIGMA(J)
IF(SIGMA(J).LT.SMI.AND.SIGMA(J).GT.0.) SMI=SIGMA(J)
211 CONTINUE
OBU=SMA/SMI
230 FORMAT(/'ЧИСЛО ОБУСЛОВЛЕННОСТИ=',E15.7)
WRITE(6,230) OBU
SIGMA1=0.
DO 30 J=1,N
IF(SIGMA(J).GT. SIGMA1) SIGMA1=SIGMA(J)
C(J)=0.
30 CONTINUE
TAU=SIGMA1*0.1E-6
DO 60 J=1,N
IF(SIGMA(J).LE.TAU) GO TO 60
S=0.
DO 40 I=1,N
S=S+U(I,J)*Y(I)
40 CONTINUE
S=S/SIGMA(J)
DO 50 I=1,N
C(I)=C(I) + S*V(I,J)
50 CONTINUE
60 CONTINUE
write (6,560)
WRITE (6,6) (C(I),I=1,3)
DO 322 J=1,N
SS=0.
DO 321 I=1,M
321 SS=A(J,I)*C(I)+SS
322 Y0(J)=SS
write (6,570)
WRITE (6,6) (Y0(I),I=1,3)
C WRITE(6,7)
C DO 4 I=1,M
C WRITE(6,6) (U(I,J),J=1,N)
C4 CONTINUE
C WRITE(6,7)
C DO 5 I=1,N
C WRITE(6,6) (V(I,J),J=1,N)
C5 CONTINUE
6 FORMAT(3E15.7)
560 format(2x,'roots')
570 format(2x,'right')
7 FORMAT(1H)
STOP
E N D
SUBROUTINE SVD(NM,M,N,A,W,MATU,U,MATV,V,IERR,RV1)
REAL A(NM,N),W(N),U(NM,N),V(NM,N),RV1(N)
LOGICAL MATU,MATV
IERR=0
DO 100 I=1,M
DO 100 J=1,N
U(I,J)=A(I,J)
100 CONTINUE
G=0.0
SCALE=0.0
ANORM=0.0
DO 300 I=1,N
L=I+1
RV1(I)=SCALE*G
G=0.0
S=0.0
SCALE=0.0
IF(I.GT.M) GO TO 210
DO 120 K=I,M
120 SCALE=SCALE+ABS(U(K,I))
IF(SCALE.EQ.0.0) GO TO 210
DO 130 K=I,M
U(K,I)=U(K,I)/SCALE
S=S+U(K,I)**2
130 CONTINUE
F=U(I,I)
G=-SIGN(SQRT(S),F)
H=F*G-S
U(I,I)=F-G
IF(I.EQ.N) GO TO 190
DO 150 J=L,N
S=0.0
DO 140 K=I,M
140 S=S+U(K,I)*U(K,J)
F=S/H
DO 150 K=I,M
U(K,J)=U(K,J)+F*U(K,I)
150 CONTINUE
190 DO 200 K=I,M
200 U(K,I)=SCALE*U(K,I)
210 W(I)=SCALE*G
G=0.0
S=0.0
SCALE=0.0
IF(I.GT.M.OR.I.EQ.N) GO TO 290
DO 220 K=L,N
220 SCALE=SCALE+ABS(U(I,K))
IF(SCALE.EQ.0.0) GO TO 290
DO 230 K=L,N
U(I,K)=U(I,K)/SCALE
S=S+U(I,K)**2
230 CONTINUE
F=U(I,L)
G=-SIGN(SQRT(S),F)
H=F*G-S
U(I,L)=F-G
DO 240 K=L,N
240 RV1(K)=U(I,K)/H
IF(I.EQ.M) GO TO 270
DO 260 J=L,M
S=0.0
DO 250 K=L,N
250 S=S+U(J,K)*U(I,K)
DO 260 K=L,N
U(J,K)=U(J,K)+S*RV1(K)
260 CONTINUE
270 DO 280 K=L,N
280 U(I,K)=SCALE*U(I,K)
290 ANORM=AMAX1(ANORM,ABS(W(I))+ABS(RV1(I)))
300 CONTINUE
IF(.NOT.MATV) GO TO 410
DO 400 II=1,N
I=N+1-II
IF(I.EQ.N) GO TO 390
IF(G.EQ.0.0) GO TO 360
DO 320 J=L,N
320 V(J,I)=(U(I,J)/U(I,L))/G
DO 350 J=L,N
S=0.0
DO 340 K=L,N
340 S=S+U(I,K)*V(K,J)
DO 350 K=L,N
V(K,J)=V(K,J)+S*V(K,I)
350 CONTINUE
360 DO 380 J=L,N
V(I,J)=0.0
V(J,I)=0.0
380 CONTINUE
390 V(I,I)=1.0
G=RV1(I)
L=I
400 CONTINUE
410 IF(.NOT.MATU) GO TO 510
MN=N
IF(M.LT.N) MN=M
DO 500 II=1,MN
I=MN+1-II
L=I+1
G=W(I)
IF(I.EQ.N) GO TO 430
DO 420 J=L,N
420 U(I,J)=0.0
430 IF(G.EQ.0.0) GO TO 475
IF(I.EQ.MN) GO TO 460
DO 450 J=L,N
S=0.0
DO 440 K=L,M
440 S=S+U(K,I)*U(K,J)
F=(S/U(I,I))/G
DO 450 K=I,M
U(K,J)=U(K,J)+F*U(K,I)
450 CONTINUE
460 DO 470 J=I,M
470 U(J,I)=U(J,I)/G
GO TO 490
475 DO 480 J=I,M
480 U(J,I)=0.0
490 U(I,I)=U(I,I)+1.0
500 CONTINUE
510 DO 700 KK=1,N
K1=N-KK
K=K1+1
ITS=0
520 DO 530 LL=1,K
L1=K-LL
L=L1+1
IF(ABS(RV1(L))+ANORM.EQ.ANORM) GO TO 565
IF(ABS(W(L1))+ANORM.EQ.ANORM) GO TO 540
530 CONTINUE
540 C=0.0
S=1.0
DO 560 I=L,K
F=S*RV1(I)
RV1(I)=C*RV1(I)
IF(ABS(F)+ANORM.EQ.ANORM) GO TO 565
G=W(I)
H=SQRT(F*F+G*G)
W(I)=H
C=G/H
S=-F/H
IF(.NOT.MATU) GO TO 560
DO 550 J=1,M
Y=U(J,L1)
Z=U(J,I)
U(J,L1)=Y*C+Z*S
U(J,I)=-Y*S+Z*C
550 CONTINUE
560 CONTINUE
565 Z=W(K)
IF(L.EQ.K) GO TO 650
IF(ITS.EQ.30) GO TO 1000
ITS=ITS+1
X=W(L)
Y=W(K1)
G=RV1(K1)
H=RV1(K)
F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y)
G=SQRT(F*F+1.0)
F=((X-Z)*(X+Z)+H*(Y/(F+SIGN(G,F))-H))/X
C=1.0
S=1.0
DO 600 I1=L,K1
I=I1+1
G=RV1(I)
Y=W(I)
H=S*G
G=C*G
Z=SQRT(F*F+H*H)
RV1(I1)=Z
C=F/Z
S=H/Z
F=X*C+G*S
G=-X*S+G*C
H=Y*S
Y=Y*C
IF(.NOT.MATV) GO TO 575
DO 570 J=1,N
X=V(J,I1)
Z=V(J,I)
V(J,I1)=X*C+Z*S
V(J,I)=-X*S+Z*C
570 CONTINUE
575 Z=SQRT(F*F+H*H)
W(I1)=Z
IF(Z.EQ.0.0) GO TO 580
C=F/Z
S=H/Z
580 F=C*G+S*Y
X=-S*G+C*Y
IF(.NOT.MATU) GO TO 600
DO 590 J=1,M
Y=U(J,I1)
Z=U(J,I)
U(J,I1)=Y*C+Z*S
U(J,I)=-Y*S+Z*C
590 CONTINUE
600 CONTINUE
RV1(L)=0.0
RV1(K)=F
W(K)=X
GO TO 520
650 IF(Z.GE.0.0) GO TO 700
W(K)=-Z
IF(.NOT.MATV) GO TO 700
DO 690 J=1,N
690 V(J,K)=-V(J,K)
700 CONTINUE
GO TO 1001
1000 IERR=K
1001 RETURN
E N D
ПРИЛОЖЕНИЕ 2. контрольный пример
Входные данные
(матрица изначально сингулярна – первая строка равна сумме второй и третьей с обратным знаком)
3 3 3
.3200000E 02.1400000E 02.7400000E 02
-0.2400000E 02 -0.1000000E 02 -0.5700000E 02
-0.8000000E 01 -0.4000000E 01 -0.1700000E 02
-0.1400000E 02 0.1300000E 02 0.1000000E 01
Полученный результат
МАТРИЦА А
.3200000E+02.1400000E+02.7400000E+02
-.2400000E+02 -.1000000E+02 -.5700000E+02
-.8000000E+01 -.4000000E+01 -.1700000E+02
ПРАВЫЕ ЧАСТИ
-.1400000E+02.1300000E+02.1000000E+01
СИНГУЛЯРНЫЕ ЧИСЛА
.1048255E+03
.7310871E-06
.1271749E+01
ЧИСЛО ОБУСЛОВЛЕННОСТИ=.1433830E+09
Корни
.1215394E+01.1821742E+01 -.1059419E+01
Правые корни после проверки
-.1400000E+02.1300000E+02.1000001E+01
Видно, что правые части соответствуют начальным данным. Решение верно.
[1] Матрица А эрмитоваесли она совпадает со своей комплексно сопряженной.
[2] Матрица А унитарная если 
, где – сопряженная матрица.
[3] Сингулярным разложением произвольной m ´ n –матрицы называется разложение вида, где U и V – ортогональные матрицы, а S – диагональная матрица с неотрицательными диагональными элементами. Диагональные элементы S (, i= 1,..., k, где k=min (m,n)) называются сингулярными числами А. Это множество чисел однозначно определяется матрицей А. Число ненулевых сингулярных чисел равно рангу А.
[4] Симметричная матрица положительно определена, если все ее собственные значения положительны. Положительно определенная матрица P обладает также тем свойством, что для всех.
[5] Симметричная матрица неотрицательно определена, если все ее собственные значения неотрицательны. Такая матрица P обладает также тем свойством, что для всех. Для произвольной m x n –матрицы А матрица симметрична и неотрицательно определена. Она положительно определена, если rankA=n.
[6] Обратной матрицей для квадратной невырожденной матрицы А называется такая матрица, для которой.
[7] Матрица перестановки - это квадратная матрица, столбцы которой получаются перестановкой столбцов единичной матрицы. Матрица перестановки ортогональна.
[8] Матрица А хессенбергова(верхняя хессенбергова) если для j<i– 1(сохраняется одна диагональ ниже главной диагонали). Если матрица симметричная то хессенбергова матрица становится трехдиагональной.
[9] Симметричная матрица А есть трехдиагональнаяпри для |i-j|> 1. Трехдиагональная матрица – это частный случай хесенберговой матрицы.