Документ взят из кэша поисковой машины. Адрес оригинального документа : http://foroff.phys.msu.ru/illposed/programs/Listing_01.txt Дата изменения: Mon Jul 7 00:04:43 2008 Дата индексирования: Mon Oct 1 22:10:45 2012 Кодировка: Поисковые слова: m 101

Listing file ptim.for


0001 SUBROUTINE PTIMR(AK,U0,A,B,C,D,N,M,Z,
*AN2,DL,H,C1,IMAX,ALFA,U,NU,IERR,
*AMU,C2,ALP,EPSF)
C * program to solve integral equations of the first kind
C * smoothing functional minimization is used
C * generalized residual principle used to get regularization parameter
C * transformation to the two-diagonal system used
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 EXTERNAL AK
0005 EXTERNAL PTICR0,PTIMRC,PTIMRS,PTIMRD,
*PTIMRP,PTIMRR,PTIMRQ,
*PTIMRN,PTIMRK,PTICR6,PTIMR0,PTIMRA,AUMINM
0006 DIMENSION U0(M),Z(N),U(NU)
0007 ICONT=0
C * ICONT- start/continue mode flag
C * ICONT=0-start mode
C * ICONT=1-continue mode
0008 110 CONTINUE
0009 DU=DSQRT(DL)
0010 DH=DSQRT(H)
C * working array mapping
C * name length what
C * AK(D) N*M system matrix or/and D=A*INV(S)
C * matrix of reflection vectors
C * C N stabilizer diagonal or/and
C * subdialonal of P'*P
C * B N over dialogal of the satbiliser
C * or/and over diagonal of P'*P
C * S1 N diagonal S, where C=S'*S
C * S2 N subdiagonal S
C * P1 N diagonal of the P-matrix
C * P2 N subdiagonal of the P-matrix
C * U1 N U1=R*D'*U0
C * U2 M U2=Q'*U0,
C * working array for PTIMR0
C * A N diagonal of P'*P
C * G M(N) working,
C * diagonal of P'*P+ALFA*E
C * AL N sweep coeffiecients
C * BT N sweep coeffiecients
C * working array starts
0011 NAK=1
0012 NC=NAK+N*M
0013 NB=NC+N
0014 NS1=NB+N
0015 NS2=NS1+N
0016 NP1=NS2+N
0017 NP2=NP1+N
0018 NU1=NP2+N
0019 NU2=NU1+N
0020 NA=NU2+M
0021 NG=NA+N
0022 NAL=NG+M
0023 NBT=NAL+N
0024 NOM=NBT+N
0025 IF(NOM-1.GT.NU) GOTO 64
C * K1,K2-iteration counters
0026 K1=0
0027 K2=0
0028 KJ=0
0029 AMU=0D0
0030 FUNC=0D0
C * RO - derivative wight in the norm
0031 RO=1.
0032 HX=(B-A)/(N-1.)
0033 HY=(D-C)/(M-1.)
0034 DD=HX/HY
C * EPS - accuracy of residual equation solution
0035 EPS=(C1-1.)*DL
0036 IF(ICONT.EQ.1) GOTO 111
C * getting matrix of the operator
0037 CALL PTICR0(AK,U(NAK),A,B,C,D,N,M)
C * getting matrix of stabilizer
0038 CALL PTIMRC(U(NC),U(NB),N,RO/HX**2)
C * transforming stabilizer matrix to the form S'*S
0039 CALL PTIMRS(U(NC),U(NB),U(NS1),U(NS2),N)
C * calculating D=A*INV(S)
0040 CALL PTIMRD(U(NS1),U(NS2),U(NAK),N,M)
C * applying QPR-algorithm to the matrix AK
0041 CALL PTIMR0(U(NAK),U(NP1),U(NP2),U(NG),N,M)
C * getting three-diaginal matrix P'*P
0042 CALL PTIMRP(U(NP1),U(NP2),U(NC),U(NB),U(NA),N,M)
C * getting U2=Q'*U0 to discrepancy calculations
0043 CALL PTIMRQ(U0,U(NAK),U(NU2),N,M,U(NG))
C * getting right-hand side of the system
C * U1=R*D'*U0=P'*U2
0044 CALL PTIMRR(U(NU1),U(NU2),U(NP1),U(NP2),N,M)
0045 111 CONTINUE
0046 NAN=NA+N-1
C * calculating measure of inconsistancy
0047 AN1=0.0
0048 IF(N.EQ.M)GOTO 13
0049 CALL PTICR6(U(NU2+N),U(NU2+N),M-N,AN1)
0050 AN1=AN1*HY
0051 13 CONTINUE
0052 IF(C2.LT.1.) GOTO 71
0053 CALL AUMINM(U(NA),ALP,IMAX,U(NC),
*U(NB),U(NG),U(NU1),U(NAL),
*U(NBT),Z,N,AN1,AN2,M,DU,DH,HX,HY,AN4,
*U(NP1),U(NP2),U(NU2),AMU,C2,KJ,EPSF)
0054 GOTO (70,71,69),KJ
0055 71 CONTINUE
C * starting regularization parameter loop
0056 DO 22 I=NA,NAN
0057 22 U(I+N)=U(I)+ALFA*DD
C * solving system by sweep method
0058 CALL PTIMRN(U(NC),U(NB),U(NG),U(NU1),U(NAL),U(NBT),Z,N)
C * calculating generalized discrepancy
0059 CALL PTIMRA(AN1,AN2,Z,N,M,DU,DH,HX,HY,
*AN4,U(NP1),U(NP2),U(NU2),FUNC,AMU,C2)
0060 IF(C1.LE.1)GOTO 100
0061 IF(ALFA.EQ.0.)GOTO 68
0062 IF(AN4.GE.EPS) GOTO 11
C * if discrepancy is less EPS then multiplay ALPHA by 2
C * while it would be greater EPS
0063 ALFA=2.*ALFA
0064 K1=K1+1
0065 IF(K1.EQ.IMAX) GOTO 65
0066 GOTO 13
C * setting start points for chord method
0067 11 F0=AN4
0068 X0=1./ALFA
0069 ALFA=ALFA*2.
0070 X=1./ALFA
0071 DO 5 I=NA,NAN
0072 5 U(I+N)=U(I)+ALFA*DD
C * sweep method
0073 CALL PTIMRN(U(NC),U(NB),U(NG),U(NU1),U(NAL),U(NBT),Z,N)
C * calculating discrepancy AN2
C * and generalized discrepancy AN4
0074 CALL PTIMRA(AN1,AN2,Z,N,M,DU,DH,HX,HY,
*AN4,U(NP1),U(NP2),U(NU2),FUNC,AMU,C2)
0075 14 CONTINUE
C * if discrepance is less EPS exit the program
0076 IF(DABS(AN4).LE.EPS) GOTO 100
0077 1 CONTINUE
C * if discrepance less -EPS - go to modified chord method
0078 IF(AN4.LT.-EPS) GOTO 2
0079 IF(ALFA.EQ.0.)GOTO 68
0080 K2=K2+1
0081 IF(K2.EQ.IMAX)GO TO 66
C * chord method formulae
0082 Y=X0-F0/(AN4-F0)*(X-X0)
0083 X0=X
0084 X=Y
0085 F0=AN4
0086 ALFA=1./X
0087 DO 7 I=NA,NAN
0088 7 U(I+N)=U(I)+ALFA*DD
C * sweep method
0089 CALL PTIMRN(U(NC),U(NB),U(NG),U(NU1),
*U(NAL),U(NBT),Z,N)
C * calculating discrepancy AN2
C * and generalized discrepance AN4
0090 CALL PTIMRA(AN1,AN2,Z,N,M,DU,DH,HX,HY,
*AN4,U(NP1),U(NP2),U(NU2),FUNC,AMU,C2)
0091 GOTO 14
0092 2 CONTINUE
C * if generalized discrepancy is negative
C * starting modified chord method
0093 F=AN4
0094 23 Y=X0+F*(X-X0)/(F-F0)
0095 ALFA=1./Y
0096 DO 15 I=NA,NAN
0097 15 U(I+N)=U(I)+ALFA*DD
C * sweep method
0098 CALL PTIMRN(U(NC),U(NB),U(NG),U(NU1),U(NAL),U(NBT),Z,N)
C * calculating discrepancy AN2
C * and generalized discrepance AN4
0099 CALL PTIMRA(AN1,AN2,Z,N,M,DU,DH,HX,HY,
*AN4,U(NP1),U(NP2),U(NU2),FUNC,AMU,C2)
C * if discrepancy >-EPS,HO 0100 IF(DABS(AN4).LE.EPS) GOTO 102
0101 K2=K2+1
0102 IF(K2.GE.IMAX) GOTO 67
0103 IF(ALFA.EQ.0.)GOTO 68
0104 IF(AN4.LE.-EPS) GOTO 37
0105 X0=Y
0106 F0=AN4
0107 GOTO 23
0108 37 CONTINUE
C * changing interval
0109 X=Y
0110 F=AN4
0111 GOTO 23
0112 ENTRY PTIMRE
C * entry to continue mode of calculations
0113 ICONT=1
0114 GOTO 110
0115 64 CONTINUE
C * working array is too short
0116 IERR=64
0117 GOTO 101
0118 65 CONTINUE
C * initial regularization parameter is too small
0119 IERR=65
0120 GOTO 101
0121 66 CONTINUE
C * IMAX iterations were made
0122 IERR=66
0123 GOTO 101
0124 67 CONTINUE
C * IMAX iterations by modified chord method were made
0125 IERR=67
0126 GOTO 101
0127 68 CONTINUE
C * ALFA=0 is set or found
0128 IERR=68
0129 GOTO 101
0130 69 CONTINUE
C * cannot localize minimum in IMAX iterations
0131 IERR=69
0132 GOTO 101
0133 70 CONTINUE
C * cannot find AMU with specified accuracy
0134 IERR=70
0135 GOTO 101
0136 102 CONTINUE
C * modified chord method was used
0137 IERR=1
0138 GOTO 101
0139 100 CONTINUE
C * solution was found by the base chord method
0140 IERR=0
0141 101 CONTINUE
C * performing inverse transformation to the old variables
0142 CALL PTIMRK(Z,U(NS1),U(NS2),U(NAK),N,M,
*Z,U(NG))
0143 RETURN
0144 END


