楼主
本人是一名大四学生，最近在做论文。老师要求我用FORTRAN做，可我对于FORTRAN语言一点都不懂，希望有贵人能指点我，并且给我一些意见，论文要求以及两个程序我都附在后面，希望达人不吝指教！或者将建议发到我的油箱cvbvc_1986@126.com，谢谢！
论文要求：
1、 研究一维地电模型正演程序代码；（1.得出结果需要哪些数据和文件？）
2、 研究反演算法（2.反演程序又和正演程序有什么关系？）
3、 实现反演算法并验证。（上面1.2两条是本人迷惘之处）
正演程序： dimension r(28),h(28),rr(50),ff(50),af(50)
complex z(28)
character*12 file1,file2,file3
pi=3.1415926
cc=0.0000001
c=8.0*pi**2*cc
write(*,*)'input m=? n=? '
read(*,*) m,n
do 5 i=1,n
write(*,*) 'number layar=',i
write(*,*) 'input depth resistance'
read(*,*) h(i),r(i)
5 continue
c write(*,*) 'please input freqence point file name'
c read(*,'(a)')
file1='fre.dat'
write(*,*)'please input model file name'
read(*,'(a)')file2
write(*,*)'please input fc result file name'
read(*,'(a)')file3
open(4,file=file2)
write(4,*)(h(i),r(i),i=1,n)
close(4)
open(8,file=file1,STATUS='OLD')
read(8,*)(ff(i),i=1,m)
close(8)
do 15 i=1,m
f=ff(i)
call fc(r,h,z,n,f)
re=real(z(1))
aim=aimag(z(1))
af(i)=atan(aim/re)*180/pi
rr(i)=(re**2+aim**2)/(c*f)
15 continue
open(3,file=file3)
write(3,*) (1/ff(k),af(k),k=1,m)
close(3)
stop
end
subroutine fc(r,h,z,n,f)
dimension h(28),r(28)
complex z(28),a,b,d
pi=3.1415926
cc=0.0000001
c=8.0*pi**2*cc
z(n)=sqrt(cmplx(0.0,c*r(n)*f))
do 10 j=1,n-1
j1=n-j
a=sqrt(cmplx(0.0,c*f*r(j1)))
b=sqrt(cmplx(0.0,c*f/r(j1)))
d=exp(-2.0*b*h(j1))
z(j1)=a*(a*(1-d)+z(j1+1)*(1+d))/(a*(1+d)+z(j1+1)*(1-d))
10 continue
return
end
反演程序： dimension r(15),h(15),rr(50),ac(50,29),q(29),u(50,29),v(29,29),
* e(29),ee(50),xz(29),ee1(50),r1(15),h1(15),ff(50),ft(29)
complex z(15),zz(15),zh(15),zr(15),z1r(15),z1h(15)
pi=3.1415926
cc=0.0000001
c=8.0*pi**2*cc
write(*,*)'************************************************'
write(*,*)'** **'
write(*,*)'** m=? n=? lmax=? **'
write(*,*)'** **'
write(*,*)'************************************************'
read(*,*) m,n,lmax
n1=2*n-1
open(8,file='fre.dat',status='old')
read(8,*) (ff(i),i=1,m)
do 2 i=1,m
2 ff(i)=1.0/ff(i)
write(*,*) 'm=',m, 'n=',n, 'lmax=',lmax
write(*,*)'do you have field data if you have input 1'
read(*,*) nn
if(nn.eq.1) goto 123
do 5 i=1,n
write(*,*) 'number layar=',i
write(*,*) 'input depth resistance'
read(*,*) h(i),r(i)
5 continue
open(6,file='rr.dat',status='new')
do 15 i=1,m
f=ff(i)
call fc(r,h,z,n,f)
re=real(z(1))
aim=aimag(z(1))
rr(i)=(re**2+aim**2)/(c*f)
15 continue
write(6,*) (rr(i),i=1,m)
close(6)
123 continue
c------ initial model---------------------------------------
open(5,file='model',status='old')
read(5,*) (h(i),r(i),i=1,n)
close(5)
c-----------------------------------------------------------
c****** input field data ******************
open(6,file='rr.dat',status='old')
read(6,*) (rr(i),i=1,m)
close(6)
c*******************************************
sum=0.0
do 600 i=1,m
f=ff(i)
call fc(r,h,z,n,f)
s=1.0/(f*c)*((real(z(1)))**2+(aimag(z(1)))**2)
sum=sum+(s-rr(i))**2
ee(i)=rr(i)-s
600 continue
q0=sum
l=0
1 arf=1.e-10
continue
l=l+1
write(*,*)'l=',l, 'q0=',q0
do 200 i=1,m
f=ff(i)
call fc(r,h,z,n,f)
call fc1(r,h,z,n,f,zz)
call fc2(r,h,z,zz,n,f,zh,zr)
call fc3(zz,zh,zr,n,f,z1h,z1r)
do 220 j=1,n
ac(i,j)=1.0/(c*f)*(z(1)*cmplx(real(z1r(j)),-aimag(z1r(j)))+
* z1r(j)*cmplx(real(z(1)),-aimag(z(1))))
220 continue
do 230 j=1,n-1
ac(i,j+n)=1.0/(c*f)*(z(1)*cmplx(real(z1h(j)),-aimag(z1h(j)))+
* z1h(j)*cmplx(real(z(1)),-aimag(z(1))))
230 continue
200 continue
call svd(m,n1,ac,q,u,v,e)
500 do 50 j=1,n1
yy=0.0
do 30 k=1,n1
if(q(k).le.td) goto 30
buk=0.0
do 40 i=1,m
buk=buk+ee(i)*u(i,k)
40 continue
qk=q(k)/(q(k)*q(k)+arf/(q(k)*q(k)))
yy=yy+qk*buk*v(j,k)
30 continue
xz(j)=yy
50 continue
do 650 i=1,n
650 r1(i)=r(i)+xz(i)
do 660 i=1,n-1
660 h1(i)=h(i)+xz(n+i)
sum1=0.0
do 680 i=1,m
f=ff(i)
call fc(r1,h1,z,n,f)
ss=1.0/(f*c)*((real(z(1)))**2+(aimag(z(1)))**2)
sum1=sum1+(rr(i)-ss)**2
ee1(i)=rr(i)-ss
680 continue
if(sum1.gt.q0) then
arf=arf*10.0
goto 500
else
do 690 j=1,n
r(j)=r1(j)
h(j)=h1(j)
690 continue
do 695 i=1,m
695 ee(i)=ee1(i)
q0=sum1
if(l.gt.lmax) goto 105
endif
c arf=arf/10.0
goto 1
105 write(*,*) 'iter gt lmax'
open(7,file='result',status='new')
write(7,*) (h(i),r(i),i=1,n)
close(7)
stop
end
subroutine fc(r,h,z,n,f)
dimension h(15),r(15)
complex z(15),a,b,d
pi=3.1415926
cc=0.0000001
c=8.0*pi**2*cc
z(n)=sqrt(cmplx(0.0,c*r(n)*f))
do 10 j=1,n-1
j1=n-j
a=sqrt(cmplx(0.0,c*f*r(j1)))
b=sqrt(cmplx(0.0,c*f/r(j1)))
d=exp(-2.0*b*h(j1))
z(j1)=a*(a*(1-d)+z(j1+1)*(1+d))/(a*(1+d)+z(j1+1)*(1-d))
10 continue
return
end
subroutine fc1(r,h,z,n,f,zz)
complex zz(15),z(15),a,b,d,e
dimension r(15),h(15)
pi=3.1415926
cc=0.0000001
c=8.0*pi**2*cc
do 10 i=1,n-1
a=sqrt(cmplx(0.0,c*f*r(i)))
b=sqrt(cmplx(0.0,c*f/r(i)))
d=exp(-2.0*b*h(i))
e=a*(1+d)+z(i+1)*(1-d)
10 zz(i)=4.0*a**2/e**2*d
return
end
subroutine fc2(r,h,z,zz,n,f,zh,zr)
complex zz(15),z(15),zh(15),zr(15),a,b
dimension r(15),h(15)
pi=3.1415926
cc=0.0000001
c=8.0*pi**2*cc
a=cmplx(0.0,f*c)
b=sqrt(cmplx(0.0,c*f*r(n)))
do 10 i=1,n-1
10 zh(i)=(a-z(i+1)**2/r(i))*zz(i)
do 20 i=1,n-1
zr(i)=1.0/(2.0*r(i))*(z(i)-z(i+1)*zz(i)-h(i)*zh(i))
20 continue
zr(n)=1.0/2.0*a/b
return
end
subroutine fc3(zz,zh,zr,n,f,z1h,z1r)
complex zz(15),zh(15),zr(15),z1h(15),z1r(15),x
do 10 i=1,n-1
x=cmplx(1.0,0.0)
if(i.eq.1) then
z1h(i)=zh(1)
z1r(i)=zr(1)
else
do 20 j=1,i-1
20 x=x*zz(j)
z1h(i)=x*zh(i)
z1r(i)=x*zr(i)
endif
10 continue
z1r(n)=x*zz(n-1)*zr(n)
return
end
C THIS IS SVD PROGRAM
C SUBROUTINE SVD(M,N,A,Q,U,V,E)
DIMENSION A(50,29),Q(29),U(50,29),V(29,29),E(29)
LOGICAL WITHU,WITHV
EPS=1.E-14
TOL=1.E-16
WITHU=.TRUE.
WITHV=.TRUE.
DO 10 I=1,M
DO 10 J=1,N
10 U(I,J)=A(I,J)
G=0.
X=0.
C HOUSHOLDER'S REDUCTION TO BIDIAGNAL FORM
DO 200 I=1,N
E(I)=G
S=0.
L=I+1
DO 20 J=I,M
20 S=S+U(J,I)**2
IF(S.GE.TOL)GOTO 30
G=0.
GOTO 90
30 F=U(I,I)
IF(F.LT.0.)GOTO 40
G=-SQRT(S)
GOTO 50
40 G=SQRT(S)
50 H=F*G-S
U(I,I)=F-G
IF(L.GT.N)GOTO 90
DO 80 J=L,N
S=0.
DO 60 K=I,M
60 S=S+U(K,I)*U(K,J)
F=S/H
DO 70 K=I,M
70 U(K,J)=U(K,J)+F*U(K,I)
80 CONTINUE
90 Q(I)=G
S=0.
IF(L.GT.N)GOTO 102
DO 100 J=L,N
100 S=S+U(I,J)**2
102 IF(S.GE.TOL)GOTO 110
G=0.0
GOTO 180
110 F=U(I,I+1)
IF(F.LT.0.)GOTO 120
G=-SQRT(S)
GOTO 130
120 G=SQRT(S)
130 H=F*G-S
U(I,I+1)=F-G
DO 140 J=L,N
140 E(J)=U(I,J)/H
DO 170 J=L,M
S=0.
DO 150 K=L,N
150 S=S+U(J,K)*U(I,K)
DO 160 K=L,N
160 U(J,K)=U(J,K)+S*E(K)
170 CONTINUE
180 Y=ABS(Q(I))+ABS(E(I))
IF(Y.LE.X)GOTO 200
X=Y
200 CONTINUE
C ACCUMULATION OF RIGHT-HAND TRANSFORMATION
IF(WITHV)GO TO 210
GOTO 290
210 DO 280 II=1,N
I=N+1-II
IF(G.EQ.0.)GOTO 260
H=U(I,I+1)*G
DO 220 J=L,N
220 V(J,I)=U(I,J)/H
DO 250 J=L,N
S=0.
DO 230 K=L,N
230 S=S+U(I,K)*V(K,J)
DO 240 K=L,N
240 V(K,J)=V(K,J)+S*V(K,I)
250 CONTINUE
260 IF(L.GT.N)GOTO 272
DO 270 J=L,N
V(I,J)=0.
270 V(J,I)=0.
272 V(I,I)=1.
G=E(I)
L=I
280 CONTINUE
C ACCUMULATION OF LEFT-HANDM TRANSFORMATION
290 IF(WITHU)GOTO 300
GOTO 400
300 DO 390 II=1,N
I=N+1-II
L=I+1
G=Q(I)
IF(L.GT.N)GOTO 311
DO 310 J=L,N
310 U(I,J)=0.
311 IF(G.EQ.0.)GOTO 360
H=U(I,I)*G
IF(L.GT.N)GOTO 341
DO 340 J=L,N
S=0.
DO 320 K=L,M
320 S=S+U(K,I)*U(K,J)
F=S/H
DO 330 K=I,M
330 U(K,J)=U(K,J)+F*U(K,I)
340 CONTINUE
341 DO 350 J=I,M
350 U(J,I)=U(J,I)/G
GOTO 380
360 DO 370 J=I,M
370 U(J,I)=0.
380 U(I,I)=U(I,I)+1.
390 CONTINUE
C DIAGONALIZATION OF THE BIDIAGONAL FORM
400 EPS=EPS*X
DO 580 KK=1,N
K=N+1-KK
410 DO 420 LL=1,K
L=K+1-LL
IF(ABS(E(L)).LE.EPS)GOTO 470
IF(ABS(Q(L-1)).LE.EPS)GOTO 430
420 CONTINUE
C CANCELLATION OF E(L) IF L>1
430 C=0.
S=1.
L1=L-1
DO 460 I=L,K
F=E(I)*S
E(I)=C*E(I)
IF(ABS(F).LE.EPS)GOTO 470
G=Q(I)
H=SQRT(F*F+G*G)
Q(I)=H
C=G/H
S=-F/H
IF(WITHU) GOTO 440
GOTO 460
440 DO 450 J=1,M
Y=U(J,L1)
Z=U(J,I)
U(J,L1)=Y*C+Z*S
450 U(J,I)=-Y*S+Z*C
460 CONTINUE
470 Z=Q(K)
IF(L.EQ.K)GOTO 550
C SHIFT FORM BOTTOM IXI MINOR
X=Q(L)
Y=Q(K-1)
G=E(K-1)
H=E(K)
F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.*H*Y)
G=SQRT(F*F+1.)
IF(F.LT.0.)GOTO 480
F=((X-Z)*(X+Z)+H*(Y/(F+G)-H))/X
GOTO 490
480 F=((X-Z)*(X+Z)+H*(Y/(F-G)-H))/X
C NEXT QR TRANSFORMATION
490 S=1.
C=1.
LO=L+1
DO 540 I=LO,K
G=E(I)
Y=Q(I)
H=S*G
G=C*G
E(I-1)=SQRT(F*F+H*H)
Z=E(I-1)
C=F/Z
S=H/Z
F=X*C+G*S
G=-X*S+G*C
H=Y*S
Y=Y*C
IF(WITHV)GOTO 495
GOTO 510
495 DO 500 J=1,N
X=V(J,I-1)
Z=V(J,I)
V(J,I-1)=X*C+Z*S
500 V(J,I)=-X*S+Z*C
510 Q(I-1)=SQRT(F*F+H*H)
Z=Q(I-1)
C=F/Z
S=H/Z
F=C*G+S*Y
X=-S*G+C*Y
IF(WITHU)GOTO 520
GOTO 540
520 DO 530 J=1,M
Y=U(J,I-1)
Z=U(J,I)
U(J,I-1)=Y*C+Z*S
530 U(J,I)=-Y*S+Z*C
540 CONTINUE
E(L)=0.
E(K)=F
Q(K)=X
GOTO 410
550 IF(Z.GE.0.)GOTO 580
C Q(K) IS MODE NON-NEGATIVE
Q(K)=-Z
IF(WITHV) GOTO 560
GO TO 580
560 DO 570 J=1,N
570 V(J,K)=-V(J,K)
580 CONTINUE
c WRITE(*,581)
581 FORMAT(1X,' SINGULE----VALUE')
c WRITE(*,590)(Q(II),II=1,N)
590 FORMAT(1X,8(2X,E13.7))
RETURN
END
xiexie !
论文要求：
1、 研究一维地电模型正演程序代码；（1.得出结果需要哪些数据和文件？）
2、 研究反演算法（2.反演程序又和正演程序有什么关系？）
3、 实现反演算法并验证。（上面1.2两条是本人迷惘之处）
正演程序： dimension r(28),h(28),rr(50),ff(50),af(50)
complex z(28)
character*12 file1,file2,file3
pi=3.1415926
cc=0.0000001
c=8.0*pi**2*cc
write(*,*)'input m=? n=? '
read(*,*) m,n
do 5 i=1,n
write(*,*) 'number layar=',i
write(*,*) 'input depth resistance'
read(*,*) h(i),r(i)
5 continue
c write(*,*) 'please input freqence point file name'
c read(*,'(a)')
file1='fre.dat'
write(*,*)'please input model file name'
read(*,'(a)')file2
write(*,*)'please input fc result file name'
read(*,'(a)')file3
open(4,file=file2)
write(4,*)(h(i),r(i),i=1,n)
close(4)
open(8,file=file1,STATUS='OLD')
read(8,*)(ff(i),i=1,m)
close(8)
do 15 i=1,m
f=ff(i)
call fc(r,h,z,n,f)
re=real(z(1))
aim=aimag(z(1))
af(i)=atan(aim/re)*180/pi
rr(i)=(re**2+aim**2)/(c*f)
15 continue
open(3,file=file3)
write(3,*) (1/ff(k),af(k),k=1,m)
close(3)
stop
end
subroutine fc(r,h,z,n,f)
dimension h(28),r(28)
complex z(28),a,b,d
pi=3.1415926
cc=0.0000001
c=8.0*pi**2*cc
z(n)=sqrt(cmplx(0.0,c*r(n)*f))
do 10 j=1,n-1
j1=n-j
a=sqrt(cmplx(0.0,c*f*r(j1)))
b=sqrt(cmplx(0.0,c*f/r(j1)))
d=exp(-2.0*b*h(j1))
z(j1)=a*(a*(1-d)+z(j1+1)*(1+d))/(a*(1+d)+z(j1+1)*(1-d))
10 continue
return
end
反演程序： dimension r(15),h(15),rr(50),ac(50,29),q(29),u(50,29),v(29,29),
* e(29),ee(50),xz(29),ee1(50),r1(15),h1(15),ff(50),ft(29)
complex z(15),zz(15),zh(15),zr(15),z1r(15),z1h(15)
pi=3.1415926
cc=0.0000001
c=8.0*pi**2*cc
write(*,*)'************************************************'
write(*,*)'** **'
write(*,*)'** m=? n=? lmax=? **'
write(*,*)'** **'
write(*,*)'************************************************'
read(*,*) m,n,lmax
n1=2*n-1
open(8,file='fre.dat',status='old')
read(8,*) (ff(i),i=1,m)
do 2 i=1,m
2 ff(i)=1.0/ff(i)
write(*,*) 'm=',m, 'n=',n, 'lmax=',lmax
write(*,*)'do you have field data if you have input 1'
read(*,*) nn
if(nn.eq.1) goto 123
do 5 i=1,n
write(*,*) 'number layar=',i
write(*,*) 'input depth resistance'
read(*,*) h(i),r(i)
5 continue
open(6,file='rr.dat',status='new')
do 15 i=1,m
f=ff(i)
call fc(r,h,z,n,f)
re=real(z(1))
aim=aimag(z(1))
rr(i)=(re**2+aim**2)/(c*f)
15 continue
write(6,*) (rr(i),i=1,m)
close(6)
123 continue
c------ initial model---------------------------------------
open(5,file='model',status='old')
read(5,*) (h(i),r(i),i=1,n)
close(5)
c-----------------------------------------------------------
c****** input field data ******************
open(6,file='rr.dat',status='old')
read(6,*) (rr(i),i=1,m)
close(6)
c*******************************************
sum=0.0
do 600 i=1,m
f=ff(i)
call fc(r,h,z,n,f)
s=1.0/(f*c)*((real(z(1)))**2+(aimag(z(1)))**2)
sum=sum+(s-rr(i))**2
ee(i)=rr(i)-s
600 continue
q0=sum
l=0
1 arf=1.e-10
continue
l=l+1
write(*,*)'l=',l, 'q0=',q0
do 200 i=1,m
f=ff(i)
call fc(r,h,z,n,f)
call fc1(r,h,z,n,f,zz)
call fc2(r,h,z,zz,n,f,zh,zr)
call fc3(zz,zh,zr,n,f,z1h,z1r)
do 220 j=1,n
ac(i,j)=1.0/(c*f)*(z(1)*cmplx(real(z1r(j)),-aimag(z1r(j)))+
* z1r(j)*cmplx(real(z(1)),-aimag(z(1))))
220 continue
do 230 j=1,n-1
ac(i,j+n)=1.0/(c*f)*(z(1)*cmplx(real(z1h(j)),-aimag(z1h(j)))+
* z1h(j)*cmplx(real(z(1)),-aimag(z(1))))
230 continue
200 continue
call svd(m,n1,ac,q,u,v,e)
500 do 50 j=1,n1
yy=0.0
do 30 k=1,n1
if(q(k).le.td) goto 30
buk=0.0
do 40 i=1,m
buk=buk+ee(i)*u(i,k)
40 continue
qk=q(k)/(q(k)*q(k)+arf/(q(k)*q(k)))
yy=yy+qk*buk*v(j,k)
30 continue
xz(j)=yy
50 continue
do 650 i=1,n
650 r1(i)=r(i)+xz(i)
do 660 i=1,n-1
660 h1(i)=h(i)+xz(n+i)
sum1=0.0
do 680 i=1,m
f=ff(i)
call fc(r1,h1,z,n,f)
ss=1.0/(f*c)*((real(z(1)))**2+(aimag(z(1)))**2)
sum1=sum1+(rr(i)-ss)**2
ee1(i)=rr(i)-ss
680 continue
if(sum1.gt.q0) then
arf=arf*10.0
goto 500
else
do 690 j=1,n
r(j)=r1(j)
h(j)=h1(j)
690 continue
do 695 i=1,m
695 ee(i)=ee1(i)
q0=sum1
if(l.gt.lmax) goto 105
endif
c arf=arf/10.0
goto 1
105 write(*,*) 'iter gt lmax'
open(7,file='result',status='new')
write(7,*) (h(i),r(i),i=1,n)
close(7)
stop
end
subroutine fc(r,h,z,n,f)
dimension h(15),r(15)
complex z(15),a,b,d
pi=3.1415926
cc=0.0000001
c=8.0*pi**2*cc
z(n)=sqrt(cmplx(0.0,c*r(n)*f))
do 10 j=1,n-1
j1=n-j
a=sqrt(cmplx(0.0,c*f*r(j1)))
b=sqrt(cmplx(0.0,c*f/r(j1)))
d=exp(-2.0*b*h(j1))
z(j1)=a*(a*(1-d)+z(j1+1)*(1+d))/(a*(1+d)+z(j1+1)*(1-d))
10 continue
return
end
subroutine fc1(r,h,z,n,f,zz)
complex zz(15),z(15),a,b,d,e
dimension r(15),h(15)
pi=3.1415926
cc=0.0000001
c=8.0*pi**2*cc
do 10 i=1,n-1
a=sqrt(cmplx(0.0,c*f*r(i)))
b=sqrt(cmplx(0.0,c*f/r(i)))
d=exp(-2.0*b*h(i))
e=a*(1+d)+z(i+1)*(1-d)
10 zz(i)=4.0*a**2/e**2*d
return
end
subroutine fc2(r,h,z,zz,n,f,zh,zr)
complex zz(15),z(15),zh(15),zr(15),a,b
dimension r(15),h(15)
pi=3.1415926
cc=0.0000001
c=8.0*pi**2*cc
a=cmplx(0.0,f*c)
b=sqrt(cmplx(0.0,c*f*r(n)))
do 10 i=1,n-1
10 zh(i)=(a-z(i+1)**2/r(i))*zz(i)
do 20 i=1,n-1
zr(i)=1.0/(2.0*r(i))*(z(i)-z(i+1)*zz(i)-h(i)*zh(i))
20 continue
zr(n)=1.0/2.0*a/b
return
end
subroutine fc3(zz,zh,zr,n,f,z1h,z1r)
complex zz(15),zh(15),zr(15),z1h(15),z1r(15),x
do 10 i=1,n-1
x=cmplx(1.0,0.0)
if(i.eq.1) then
z1h(i)=zh(1)
z1r(i)=zr(1)
else
do 20 j=1,i-1
20 x=x*zz(j)
z1h(i)=x*zh(i)
z1r(i)=x*zr(i)
endif
10 continue
z1r(n)=x*zz(n-1)*zr(n)
return
end
C THIS IS SVD PROGRAM
C SUBROUTINE SVD(M,N,A,Q,U,V,E)
DIMENSION A(50,29),Q(29),U(50,29),V(29,29),E(29)
LOGICAL WITHU,WITHV
EPS=1.E-14
TOL=1.E-16
WITHU=.TRUE.
WITHV=.TRUE.
DO 10 I=1,M
DO 10 J=1,N
10 U(I,J)=A(I,J)
G=0.
X=0.
C HOUSHOLDER'S REDUCTION TO BIDIAGNAL FORM
DO 200 I=1,N
E(I)=G
S=0.
L=I+1
DO 20 J=I,M
20 S=S+U(J,I)**2
IF(S.GE.TOL)GOTO 30
G=0.
GOTO 90
30 F=U(I,I)
IF(F.LT.0.)GOTO 40
G=-SQRT(S)
GOTO 50
40 G=SQRT(S)
50 H=F*G-S
U(I,I)=F-G
IF(L.GT.N)GOTO 90
DO 80 J=L,N
S=0.
DO 60 K=I,M
60 S=S+U(K,I)*U(K,J)
F=S/H
DO 70 K=I,M
70 U(K,J)=U(K,J)+F*U(K,I)
80 CONTINUE
90 Q(I)=G
S=0.
IF(L.GT.N)GOTO 102
DO 100 J=L,N
100 S=S+U(I,J)**2
102 IF(S.GE.TOL)GOTO 110
G=0.0
GOTO 180
110 F=U(I,I+1)
IF(F.LT.0.)GOTO 120
G=-SQRT(S)
GOTO 130
120 G=SQRT(S)
130 H=F*G-S
U(I,I+1)=F-G
DO 140 J=L,N
140 E(J)=U(I,J)/H
DO 170 J=L,M
S=0.
DO 150 K=L,N
150 S=S+U(J,K)*U(I,K)
DO 160 K=L,N
160 U(J,K)=U(J,K)+S*E(K)
170 CONTINUE
180 Y=ABS(Q(I))+ABS(E(I))
IF(Y.LE.X)GOTO 200
X=Y
200 CONTINUE
C ACCUMULATION OF RIGHT-HAND TRANSFORMATION
IF(WITHV)GO TO 210
GOTO 290
210 DO 280 II=1,N
I=N+1-II
IF(G.EQ.0.)GOTO 260
H=U(I,I+1)*G
DO 220 J=L,N
220 V(J,I)=U(I,J)/H
DO 250 J=L,N
S=0.
DO 230 K=L,N
230 S=S+U(I,K)*V(K,J)
DO 240 K=L,N
240 V(K,J)=V(K,J)+S*V(K,I)
250 CONTINUE
260 IF(L.GT.N)GOTO 272
DO 270 J=L,N
V(I,J)=0.
270 V(J,I)=0.
272 V(I,I)=1.
G=E(I)
L=I
280 CONTINUE
C ACCUMULATION OF LEFT-HANDM TRANSFORMATION
290 IF(WITHU)GOTO 300
GOTO 400
300 DO 390 II=1,N
I=N+1-II
L=I+1
G=Q(I)
IF(L.GT.N)GOTO 311
DO 310 J=L,N
310 U(I,J)=0.
311 IF(G.EQ.0.)GOTO 360
H=U(I,I)*G
IF(L.GT.N)GOTO 341
DO 340 J=L,N
S=0.
DO 320 K=L,M
320 S=S+U(K,I)*U(K,J)
F=S/H
DO 330 K=I,M
330 U(K,J)=U(K,J)+F*U(K,I)
340 CONTINUE
341 DO 350 J=I,M
350 U(J,I)=U(J,I)/G
GOTO 380
360 DO 370 J=I,M
370 U(J,I)=0.
380 U(I,I)=U(I,I)+1.
390 CONTINUE
C DIAGONALIZATION OF THE BIDIAGONAL FORM
400 EPS=EPS*X
DO 580 KK=1,N
K=N+1-KK
410 DO 420 LL=1,K
L=K+1-LL
IF(ABS(E(L)).LE.EPS)GOTO 470
IF(ABS(Q(L-1)).LE.EPS)GOTO 430
420 CONTINUE
C CANCELLATION OF E(L) IF L>1
430 C=0.
S=1.
L1=L-1
DO 460 I=L,K
F=E(I)*S
E(I)=C*E(I)
IF(ABS(F).LE.EPS)GOTO 470
G=Q(I)
H=SQRT(F*F+G*G)
Q(I)=H
C=G/H
S=-F/H
IF(WITHU) GOTO 440
GOTO 460
440 DO 450 J=1,M
Y=U(J,L1)
Z=U(J,I)
U(J,L1)=Y*C+Z*S
450 U(J,I)=-Y*S+Z*C
460 CONTINUE
470 Z=Q(K)
IF(L.EQ.K)GOTO 550
C SHIFT FORM BOTTOM IXI MINOR
X=Q(L)
Y=Q(K-1)
G=E(K-1)
H=E(K)
F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.*H*Y)
G=SQRT(F*F+1.)
IF(F.LT.0.)GOTO 480
F=((X-Z)*(X+Z)+H*(Y/(F+G)-H))/X
GOTO 490
480 F=((X-Z)*(X+Z)+H*(Y/(F-G)-H))/X
C NEXT QR TRANSFORMATION
490 S=1.
C=1.
LO=L+1
DO 540 I=LO,K
G=E(I)
Y=Q(I)
H=S*G
G=C*G
E(I-1)=SQRT(F*F+H*H)
Z=E(I-1)
C=F/Z
S=H/Z
F=X*C+G*S
G=-X*S+G*C
H=Y*S
Y=Y*C
IF(WITHV)GOTO 495
GOTO 510
495 DO 500 J=1,N
X=V(J,I-1)
Z=V(J,I)
V(J,I-1)=X*C+Z*S
500 V(J,I)=-X*S+Z*C
510 Q(I-1)=SQRT(F*F+H*H)
Z=Q(I-1)
C=F/Z
S=H/Z
F=C*G+S*Y
X=-S*G+C*Y
IF(WITHU)GOTO 520
GOTO 540
520 DO 530 J=1,M
Y=U(J,I-1)
Z=U(J,I)
U(J,I-1)=Y*C+Z*S
530 U(J,I)=-Y*S+Z*C
540 CONTINUE
E(L)=0.
E(K)=F
Q(K)=X
GOTO 410
550 IF(Z.GE.0.)GOTO 580
C Q(K) IS MODE NON-NEGATIVE
Q(K)=-Z
IF(WITHV) GOTO 560
GO TO 580
560 DO 570 J=1,N
570 V(J,K)=-V(J,K)
580 CONTINUE
c WRITE(*,581)
581 FORMAT(1X,' SINGULE----VALUE')
c WRITE(*,590)(Q(II),II=1,N)
590 FORMAT(1X,8(2X,E13.7))
RETURN
END
xiexie !