0001 SUBROUTINE PTIMRC(C,B,N,RO)
C * gets stabilizer
C * C-diagonal, B-subdiagonal
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION C(N),B(N)
0005 DO 1 I=1,N
0006 C(I)=1.+2.*RO
0007 B(I)=-RO
0008 IF(I.EQ.1.OR.I.EQ.N) C(I)=1.+RO
0009 1 CONTINUE
0010 RETURN
0011 END


0001 SUBROUTINE PTIMRS(C,B,S1,S2,N)
C * transforms stabilizer to the form S'*S,
C * S1-diagonal, S2-subdiagonal
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION C(N),B(N),S1(N),S2(N)
0005 S1(1)=DSQRT(C(1))
0006 S2(1)=B(1)/S1(1)
0007 DO 1 I=2,N
0008 S1(I)=DSQRT(C(I)-(S2(I-1))**2)
0009 1 S2(I)=B(I)/S1(I)
0010 RETURN
0011 END


0001 SUBROUTINE PTIMRD(S1,S2,AA,N,M)
C * forms matrix D=A*INV(S)
C * S1 and S2 are diagonal and subdiagonal of the S
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION S1(N),S2(N),AA(M,N)
0005 DO 1 J=1,M
0006 1 AA(J,1)=AA(J,1)/S1(1)
0007 DO 2 I=2,N
0008 DO 3 J=1,M
0009 AA(J,I)=(AA(J,I)-S2(I-1)*AA(J,I-1))/S1(I)
0010 3 CONTINUE
0011 2 CONTINUE
0012 RETURN
0013 END


0001 SUBROUTINE PTIMR0(D,P1,P2,A,N,M)
C * QPR-algorithm of the matrix transformation
C * D(M,N)=Q*P*R
C * P1,P2 -diagonal and subdiagonal of the matrix P
C * matrix D(M,N) will contain:
C * in columns under diagonal -
C * left reflection vectors (Q)
C * in rows above diagonal -
C * right reflection vectors (R)
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 EXTERNAL PTIMR1
0005 DIMENSION D(M,N),P1(N),P2(N),A(M)
0006 DO 6 K=1,N
C * getting K-th column
0007 DO 11 J=K,M
0008 11 A(J)=D(J,K)
C * forming reflection vector
0009 CALL PTIMR1(A(K),M-K+1,B)
0010 A(K)=A(K)-B
0011 CALL PTIMR1(A(K),M-K+1,T)
0012 IF(T.EQ.0.)GOTO 12
0013 DO 13 J=K,M
0014 13 A(J)=A(J)/T
0015 GOTO 14
0016 12 CONTINUE
0017 DO 15 J=K,M
0018 15 A(J)=0.
0019 14 CONTINUE
C * A(K) - A(M) - reflection vector
C * writing reflection vector to the matrix D
0020 DO 16 J=K,M
0021 16 D(J,K)=A(J)
C * getting diagonal element - P1
0022 P1(K)=B
C * reflect columns from K+1 to N
0023 K1=K+1
0024 IF(K1.GT.N)GOTO 20
0025 DO 17 I=K1,N
0026 S=0.
0027 DO 18 J=K,M
0028 18 S=S+D(J,I)*A(J)
0029 DO 19 J=K,M
0030 19 D(J,I)=D(J,I)-2.*S*A(J)
0031 17 CONTINUE
0032 20 CONTINUE
0033 IF(K1.LT.N)GOTO 60
0034 IF(K1.EQ.N)P2(K)=D(K,K1)
0035 GOTO 6
0036 60 CONTINUE
C * getting K-th row of matrix D
0037 DO 21 I=K1,N
0038 21 A(I)=D(K,I)
C * form reflection vector
0039 CALL PTIMR1(A(K1),N-K1+1,B)
0040 A(K1)=A(K1)-B
0041 CALL PTIMR1(A(K1),N-K1+1,T)
0042 IF(T.EQ.0.)GOTO 22
0043 DO 23 I=K1,N
0044 23 A(I)=A(I)/T
0045 GOTO 24
0046 22 CONTINUE
0047 DO 25 I=K1,N
0048 25 A(I)=0.
0049 24 CONTINUE
C * B A(K1) - A(N) - reflection vector
C * writing reflection vector to K-th row of matrix
0050 DO 26 I=K1,N
0051 26 D(K,I)=A(I)
C * getting subdiagonal element P2
0052 P2(K)=B
C * reflecting rows from K+1 to M
0053 DO 27 J=K1,M
0054 S=0.
0055 DO 28 I=K1,N
0056 28 S=S+D(J,I)*A(I)
0057 DO 29 I=K1,N
0058 29 D(J,I)=D(J,I)-2.*S*A(I)
0059 27 CONTINUE
0060 6 CONTINUE
0061 RETURN
0062 END


0001 SUBROUTINE PTIMR1(A,N,S)
C * gets norm S of the vector A of length N
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION A(N)
0005 R=0.
0006 DO 1 I=1,N
0007 R1=DABS(A(I))
0008 IF(R1.GE.R)R=R1
0009 1 CONTINUE
0010 S=0.
0011 IF(R.EQ.0.)GOTO 999
0012 R1=0.
0013 DO 2 I=1,N
0014 R1=R1+(A(I)/R)**2
0015 2 CONTINUE
0016 S=DSQRT(R1)*R
0017 999 CONTINUE
0018 RETURN
0019 END


0001 SUBROUTINE PTIMRP(P1,P2,A,B,C,N,M)
C * gets matrix P'*P
C * C-diagonal, A-subdiagonal,
C * B-above diagonal of the result
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION P1(N),P2(N),A(N),B(N),C(N)
0005 C(1)=P1(1)**2
0006 B(1)=P1(1)*P2(1)
0007 A(2)=B(1)
0008 C(N)=P1(N)**2+P2(N-1)**2
0009 B(N)=0.
0010 A(1)=0.
0011 N1=N-1
0012 IF(2.GT.N1)GOTO 2
0013 DO 1 I=2,N1
0014 C(I)=P1(I)**2+P2(I-1)**2
0015 B(I)=P1(I)*P2(I)
0016 A(I+1)=B(I)
0017 1 CONTINUE
0018 2 CONTINUE
0019 RETURN
0020 END


0001 SUBROUTINE PTIMRQ(U0,D,U2,N,M,A)
C * gets vector U2=Q'*U0 for discrepancy calculation
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 EXTERNAL PTICR6,PTICR1
0005 DIMENSION U0(M),D(M,N),U2(M),A(M)
C * moving right-hand side U0
0006 DO 1 J=1,M
0007 1 U2(J)=U0(J)
C * reflecting vector U2
0008 DO 3 K=1,N
0009 DO 4 J=K,M
0010 4 A(J)=D(J,K)
C * in A - reflection vector
C * reflecting vector U2
0011 CALL PTICR6(A(K),U2(K),M-K+1,S)
0012 CALL PTICR1(U2(K),A(K),U2(K),M-K+1,-2.*S)
0013 3 CONTINUE
0014 RETURN
0015 END


0001 SUBROUTINE PTIMRR(U1,U2,P1,P2,N,M)
C * gets right-hand side of equation
C * U1=P'*U2 for sweep method
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION U2(M),P1(N),P2(N),U1(N)
0005 U1(1)=P1(1)*U2(1)
0006 DO 1 I=2,N
0007 U1(I)=U2(I-1)*P2(I-1)+U2(I)*P1(I)
0008 1 CONTINUE
0009 RETURN
0010 END


0001 SUBROUTINE PTIMRN(A,B,C,F,AL,BT,X,N)
C * finds solution of the system
C * A(I)*X(I-1)+C(I)*X(I)+B(I)*X(I+1)=F(I)
C * by sweep method
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION A(N),B(N),C(N),F(N),AL(N),BT(N)
0005 DIMENSION X(N)
0006 AL(1)=0.
0007 BT(1)=0.
0008 DO 1 I=2,N
0009 Y=-(C(I-1)+A(I-1)*AL(I-1))
0010 AL(I)=B(I-1)/Y
0011 BT(I)=(A(I-1)*BT(I-1)-F(I-1))/Y
0012 1 CONTINUE
C * sweep coeffiecients are found
0013 X(N)=(A(N)*BT(N)-F(N))/(-C(N)-A(N)*AL(N))
0014 DO 2 I=2,N
0015 J=N-I+2
0016 X(J-1)=AL(J)*X(J)+BT(J)
0017 2 CONTINUE
0018 RETURN
0019 END


0001 SUBROUTINE PTIMRA(AN1,AN2,Z,N,M,DL,
*H,HX,HY,AN4,P1,P2,U2,FUNC,AMU,C2)
C * gets discrepancy AN2
C * and generalized discrepancy AN4
C * AN1 - measure of inconsistency of the system
C * DL,H - errors (not a square ones)
C * in right side and operator
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 EXTERNAL PTICR6
0005 DIMENSION P1(N),P2(N),Z(N),U2(M)
C * calculating discrepancy
0006 SUM=0.
0007 DO 46 I=2,N
0008 46 SUM=SUM+(P1(I-1)*Z(I-1)+
+P2(I-1)*Z(I)-U2(I-1))**2
0009 SUM=SUM+(P1(N)*Z(N)-U2(N))**2
0010 IF (C2.GE.1.) GOTO 6
0011 AN2=SUM*HY+AN1
C * calculating generalized discrepancy
0012 AN4=AN2-DL**2-AN1
0013 IF(H.EQ.0.)GOTO 5
C * if h#0 then calculating norm of solution
0014 CALL PTICR6(Z,Z,N,S)
0015 S=DSQRT(S*HX)
C * S - W21 norm of solution
0016 AN4=AN2-(DL+H*S)**2-AN1
0017 5 CONTINUE
0018 GOTO 7
0019 6 CONTINUE
C * using inconsistence measure
0020 AN2=SUM*HY
0021 AN3=DSQRT(AN2)
0022 CALL PTICR6(Z,Z,N,S)
0023 S=DSQRT(S*HX)
0024 FUNC=AN3+H*S
0025 AN4=AN3-H*S-DL-AMU
0026 7 CONTINUE
0027 RETURN
0028 END


0001 SUBROUTINE PTIMRK(X,S1,S2,D,N,M,Z,A)
C * inverse transformation
C * Z=INV(S)*R'*X
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 EXTERNAL PTICR6,PTICR1
0005 DIMENSION X(N),S1(N),S2(N),D(M,N),
*Z(N),A(N)
0006 IF(3.GT.N)GOTO 3
0007 DO 1 II=3,N
0008 K=N-II+1
0009 K1=K+1
0010 DO 2 I=K1,N
0011 2 A(I)=D(K,I)
C * in A - reflection vector
C * reflecting vector X
0012 CALL PTICR6(A(K1),X(K1),N-K1+1,S)
0013 CALL PTICR1(X(K1),A(K1),X(K1),N-K1+1,-2.*S)
0014 1 CONTINUE
0015 3 CONTINUE
C * multiplaying by matrix inverse to S (solving the system)
0016 N1=N-1
0017 Z(N)=X(N)/S1(N)
0018 DO 5 I=1,N1
0019 J=N-I
0020 Z(J)=(X(J)-S2(J)*Z(J+1))/S1(J)
0021 5 CONTINUE
0022 RETURN
0023 END


0001 SUBROUTINE AUMINM(ZV,ALP,KMAX,A,B,C,F,
*AL,BT,X,N,AN1,AN2,M,DL,H,HX,HY,
*AN4,P1,P2,U2,AMU,C2,K2,EPSF)
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION ZV(N),A(N),B(N),C(N),F(N),AL(N)
0005 DIMENSION BT(N),X(N),P1(N),P2(N),U2(M)
0006 EXTERNAL PTIMRN,PTIMRA
0007 LI=0
0008 DD=HX/HY
0009 K3=0
0010 K2=0
0011 X1=0D0
0012 1 CONTINUE
0013 DO 13 I=1,N
0014 13 C(I)=ZV(I)+ALP*DD
0015 CALL PTIMRN(A,B,C,F,AL,BT,X,N)
0016 CALL PTIMRA(AN1,AN2,X,N,M,DL,H,HX,HY,
*AN4,P1,P2,U2,FUNC,AMU,C2)
0017 GOTO (2,5,6),LI
0018 F2=FUNC
0019 X2=ALP
0020 ALP=ALP*2D0
0021 LI=1
0022 GOTO 1
0023 2 F3=FUNC
0024 X3=ALP
0025 IF(F3-F2) 3,3,4
0026 3 F1=F2
0027 X1=X2
0028 F2=F3
0029 X2=X3
0030 ALP=ALP*2D0
0031 K3=K3+1
0032 IF(K3.GE.KMAX) GOTO 15
0033 GOTO 1
0034 4 X4=X3
0035 F4=F3
0036 AU=(DSQRT(5D0)-1D0)/2D0
0037 10 X3=X1+AU*(X4-X1)
0038 X2=X4-AU*(X4-X1)
0039 ALP=X2
0040 LI=2
0041 GOTO 1
0042 5 F2=FUNC
0043 ALP=X3
0044 LI=3
0045 GOTO 1
0046 6 F3=FUNC
0047 IF(DABS(X3-X2).LE.1D0/C2) GOTO 12
0048 IF(X1.EQ.0D0) GOTO 322
0049 IF(DABS(DMIN1(F2,F3)-DMAX1(F1,F4)).LE.EPSF) GOTO 12
0050 GOTO 321
0051 322 IF(F4-DMIN1(F2,F3).LE.EPSF) GOTO 12
0052 321 K2=K2+1
0053 IF(K2.GE.KMAX) GOTO 14
0054 IF(F3-F2) 7,8,9
0055 7 X1=X2
0056 F1=F2
0057 GOTO 10
0058 9 X4=X3
0059 F4=F3
0060 GOTO 10
0061 8 X4=(X4+X3)/2D0
0062 GOTO 10
0063 14 K2=1
0064 12 ALP=(X3+X2)/2D0
0065 AMU=(F3+F2)/2D0+DL
0066 K2=2
0067 GOTO 16
0068 15 K2=3
0069 16 RETURN
0070 END