Документ взят из кэша поисковой машины. Адрес оригинального документа : http://foroff.phys.msu.ru/illposed/programs/all_listings.txt Дата изменения: Mon Jul 7 00:04:43 2008 Дата индексирования: Mon Oct 1 22:15:00 2012 Кодировка: IBM-866 Поисковые слова: 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

0001 SUBROUTINE PTIZR(AK,U0,A,B,C,D,N,M,Z,
*AN2,DL,H,C1,IMAX,ALFA,U,NU,IERR)
C * Program to solve integral equations
C * of the first type by Tikhonov's method with
C * generalized residual principle.
C * The conjugation method algorithm is used
C * for functional minimization
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION U0(M),U(NU),Z(N)
0005 EXTERNAL AK
0006 EXTERNAL PTICR0,PTIZR1,PTIZRA
C * working array mapping
C * name length what
C * A N*M operator matrix
C * GRD N gradient in B PTIZR1
C * S N descent direction in PTIZR1
C * P M working array for PTIZR1
0007 ICONT=0
C * ICONT - start/continue mode flag
C * ICONT=0 - start mode
C * ICONT=1 - continue mode
0008 110 CONTINUE
0009 DU=SQRT(DL)
0010 DH=SQRT(H)
C * forming array starts
0011 NA=1
0012 NGRD=NA+N*M
0013 NS=NGRD+N
0014 NP=NS+N
0015 NU1=NP+M
0016 IF(NU1-1.GT.NU) GOTO 64
C * K1,K2 - iterations counters
0017 K1=0
0018 K2=0
0019 N1=N+1
0020 HX=(B-A)/(N-1.)
0021 HY=(D-C)/(M-1.)
0022 DD=HX/HY
C * RO - derivative wight in norm in W21,
C * divied by grid step square
0023 RO=1./HX**2
C * EPS - accuracy for discrepancy equation
0024 EPS=(C1-1.)*DL
0025 IF(ICONT.EQ.1) GOTO 111
C * forming matrix operator
0026 CALL PTICR0(AK,U(NA),A,B,C,D,N,M)
0027 111 CONTINUE
C * getting initial ALFA so that
C * generalized discrepancy will be positive
0028 13 CONTINUE
C * minimizing functional by conjugate gradient method
0029 CALL PTIZR1(U(NA),Z,U0,N,M,ITER,
*0D0,0D0,N,AN2,ALFA*DD,RO,Z,
*U(NGRD),U(NS),U(NP),IER)
C * calculating generalized discrepancy
0030 CALL PTIZRA(AN2,Z,N,DU,DH,HX,HY,RO,AN4)
0031 IF(C1.LE.1.)GOTO 100
0032 IF(ALFA.EQ.0.)GOTO 68
0033 IF(AN4.GT.EPS) GOTO 11
C * if 0034 K1=K1+1
0035 IF(K1.EQ.IMAX) GOTO 65
0036 ALFA=2.*ALFA
0037 GOTO 13
C * setting two initial points for chord method
0038 11 CONTINUE
0039 F0=AN4
0040 X0=1./ALFA
0041 ALFA=ALFA*2.
0042 X=1./ALFA
C * conjugate gradient method
0043 CALL PTIZR1(U(NA),Z,U0,N,M,ITER,
*0D0,0D0,N,AN2,ALFA*DD,RO,Z,
*U(NGRD),U(NS),U(NP),IER)
C * calculating generalized discrepancy
0044 CALL PTIZRA(AN2,Z,N,DU,DH,HX,HY,RO,AN4)
0045 14 CONTINUE
C * if accuracy achieved exit the program
0046 IF(DABS(AN4).LE.EPS) GOTO 100
C * if discrepancy < -EPS starting modified chord method
0047 IF(AN4.LE.-EPS)GOTO 2
0048 IF(ALFA.EQ.0.)GOTO 68
0049 K2=K2+1
0050 IF(K2.EQ.IMAX)GOTO 66
C * chord method formulae
0051 Y=X0-F0/(AN4-F0)*(X-X0)
0052 X0=X
0053 X=Y
0054 F0=AN4
0055 ALFA=1./X
C * conjugate gradient method
0056 CALL PTIZR1(U(NA),Z,U0,N,M,ITER,
*0D0,0D0,N,AN2,ALFA*DD,RO,Z,
*U(NGRD),U(NS),U(NP),IER)
C * calculating generalized discrepancy
0057 CALL PTIZRA(AN2,Z,N,DU,DH,HX,HY,RO,AN4)
0058 GOTO 14
0059 2 CONTINUE
C * starting modified chord method
0060 F=AN4
0061 23 CONTINUE
0062 Y=X0+F*(X-X0)/(F-F0)
0063 ALFA=1./Y
C * conjugate gradient method
0064 CALL PTIZR1(U(NA),Z,U0,N,M,ITER,
*0D0,0D0,N,AN2,ALFA*DD,RO,Z,
*U(NGRD),U(NS),U(NP),IER)
C * calculating generalized discrepancy
0065 CALL PTIZRA(AN2,Z,N,DY,DH,HX,HY,RO,AN4)
C * if accuracy achieved exit the program
0066 IF(DABS(AN4).LE.EPS) GOTO 101
0067 IF(AN4.LE.-EPS) GOTO 37
0068 IF(ALFA.EQ.0.)GOTO 68
0069 K2=K2+1
0070 IF(K2.EQ.IMAX)GOTO 67
0071 X0=Y
0072 F0=AN4
0073 GOTO 23
0074 37 CONTINUE
C * changing interval
0075 X=Y
0076 F=AN4
0077 GOTO 23
0078 ENTRY PTIZRE
C * entry to continue calculations
0079 ICONT=1
0080 GOTO 110
0081 64 CONTINUE
C * working array is too short
0082 IERR=64
0083 GOTO 999
0084 65 CONTINUE
C * initial regalarization parameter is two small
0085 IERR=65
0086 GOTO 999
0087 66 CONTINUE
C * IMAX iterations of the chord method made
0088 IERR=66
0089 GOTO 999
0090 67 CONTINUE
C * IMAX iterations by the modified chord method made
0091 IERR=67
0092 GOTO 999
0093 68 CONTINUE
C * ALFA=0 is specified or found
0094 IERR=68
0095 GOTO 999
0096 100 CONTINUE
C * solution is found by chord method
0097 IERR=0
0098 GOTO 999
0099 101 CONTINUE
C * solution is found by modified chord method
0100 IERR=1
0101 999 CONTINUE
0102 RETURN
0103 END


0001 SUBROUTINE PTIZR1(A,Z0,U0,N,M,ITER,DL2,
*ANGRD,IMAX,AN2,ALF,RO,Z,GRAD,S,U,IERR)
C * minimizes Tikhonov's functional
C * by conjugate gradient method.
C * AK(M,N) - operator matrix in equation AZ=U
C * Z0 - start point
C * U0 - right-hand side of th equation AZ=U
C * DL2 - discrepancy level to stop iterations
C * ANGRD - gradient norm level
C * AN2 - discrepancy achieved
C * ALF - regularization parameter
C * RO - differential derivative wight
C * in stabilizer
C * solution - extremal in the array Z
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION A(M,N),Z0(N),U0(M),
*Z(N),GRAD(N),S(N),U(M)
0005 EXTERNAL PTICR1,PTICR4,PTICR8,
*PTICR6,PTICR3,PTICR5,PTICRO
0006 ITER=0
0007 CALL PTICR1(Z0,Z0,Z,N,0D0)
C * getting functional gradient
0008 CALL PTICR3(A,Z,U,N,M)
0009 CALL PTICR4(GRAD,U,U0,A,N,M)
0010 CALL PTICR8(GRAD,Z,N,ALF,RO)
C * getting conjugate direction S
0011 CALL PTICR1(GRAD,GRAD,S,N,0D0)
C * getting norm of S
0012 CALL PTICR6(S,S,N,RNORM1)
C * Attention! Machine constant
0013 ALM=1.E+18
C * solving one-dimensional minimization problem
C * getting descent step
0014 7 CONTINUE
0015 CALL PTICRO(A,Z,GRAD,U,S,ALM,BETA,
*N,M,ALF,RO,IED)
0016 IF(BETA.EQ.0.)GOTO 14
0017 BETA=-BETA
C * getting new approximation Z
0018 CALL PTICR1(Z,S,Z,N,BETA)
C * getting discrepancy
0019 CALL PTICR3(A,Z,U,N,M)
0020 CALL PTICR5(U,U0,M,AN2)
C * if accuracy achieved exit the program
0021 IF(AN2.LE.DL2)GOTO 10
0022 ITER=ITER+1
C * getting gradient
0023 CALL PTICR4(GRAD,U,U0,A,N,M)
0024 CALL PTICR8(GRAD,Z,N,ALF,RO)
C * getting gradient norm
0025 CALL PTICR6(GRAD,GRAD,N,RNORM)
0026 RNM=RNORM/RNORM1
C * getting conjugate direction
0027 CALL PTICR1(GRAD,S,S,N,RNM)
0028 RNORM1=RNORM
C * if gradient norm is too small, exit the program
0029 IF(RNORM1.LE.ANGRD) GOTO 11
0030 IF(ITER.GE.IMAX) GOTO 13
0031 GOTO 7
0032 10 CONTINUE
C * accuracy achieved
0033 IERR=0
0034 GOTO 999
0035 11 CONTINUE
C * gradient norm is small
0036 IERR=1
0037 GOTO 999
0038 13 CONTINUE
C * IMAX iterations made
0039 IERR=2
0040 GOTO 999
0041 14 CONTINUE
C * minimization step equals zero
0042 IERR=3
0043 999 CONTINUE
0044 RETURN
0045 END


0001 SUBROUTINE PTIZRA(AN2,Z,N,DL,DH,
*HX,HY,RO,AN4)
C * calculates generalized discrepancy
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION Z(N)
0005 EXTERNAL PTICR6
C * norm discrepancy
0006 AN2=AN2*HY
0007 AN4=AN2-DL**2
C * if H=0 we don't need to calculate solution norm
0008 IF(DH.EQ.0.)GOTO 999
C * calculating solution norm in W21
0009 CALL PTICR6(Z,Z,N,S)
0010 S=S*HX
C * S - solution norm square in L2
0011 S1=0.
0012 DO 1 I=2,N
0013 S1=S1+(Z(I)-Z(I-1))**2
0014 1 CONTINUE
0015 S1=S1*HX
C * S1 - solution derivative norm square in L2
0016 S=DSQRT(S+RO*S1)
C * S - solution norm in W21
0017 AN4=AN2-(DL+DH*S)**2
0018 999 CONTINUE
0019 RETURN
0020 END


Listing file ptip.for


0001 SUBROUTINE PTIPR(AK,U0,A,B,C,D,N,M,Z,
*IC,AN2,DL,H,C1,ANGRD,IMAX,
*ALFA,R,NR,IERR)
C * Program to solve integral equations
C * of the first kind by regularization method
C * with generalized principle of residual to get
C * regularization parameter.
C * To minimize smoothing functional the conjugate
C * gradient methos is used
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION U0(M),R(NR),Z(N)
0005 EXTERNAL AK
0006 EXTERNAL PTICR0,PTICR1,
*PTISR1,PTISR2,PTISR3
C * Working array mapping
C * name length what
C * A: N*M operator matrix
C * H: N descent direction
C * G: N gradient
C * U: M operator value
C * U1: M working array
C * S: N working array
C *
C * NR=N*M+3N+2M
0007 ICONT=0
C * ICONT - start/continue mode flag
C * ICONT=0 start mode
C * ICONT=1 continue mode
0008 110 CONTINUE
0009 IF(IC.NE.0.AND.IC.NE.1)GOTO 69
C * getting arrays starts
0010 NA=1
0011 NH=N*M+1
0012 NG=NH+N
0013 NU=NG+N
0014 NU1=NU+M
0015 NS=NU1+M
0016 NMAX=NS+N
0017 IF(NMAX-1.GT.NR)GOTO 64
0018 DU=DSQRT(DL)
0019 DH=DSQRT(H)
C * K1,K2 - iteration numbers
0020 K1=0
0021 K2=0
0022 N1=N+1
0023 HX=(B-A)/(N-1.)
0024 HY=(D-C)/(M-1.)
0025 DD=HX/HY
C * RO - derivative wight in the norm
C * divided by grid step square
0026 IF(IC.EQ.0)RO=1./HX**2
0027 IF(IC.EQ.1)RO=0.0
C * EPS - accuracy of residual equation solution
0028 EPS=(C1-1.)*DL
0029 IF(ICONT.EQ.1) GOTO 111
C * getting operator matrix A
0030 CALL PTICR0(AK,R(NA),A,B,C,D,N,M)
0031 111 CONTINUE
C * make transformation to П-plus
0032 CALL PTISR2(R(NA),Z,N,M,IC,R(NS))
0033 CALL PTICR1(Z,Z,R(NH),N,0D0)
C * finding regularization parameter with
C * positive generalized discrepancy
0034 13 CONTINUE
C * getting functional minimun by conjugate gradient method
0035 CALL PTISR1(R(NA),Z,U0,N,M,ITER,
*0D0,0D0,N,AN2,ALFA*DD,RO,Z,R(NU),
*R(NU1),R(NH),R(NG),R(NS),IER)
C * calculating generalized dicrepancy
0036 CALL PTIZRA(AN2,Z,N,DU,DH,HX,HY,RO,AN4)
0037 IF(C1.LE.1.)GOTO 100
0038 IF(ALFA.EQ.0.)GOTO 68
0039 IF(AN4.GT.EPS) GOTO 11
C * multiplay ALFA by 2 while discrepancy is less EPS
0040 K1=K1+1
0041 IF(K1.EQ.IMAX) GOTO 65
0042 ALFA=2.*ALFA
0043 GOTO 13
C * setting initial points of the chird method
0044 11 CONTINUE
0045 F0=AN4
0046 X0=1./ALFA
0047 ALFA=ALFA*2.
0048 X=1./ALFA
C * using conjugate gradien method to minimize functional
0049 CALL PTISR1(R(NA),Z,U0,N,M,ITER,
*0D0,0D0,N,AN2,ALFA*DD,RO,Z,R(NU),
*R(NU1),R(NH),R(NG),R(NS),IER)
C * getting generalized discrepancy
0050 CALL PTIZRA(AN2,Z,N,DU,DH,HX,HY,RO,AN4)
0051 14 CONTINUE
C * if accuracy achieved, exit the program
0052 IF(DABS(AN4).LE.EPS) GOTO 100
C * if discrepancy < 0, switch to modified chord method
0053 IF(AN4.LE.-EPS)GOTO 2
0054 IF(ALFA.EQ.0.)GOTO 68
0055 K2=K2+1
0056 IF(K2.EQ.IMAX)GOTO 66
C * chord method formulae
0057 Y=X0-F0/(AN4-F0)*(X-X0)
0058 X0=X
0059 X=Y
0060 F0=AN4
0061 ALFA=1./X
C * using conjugate gradient method
0062 CALL PTISR1(R(NA),Z,U0,N,M,ITER,
*0D0,0D0,N,AN2,ALFA*DD,RO,Z,R(NU),
*R(NU1),R(NH),R(NG),R(NS),IER)
C * calculating generalized discrepancy
0063 CALL PTIZRA(AN2,Z,N,DU,DH,HX,HY,RO,AN4)
0064 GOTO 14
0065 2 CONTINUE
C * starting modified chord method
0066 F=AN4
0067 23 CONTINUE
0068 Y=X0+F*(X-X0)/(F-F0)
0069 ALFA=1./Y
C * using conjugate gradient method
0070 CALL PTISR1(R(NA),Z,U0,N,M,ITER,
*0D0,0D0,N,AN2,ALFA*DD,RO,Z,R(NU),
*R(NU1),R(NH),R(NG),R(NS),IER)
C * getting generalized discrepancy
0071 CALL PTIZRA(AN2,Z,N,DY,DH,HX,HY,RO,AN4)
C * if accuracy achieved, exit the problem
0072 IF(DABS(AN4).LE.EPS) GOTO 101
0073 IF(AN4.LE.-EPS) GOTO 37
0074 IF(ALFA.EQ.0.)GOTO 68
0075 K2=K2+1
0076 IF(K2.EQ.IMAX)GOTO 67
0077 X0=Y
0078 F0=AN4
0079 GOTO 23
0080 37 CONTINUE
C * changing the interval
0081 X=Y
0082 F=AN4
0083 GOTO 23
0084 ENTRY PTIPRE
C * entry to continue calculations
0085 ICONT=1
0086 GOTO 110
0087 64 CONTINUE
C * working array is too short
0088 IERR=64
0089 GOTO 9999
0090 65 CONTINUE
C * initial regularization parameter is too small
0091 IERR=65
0092 GOTO 999
0093 66 CONTINUE
C * IMAX iterations of chord method made
0094 IERR=66
0095 GOTO 999
0096 67 CONTINUE
C * IMAX iterations of modified chord method made
0097 IERR=67
0098 GOTO 999
0099 68 CONTINUE
C * ALFA=0 is set or found
0100 IERR=68
0101 GOTO 999
0102 69 CONTINUE
C * incorrect set type is specified
0103 IERR=69
0104 GOTO 9999
0105 100 CONTINUE
C * solution is found
0106 IERR=0
0107 GOTO 999
0108 101 CONTINUE
C * solution is found by modified chord method
0109 IERR=1
0110 999 CONTINUE
C * return from П-plus to original coordinates
0111 CALL PTICR1(Z,Z,R(NH),N,0D0)
0112 CALL PTISR3(Z,N,IC,R(NS))
0113 9999 CONTINUE
0114 RETURN
0115 END


0001 SUBROUTINE PTIZRA(AN2,Z,N,DL,DH,
*HX,HY,RO,AN4)
C * calculates generalized discrepancy
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION Z(N)
0005 EXTERNAL PTICR6
C * norm discrepancy
0006 AN2=AN2*HY
0007 AN4=AN2-DL**2
C * if H=0 we don't need to calculate solution norm
0008 IF(DH.EQ.0.)GOTO 999
C * calculating solution norm in W21
0009 CALL PTICR6(Z,Z,N,S)
0010 S=S*HX
C * S - solution norm square in L2
0011 S1=0.
0012 DO 1 I=2,N
0013 S1=S1+(Z(I)-Z(I-1))**2
0014 1 CONTINUE
0015 S1=S1*HX
C * S1 - solution derivative norm square in L2
0016 S=DSQRT(S+RO*S1)
C * S - solution norm in W21
0017 AN4=AN2-(DL+DH*S)**2
0018 999 CONTINUE
0019 RETURN
0020 END




0001 SUBROUTINE PTIKR(AK,U0,A,B,C,D,
*L1,L2,N,Z,AN,DL,HH,C1,IMAX,ALPHA,
*U,NU,IERR)
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 REAL*8 L1,L2
0005 DIMENSION U(NU),U0(N),Z(N)
0006 EXTERNAL AK
0007 EXTERNAL PTICR1,PTICR2,PTIKR1
C * working array mapping
C * name length what it is
C * ARE N real part of kernel Fourier image
C * AIM N imagenary part of kernel Fourier image
C * URE N real part of right side Fourier image
C * UIM N imagenary part of right side Fourier image
C * ZIM N imagenary part of solution
C * W N stabilizator
0008 IP=1
C * IP - mode:
C * on the first call IP=+1
C * on the next calls through PTIKRE IP=-1
0009 EPRO=0.
0010 IF(C1.GT.1.)EPRO=(C1-1.)*DL
C * EPRO - accuracy for discrepancy equation
C * if EPRO=0, the only extremal is calculated with the
C * fixed regularization parameter ALFA
0011 100 CONTINUE
C * calculating support of solution
0012 A=C-.5*(L1+L2)
0013 B=D-.5*(L1+L2)
C * T - period
0014 T=D-C
C * getting arrays starts
0015 NAR=1
0016 NAI=NAR+N
0017 NW=NAI+N
0018 NZI=NW+N
0019 NUR=NZI+N
0020 NUI=NUR+N
0021 NMAX=NUI+N
C * check if the working array size if sufficient
0022 IF(NMAX-1.GT.NU) GO TO 64
0023 IF(IP.EQ.-1) GO TO 101
C * setting right side
0024 CALL PTICR1(U0,U0,U(NUR),N,0D0)
0025 CALL PTICR2(U(NUI),0D0,N)
0026 101 CONTINUE
C * solving equation
0027 CALL PTIKR1(AK,U(NAR),U(NAI),Z,
*U(NZI),U(NUR),U(NUI),U(NW),N,ALPHA,
*L1,L2,AN,OM,T,DSQRT(DL),DSQRT(HH),
*IP,EPRO,IQ,IMAX,IERR)
0028 GO TO 999
0029 64 CONTINUE
C * working array is too short
0030 IERR=64
0031 GO TO 999
C * the entry PTIKRE is for the repeated solving
C * of the same equation and right side
0032 ENTRY PTIKRE
0033 IP=-1
0034 GO TO 100
0035 999 CONTINUE
0036 RETURN
0037 END


0001 SUBROUTINE PTIKR1(AK,ARE,AIM,ZRE,ZIM,
*URE,UIM,W,N,ALP,L1,L2,BETA,RO,A,
*DEL,HH,IPAR,EPRO,IQ,IMAX,IERR)
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 REAL*8 L1,L2
0005 DIMENSION ARE(N),ZRE(N),URE(N),
*AIM(N),ZIM(N),UIM(N),W(N)
0006 EXTERNAL FTF1C
0007 H=A/N
0008 HA=H/N
C * checking if it is the first call
0009 IF(IPAR.EQ.-1) GO TO 2
C * getting equation kernel
0010 DO 1 I=1,N
0011 ARG=(I-N/2-1)*H+0.5*(L1+L2)
0012 ARE(I)=AK(ARG)*H
0013 AIM(I)=0.0
0014 IF(ARG.LT.L1.OR.ARG.GT.L2) ARE(I)=0.0
0015 W(I)=(2.0/H*DSIN(3.14159265358D0/
/N*(I-1)))**2
0016 1 CONTINUE
0017 P=1.0
C * Fourier transformation of the kernel and right side
0018 CALL FTF1C(ARE,AIM,N,1,1,P)
0019 CALL FTF1C(URE,UIM,N,1,1,P)
0020 2 CONTINUE
C * getting ALP so, that RO(ALP) > 0
C * IQ - iteration number
0021 IQ=0
0022 77 CONTINUE
0023 IQ=IQ+1
C * discrepancy calculation
0024 F1=0.0
0025 F2=0.0
0026 DO 44 M=1,N
0027 U2=URE(M)**2+UIM(M)**2
0028 A2=ARE(M)**2+AIM(M)**2
0029 X=1.0+W(M)
0030 IF(A2.EQ.0.0) GO TO 42
0031 BA=X/(A2+ALP*X)
0032 AB=1.0-ALP*BA
0033 C1=HA*U2*(BA*ALP)**2
0034 C2=HA*U2*AB*BA
0035 F1=F1+C1
0036 F2=F2+C2
0037 GO TO 44
0038 42 F1=F1+HA*U2
0039 44 CONTINUE
C * calculating generalized discrepancy
C * to check if RO(ALP)>0
0040 BETA=F1
0041 RO=F1-(DEL+HH*DSQRT(F2))**2
C * if EPRO=0 is set then finish
0042 IF(EPRO.EQ.0.0) GO TO 10
C * if ALP=0.0 is set then finish
0043 IF(ALP.EQ.0.0) GO TO 68
0044 IF(RO.GT.0.0) GO TO 33
C * go to solution calculation if there is no
C * parameter ALP with RO(ALP) > 0
0045 IF(IQ.GT.IMAX) GO TO 65
0046 ALP=2.0*ALP
0047 GO TO 77
C * starting Newton method calculations
0048 33 CONTINUE
0049 IQ=0
0050 3 CONTINUE
0051 IQ=IQ+1
C * discrepancy calculations
0052 F1=0.0
0053 F2=0.0
0054 F3=0.0
0055 DO 4 M=1,N
0056 U2=URE(M)**2+UIM(M)**2
0057 A2=ARE(M)**2+AIM(M)**2
0058 X=1.0+W(M)
0059 IF(A2.EQ.0.0) GO TO 41
0060 BA=X/(A2+ALP*X)
0061 AB=1.0-ALP*BA
0062 C1=HA*U2*(BA*ALP)**2
0063 C2=HA*U2*AB*BA
0064 C3=2.0*C1*AB
0065 F1=F1+C1
0066 F2=F2+C2
0067 F3=F3+C3
0068 GO TO 4
0069 41 F1=F1+HA*U2
0070 4 CONTINUE
C * calculation discrepancy -BETA,
C * solution norm in W21 -ZNOR,
C * generalized discrepancy - RO
C * and its derivative - DR
0071 BETA=F1
0072 ZNOR=DSQRT(F2)
0073 RO=BETA-(DEL+HH*ZNOR)**2
0074 IF(ALP.EQ.0.0) GO TO 68
0075 DR=-F3*ALP-(DEL+HH*ZNOR)*HH*F3/ZNOR
C * calculate solution if accuracy acheaved
0076 IF(DABS(RO).LT.EPRO) GO TO 10
C * starting chord method if generalized discrepancy
C* is not convex
0077 IF(RO.LT.0.0) GO TO 61
C * getting new argument by Newton method
0078 DQ=-RO/DR
0079 A1=ALP
0080 R1=RO
0081 ALP=1.0/(1.0/ALP+DQ)
0082 IF(IQ.GE.IMAX) GO TO 66
0083 GO TO 3
C * chord method
0084 61 CONTINUE
C * change interval if RO<0
0085 6 CONTINUE
0086 A0=ALP
0087 R0=RO
0088 7 CONTINUE
0089 IQ=IQ+1
C * getting new argument by chord method
0090 ALP=A0*A1*(R0-R1)/(A0*R0-A1*R1)
C * calculation discrepancies
0091 8 CONTINUE
0092 F1=0.0
0093 F2=0.0
0094 DO 9 M=1,N
0095 U2=URE(M)**2+UIM(M)**2
0096 A2=ARE(M)**2+AIM(M)**2
0097 X=1.0+W(M)
0098 IF(A2.EQ.0.0) GO TO 91
0099 BA=X/(A2+ALP*X)
0100 AB=1.0-ALP*BA
0101 C1=HA*U2*(BA*ALP)**2
0102 C2=HA*U2*AB*BA
0103 F1=F1+C1
0104 F2=F2+C2
0105 GO TO 9
0106 91 F1=F1+HA*U2
0107 9 CONTINUE
C * generalized discrepancy and solution norm
0108 BETA=F1
0109 ZNOR=DSQRT(F2)
0110 RO=BETA-(DEL+HH*ZNOR)**2
0111 IF(ALP.EQ.0.0) GO TO 68
C * calculation solution if accuracy achieved
0112 IF(DABS(RO).LT.EPRO) GO TO 11
0113 IF(IQ.EQ.IMAX) GO TO 67
0114 IF(RO.LT.0.0) GO TO 6
C * change interval if RO>0
0115 A1=ALP
0116 R1=RO
0117 GO TO 7
0118 65 CONTINUE
C * RO(ALP) <= 0 for all regularization parameters
0119 IERR=65
0120 GO TO 999
0121 66 CONTINUE
C * IMAX iterations by Newton method performed
0122 IERR=66
0123 GO TO 999
0124 67 CONTINUE
C * IMAX iterations by chord method performed
0125 IERR=67
0126 GO TO 999
0127 68 CONTINUE
C * ALP=0.0 is set or found
0128 IERR=68
0129 GO TO 999
0130 11 CONTINUE
C * solution is found by means of chord method
0131 IERR=1
0132 GO TO 999
0133 10 CONTINUE
C * solution is found by Newton method
0134 IERR=0
0135 999 CONTINUE
C * calculation Fourier image of solution
0136 SSI=-1.0
0137 DO 12 M=1,N
0138 SSI=-SSI
0139 ZZ=N*(ARE(M)**2+AIM(M)**2+ALP*(1.0+W(M)))
0140 IF(ZZ.NE.0.0) GO TO 111
0141 ZRE(M)=0.0
0142 ZIM(M)=0.0
0143 GO TO 12
0144 111 ZRE(M)=SSI*(ARE(M)*URE(M)+AIM(M)*UIM(M))/ZZ
0145 ZIM(M)=SSI*(ARE(M)*UIM(M)-AIM(M)*URE(M))/ZZ
0146 12 CONTINUE
C * inverse Fourier transformation of solution
0147 P=-1.0
0148 CALL FTF1C(ZRE,ZIM,N,1,1,P)
0149 RETURN
0150 END


0001 SUBROUTINE FTF1C(ARE,AIM,N,IN,K,P)
C * Fast Fourier Transformation
C * ARE(N) - real part of function and its image
C * AIM(N) - imagenary part of function and its image
C * if P>0 - direct Fourier transformation
C * else - inverse Fourier transformation
C * formal parameter is to be called by name
C * N - the power of 2
C * start point - IN
C * step - K
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION ARE(N),AIM(N)
0005 N1=N/2
0006 MM=N1/2
0007 N2=N1+K
0008 J=IN
0009 JJ=J
0010 1 J=J+K
0011 IF(J-N1)2,2,10
0012 2 II=JJ+N1
0013 R=ARE(J)
0014 ARE(J)=ARE(II)
0015 ARE(II)=R
0016 R=AIM(J)
0017 AIM(J)=AIM(II)
0018 AIM(II)=R
0019 J=J+K
0020 M=MM
0021 3 IF(JJ-M)5,5,4
0022 4 JJ=JJ-M
0023 M=M/2
0024 GO TO 3
0025 5 JJ=JJ+M
0026 IF(JJ-J)1,1,6
0027 6 R=ARE(J)
0028 ARE(J)=ARE(JJ)
0029 ARE(JJ)=R
0030 R=AIM(J)
0031 AIM(J)=AIM(JJ)
0032 AIM(JJ)=R
0033 I=J+N2
0034 II=JJ+N2
0035 R=ARE(I)
0036 ARE(I)=ARE(II)
0037 ARE(II)=R
0038 R=AIM(I)
0039 AIM(I)=AIM(II)
0040 AIM(II)=R
0041 GO TO 1
0042 10 I=K
0043 T=3.14159265359
0044 IF(P)13,17,11
0045 11 T=-T
0046 13 P=-T
0047 14 SI=0.
0048 CO=1.
0049 S=DSIN(T)
0050 C=DCOS(T)
0051 T=0.5*T
0052 II=I
0053 I=I+I
0054 DO 16 M=IN,II,K
0055 DO 15 J=M,N,I
0056 JJ=J+II
0057 A=ARE(JJ)
0058 B=AIM(JJ)
0059 R=A*CO-B*SI
0060 ARE(JJ)=ARE(J)-R
0061 ARE(J)=ARE(J)+R
0062 R=B*CO+A*SI
0063 AIM(JJ)=AIM(J)-R
0064 15 AIM(J)=AIM(J)+R
0065 R=C*CO-S*SI
0066 SI=C*SI+S*CO
0067 16 CO=R
0068 IF(I-N)14,17,17
0069 17 RETURN
0070 END


0001 SUBROUTINE PTITR(AK,U0,ALIM,N1,N2,Z,
*DL,HH,C1,ALPHA,AN,U,NU,IMAX,IERR)
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION U(NU),U0(N1,N2)
0005 DIMENSION Z(N1,N2),ALIM(12)
0006 EXTERNAL AK
C * working array mapping
C * name length what
C * ARE N1*N2 real part of kernel Fourier image
C * AMIN N1*N2 imaginary part of kernel Fourier image
C * URE N1*N2 real part of right side Fourier image
C * UIM N1*N2 imaginary part of right side F. image
C * ZIM N1*N2 imaginary part of solition
C * W1 И W2 N1 И N2 stbilizer
0007 IP=1
C * IP - start/continue mode flag
C * start mode IP=+1
C * reenter mode through PTIKRE IP=-1
0008 EPRO=0.
0009 IF(C1.GT.1.)EPRO=(C1-1.)*DL
C * EPRO - discrepancy equation accuracy
C * if EPRO=0, the only extremal with the specified ALFA
0010 100 CONTINUE
C * calculating solution support
0011 ALIM(1)=ALIM(5)-.5*(ALIM(9)+ALIM(10))
0012 ALIM(2)=ALIM(6)-.5*(ALIM(9)+ALIM(10))
0013 ALIM(3)=ALIM(7)-.5*(ALIM(11)+ALIM(12))
0014 ALIM(4)=ALIM(8)-.5*(ALIM(11)+ALIM(12))
C * T1, T2 - periods
0015 T1=ALIM(6)-ALIM(5)
0016 T2=ALIM(8)-ALIM(7)
C * arrays starts formimg
0017 NAR=1
0018 NQU=N1*N2
0019 NAI=NAR+NQU
0020 NW1=NAI+NQU
0021 NW2=NW1+N1
0022 NZI=NW2+N2
0023 NUR=NZI+NQU
0024 NUI=NUR+NQU
0025 NMAX=NUI+NQU
C * check the length of working array
0026 IF(NMAX-1.GT.NU) GO TO 64
0027 IF(IP.EQ.-1) GO TO 101
C * create right-hand side
0028 CALL PTICR1(U0,U0,U(NUR),NQU,0D0)
0029 CALL PTICR2(U(NUI),0D0,NQU)
0030 101 CONTINUE
C * solving the equation
0031 CALL PTIKR3(AK,U(NAR),U(NAI),Z,U(NZI),
*U(NUR),U(NUI),U(NW1),U(NW2),N1,N2,
*ALPHA,ALIM,AN,OM,T1,T2,DSQRT(DL),
*DSQRT(HH),IP,EPRO,IQ,IMAX,IERR)
0032 GOTO 999
0033 64 CONTINUE
C * working array is too short
0034 IERR=64
0035 GOTO 999
C * PTITRE entry for the reentering to solve
C * equation with the same kernel and right-hand side
0036 ENTRY PTITRE
0037 IP=-1
0038 GOTO 100
0039 999 CONTINUE
0040 RETURN
0041 END


0001 SUBROUTINE PTIKR3(AK,ARE,AIM,ZRE,ZIM,
*URE,UIM,W1,W2,N1,N2,ALP,ALIM,BETA,RO,
*T1,T2,DEL,HH,IPAR,EPRO,IQ,IMAX,IERR)
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION ARE(N1,N2),ZRE(N1,N2),
*URE(N1,N2),AIM(N1,N2),ZIM(N1,N2),
*UIM(N1,N2),W1(N1),W2(N2),ALIM(12)
0005 H1=T1/N1
0006 H2=T2/N2
0007 HA=H1/N1*H2/N2
C * check if it is the first entry
0008 IF(IPAR.EQ.-1) GOTO 2
C * setting equation kernel
0009 DO 1 I=1,N1
0010 DO 1 J=1,N2
0011 ARG1=(I-N1/2-1)*H1+0.5*(ALIM(9)+ALIM(10))
0012 ARG2=(J-N2/2-1)*H2+0.5*(ALIM(11)+ALIM(12))
0013 ARE(I,J)=AK(ARG1,ARG2)*H1*H2
0014 AIM(I,J)=0.
0015 IF(ARG1.LT.ALIM(9).OR.ARG1.GT.ALIM(10).
.OR.ARG2.LT.ALIM(11).OR.ARG2.GT.ALIM(12))
* ARE(I,J)=0.
0016 W1(I)=(2.0/H1*DSIN(3.14159265359D0/N1*(I-1)))**2
0017 W2(J)=(2.0/H2*DSIN(3.14159265359D0/N2*(J-1)))**2
0018 1 CONTINUE
C * Fourier transformation of the kernel and right side
0019 P=1.0
0020 CALL FTFTC(ARE,AIM,N1,N2,P)
0021 P=1.0
0022 CALL FTFTC(URE,UIM,N1,N2,P)
0023 2 CONTINUE
C * find ALP so PO(ALP) > 0
0024 IQ=0
0025 77 CONTINUE
0026 IQ=IQ+1
C * calculating discrepancies
0027 F1=0.
0028 F2=0.
0029 DO 44 M=1,N1
0030 DO 44 N=1,N2
0031 U2=URE(M,N)**2+UIM(M,N)**2
0032 A2=ARE(M,N)**2+AIM(M,N)**2
0033 X=1.+(W1(M)+W2(N))**2
0034 IF(A2.EQ.0.) GOTO 42
0035 BA=X/(A2+ALP*X)
0036 AB=1.-ALP*BA
0037 C1=HA*U2*(BA*ALP)**2
0038 C2=HA*U2*AB*BA
0039 F1=F1+C1
0040 F2=F2+C2
0041 GOTO 44
0042 42 F1=F1+HA*U2
0043 44 CONTINUE
C * calculating generalized discrepancy
C * to check if RO(ALP) > 0
0044 BETA=F1
0045 RO=F1-(DEL+HH*DSQRT(F2))**2
C * if EPRO=0 exit the program
0046 IF(EPRO.EQ.0.) GOTO 10
C * if ALP=0.0 then exit the program
0047 IF(ALP.EQ.0.) GOTO 68
0048 IF(RO.GT.0.) GOTO 33
C * cannot find initial regularization parameter
C * RO(ALP) > 0
0049 IF(IQ.GT.IMAX) GOTO 65
0050 ALP=2.*ALP
0051 GOTO 77
C * starting the Newton method
0052 33 CONTINUE
0053 IQ=0
0054 3 CONTINUE
0055 IQ=IQ+1
C * discrepancies calculation
0056 F1=0.
0057 F2=0.
0058 F3=0.
0059 DO 4 M=1,N1
0060 DO 4 N=1,N2
0061 U2=URE(M,N)**2+UIM(M,N)**2
0062 A2=ARE(M,N)**2+AIM(M,N)**2
0063 X=1.+(W1(M)+W2(N))**2
0064 IF(A2.EQ.0.) GOTO 41
0065 BA=X/(A2+ALP*X)
0066 AB=1.-ALP*BA
0067 C1=HA*U2*(BA*ALP)**2
0068 C2=HA*U2*AB*BA
0069 C3=2.*C1*AB
0070 F1=F1+C1
0071 F2=F2+C2
0072 F3=F3+C3
0073 GOTO 4
0074 41 F1=F1+HA*U2
0075 4 CONTINUE
C * discrepancy -ВЕТА
C * solution norm in W21 -ZNOR,
C * generalized discrepancy -RO
C * and its derivative - DR
0076 BETA=F1
0077 ZNOR=DSQRT(F2)
0078 RO=BETA-(DEL+HH*ZNOR)**2
0079 IF(ALP.EQ.0.) GOTO 68
0080 DR=-F3*ALP-(DEL+HH*ZNOR)*HH*F3/ZNOR
C * accuracy achieved
0081 IF(DABS(RO).LT.EPRO) GOTO 10
C * switch to the chord method
0082 IF(RO.LT.0.) GOTO 61
C * calculating new argument
C * Newton method formulae
0083 DQ=-RO/DR
0084 A1=ALP
0085 R1=RO
0086 ALP=1./(1./ALP+DQ)
0087 IF(IQ.GE.IMAX) GOTO 66
0088 GOTO 3
C * chord method
0089 61 CONTINUE
C * changing interval if R0 < 0
0090 6 CONTINUE
0091 A0=ALP
0092 R0=RO
0093 7 CONTINUE
0094 IQ=IQ+1
C * getting new argument by chord method
0095 ALP=A0*A1*(R0-R1)/(A0*R0-A1*R1)
C * discrepancy calculations
0096 F1=0.
0097 F2=0.
0098 DO 9 M=1,N1
0099 DO 9 N=1,N2
0100 U2=URE(M,N)**2+UIM(M,N)**2
0101 A2=ARE(M,N)**2+AIM(M,N)**2
0102 X=1.+(W1(M)+W2(N))**2
0103 IF(A2.EQ.0.) GOTO 91
0104 BA=X/(A2+ALP*X)
0105 AB=1.-ALP*BA
0106 C1=HA*U2*(BA*ALP)**2
0107 C2=HA*U2*AB*BA
0108 F1=F1+C1
0109 F2=F2+C2
0110 GOTO 9
0111 91 F1=F1+HA*U2
0112 9 CONTINUE
C * getting discrepancy, generalized discrepancy
C * and solution norm
0113 BETA=F1
0114 ZNOR=DSQRT(F2)
0115 RO=BETA-(DEL+HH*ZNOR)**2
0116 IF(ALP.EQ.0.) GOTO 68
C * accuracy achieved
0117 IF(DABS(RO).LT.EPRO) GOTO 11
0118 IF(IQ.EQ.IMAX) GOTO 67
0119 IF(RO.LT.0.) GOTO 6
C * changing interval if RO > 0
0120 A1=ALP
0121 R1=RO
0122 GOTO 7
0123 65 CONTINUE
C * cannot find ALP so that RO(ALP) be greater zero
0124 IERR=65
0125 GOTO 999
0126 66 CONTINUE
C * IMAX Newton iterations made
0127 IERR=66
0128 GOTO 999
0129 67 CONTINUE
C * IMAX chord iterations made
0130 IERR=67
0131 GOTO 999
0132 68 CONTINUE
C * ALP=0.0 is specified or found
0133 IERR=68
0134 GOTO 999
0135 11 CONTINUE
C * solution is found by chord method
0136 IERR=1
0137 GOTO 999
0138 10 CONTINUE
C * solution is found by Newton method
0139 IERR=0
0140 999 CONTINUE
C * getting solution Fourier image
0141 SSI=1.0
0142 DO 12 M=1,N1
0143 SSI=-SSI
0144 DO 12 N=1,N2
0145 SSI=-SSI
0146 ZZ=N1*N2*(ARE(M,N)**2+AIM(M,N)**2+
+ALP*(1.0+(W1(M)+W2(N))**2))
0147 IF(ZZ.NE.0.) GOTO 111
0148 ZRE(M,N)=0.
0149 ZIM(M,N)=0.
0150 GOTO 12
0151 111 ZRE(M,N)=SSI*(ARE(M,N)*URE(M,N)+
+AIM(M,N)*UIM(M,N))/ZZ
0152 ZIM(M,N)=SSI*(ARE(M,N)*UIM(M,N)-
-AIM(M,N)*URE(M,N))/ZZ
0153 12 CONTINUE
C * inverse Fourier transformation
0154 P=-1.0
0155 CALL FTFTC(ZRE,ZIM,N1,N2,P)
0156 RETURN
0157 END


0001 SUBROUTINE FTFTC(ARE,AIM,N1,N2,P)
C * Fast two-dimensional Fourier transformation
C * if P > 0 - direct transformation
C * if P < 0 - inverse transformation
C * P - must be called by name
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION ARE(N1,N2),AIM(N1,N2)
0005 DO 1 I=1,N1
0006 CALL FTF1C(ARE,AIM,N1*N2,I,N1,P)
0007 1 CONTINUE
0008 DO 2 J=1,N2
0009 CALL FTF1C(ARE(1,J),AIM(1,J),N1,1,1,P)
0010 2 CONTINUE
0011 RETURN
0012 END


0001 SUBROUTINE FTF1C(XARE,XAIM,N,IN,K,P)
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION XARE(N),XAIM(N)
0005 N1=N/2
0006 MM=N1/2
0007 N2=N1+K
0008 J=IN
0009 JJ=J
0010 1 J=J+K
0011 IF(J-N1)2,2,10
0012 2 II=JJ+N1
0013 R=XARE(J)
0014 XARE(J)=XARE(II)
0015 XARE(II)=R
0016 R=XAIM(J)
0017 XAIM(J)=XAIM(II)
0018 XAIM(II)=R
0019 J=J+K
0020 M=MM
0021 3 IF(JJ-M)5,5,4
0022 4 JJ=JJ-M
0023 M=M/2
0024 GOTO 3
0025 5 JJ=JJ+M
0026 IF(JJ-J)1,1,6
0027 6 R=XARE(J)
0028 XARE(J)=XARE(JJ)
0029 XARE(JJ)=R
0030 R=XAIM(J)
0031 XAIM(J)=XAIM(JJ)
0032 XAIM(JJ)=R
0033 I=J+N2
0034 II=JJ+N2
0035 R=XARE(I)
0036 XARE(I)=XARE(II)
0037 XARE(II)=R
0038 R=XAIM(I)
0039 XAIM(I)=XAIM(II)
0040 XAIM(II)=R
0041 GOTO 1
0042 10 I=K
0043 T=3.14159265359
0044 IF(P)13,17,11
0045 11 T=-T
0046 13 P=-T
0047 14 SI=0.
0048 CO=1.
0049 S=DSIN(T)
0050 C=DCOS(T)
0051 T=0.5*T
0052 II=I
0053 I=I+I
0054 DO 16 M=IN,II,K
0055 DO 15 J=M,N,I
0056 JJ=J+II
0057 A=XARE(JJ)
0058 B=XAIM(JJ)
0059 R=A*CO-B*SI
0060 XARE(JJ)=XARE(J)-R
0061 XARE(J)=XARE(J)+R
0062 R=B*CO+A*SI
0063 XAIM(JJ)=XAIM(J)-R
0064 15 XAIM(J)=XAIM(J)+R
0065 R=C*CO-S*SI
0066 SI=C*SI+S*CO
0067 16 CO=R
0068 IF(I-N)14,17,17
0069 17 RETURN
0070 END


0001 SUBROUTINE PTIGR(AK,U0,X1,X2,Y1,Y2,N,M,
*Z,AN2,ITER,DL,IMAX,C1,C2,IC,R,NR,IERR)
C * Solving integral equations of the first kind
C * using the conditional gradient method
C * DISPATCHER
C * AK - kernel function
C * Working array mapping
C * name length what it is
C * A: N*M oparator matrix
C * H: N descent direction
C * G: N gradient
C * U: M operator value
C * U1: M working array
C *
C * NR=N*M+2(N+M)
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 EXTERNAL AK
0005 EXTERNAL PTICR0,PTIGR1
0006 DIMENSION U0(M),Z(N),R(NR)
C * get arrays indexes
0007 NA=1
0008 NH=N*M+1
0009 NG=NH+N
0010 NU=NG+N
0011 NU1=NU+M
0012 NMAX=NU1+M
0013 IF(NMAX-1.GT.NR)GOTO 64
C * getting operator matrix
0014 CALL PTICR0(AK,R(NA),X1,X2,Y1,Y2,N,M)
C * minimizing discrepancy
0015 CALL PTIGR1(R(NA),Z,U0,C1,C2,IC,N,M,
*ITER,DL*(M-1.)/(Y2-Y1),0D0,IMAX,AN2,
*Z,R(NU),R(NU1),R(NH),R(NG),IERR)
0016 AN2=AN2*(Y2-Y1)/(M-1.)
0017 GOTO 999
0018 64 CONTINUE
C * working array is too short
0019 IERR=64
0020 GOTO 999
0021 999 CONTINUE
0022 RETURN
0023 END


0001 SUBROUTINE PTIGR1(A,Z0,U0,C1,C2,IC,N,M,
*ITER,DL2,ANGRD,IMAX,AN2,
*Z,U,U1,H,G,IERR)
C * conditional gradient method for
C * discrepancy minimization
C * A(M,N) - operator matrix
C * Z0(N) - start approach
C * U0(M) - right side of equation
C * C1,C2 - constraints, restricting solution
C * IC - type of set of correctness
C * ITER - number of iterations made
C * DL2 - discrepancy level to stop iterations
C * ANGRD - gradient norm level to stop iterations
C * IMAX - maximum number of iteraqtions
C * Z(N) - solution found
C * U(M) - operator value on the solution
C * AN2 - discrepancy value on solution
C * G(N),H(N) - working arrays
C * U1(M) - working array for one-dimensional minimization
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 EXTERNAL PTIGR2,PTICR1,PTICR3
0005 EXTERNAL PTICR5,PTICR4,PTICR6
0006 DIMENSION U1(M),U(M),G(N),H(N),Z(N)
0007 DIMENSION A(M,N),Z0(N),U0(M)
0008 ITER=0
0009 ALMAX=1.
0010 CALL PTICR1(Z0,Z0,Z,N,0D0)
0011 CALL PTICR3(A,Z,U,N,M)
0012 CALL PTICR5(U,U0,M,AN2O)
0013 CALL PTICR4(G,U,U0,A,N,M)
0014 CALL PTICR6(G,G,N,ANGR)
0015 IF(AN2O.LE.DL2.OR.ANGR.LE.ANGRD.OR.
.ITER.GE.IMAX)GOTO 20
C * start of iterations
0016 14 CONTINUE
0017 ITER=ITER+1
0018 CALL PTIGR2(H,G,N,C1,C2,IC)
0019 CALL PTICR1(Z,H,H,N,-1D0)
0020 CALL PTICRO(A,Z,G,U1,H,ALMAX,AL,N,M,0.D0,0.D0,IED)
0021 CALL PTICR1(Z,H,Z,N,-AL)
0022 CALL PTICR3(A,Z,U,N,M)
0023 CALL PTICR4(G,U,U0,A,N,M)
0024 CALL PTICR6(G,G,N,ANGR)
0025 CALL PTICR5(U,U0,M,AN2)
0026 IF(AN2.LE.DL2.OR.ANGR.LE.ANGRD.OR.ITER.GE.IMAX)
+ GOTO 20
0027 IF(AN2.GE.AN2O)GOTO 21
C * storing discrepancy
0028 AN2O=AN2
0029 GOTO 14
0030 20 CONTINUE
C * normal end
0031 IERR=0
0032 GOTO 999
0033 21 CONTINUE
C * discrepancy hasn't decreased
0034 IERR=1
0035 999 CONTINUE
0036 RETURN
0037 END


0001 SUBROUTINE PTIGR2(TOP,G,N,C1,C2,IC)
C * finding optimal vertex
C * G - gradient array
C * TOP - vertex found, array
C * C1,C2 - restriction constants
C * IC=1 - monotoneous functions
C * IC=2 - monotoneous convex functions
C * IC=3 - convex functions
C * IC=-1,C1=0. - finit variation functions
C * Var F <= 2*C2
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 EXTERNAL PTICR1,PTICR2
0005 DIMENSION TOP(N),G(N)
0006 AL=0.
0007 S1=0.
0008 S2=0.
0009 B=1./(N-1.)
0010 A=1.-B
0011 SS=+1.
C * in K1 - number of the current optimal vertex
0012 K1=0
0013 ICM=IABS(IC)
0014 IF(ICM.EQ.1)GOTO 101
0015 DO 1 I=2,N
0016 S2=S2+G(I)*A
0017 1 A=A-B
C * in S2 - value of linear functional
C * without the left point
0018 101 CONTINUE
0019 DO 2 K=1,N
0020 SG=S1+S2+G(K)
C * SG - value of functional
0021 IF(AL.LE.SG)GOTO 3
0022 AL=SG
0023 K1=K
0024 SS=+1.
0025 3 CONTINUE
0026 IF(IC.GE.0)GOTO 8
C * checking symmetrical vertex, in need (IC<0)
0027 IF(AL.LE.-SG)GOTO 4
0028 AL=-SG
0029 K1=K
0030 SS=-1.
C * SS=-1 - setting the flag - symmetric value used
0031 4 CONTINUE
0032 8 CONTINUE
C * in S1 - first half of functional
C * and in S2 - the second
0033 IF(ICM.EQ.1.OR.ICM.EQ.2)S1=S1+G(K)
0034 IF(ICM.EQ.3)S1=(S1+G(K))*(K-1.)/K
0035 IF(K.GE.N-1.OR.ICM.EQ.1)GOTO 7
0036 S2=S2*(N-K)/(N-K-1.)-G(K+1)
0037 GOTO 6
0038 7 S2=0.
0039 6 CONTINUE
0040 2 CONTINUE
C * number of optimal vertex is found - K1
C * form this vertex vector
0041 IF(K1.NE.0)GOTO 9
0042 CALL PTICR2(TOP,C1,N)
0043 GOTO 999
0044 9 CONTINUE
0045 DO 5 I=1,N
0046 GOTO (801,802,803),ICM
0047 801 IF(I.LE.K1)TOP(I)=C2
0048 IF(I.GT.K1)TOP(I)=C1
0049 GOTO 5
0050 802 IF(I.LE.K1)TOP(I)=C2
0051 IF(I.GT.K1)
* TOP(I)=TOP(I-1)-(C2-C1)/(N-K1)
0052 GOTO 5
0053 803 IF(I.LE.K1.AND.K1.NE.1)
* TOP(I)=C1+(C2-C1)/(K1-1.)*(I-1.)*2.
0054 IF(I.GE.K1.AND.K1.NE.N)
* TOP(I)=C1+(C2-C1)/(N-K1)*(N-I+0.)*2.
0055 5 CONTINUE
0056 IF(IC.GE.0)GOTO 999
0057 CALL PTICR1(TOP,TOP,TOP,N,SS-1.)
0058 GOTO 999
0059 999 CONTINUE
0060 RETURN
0061 END

0001 SUBROUTINE PTILR(AK,U0,X1,X2,Y1,Y2,N,M,
*Z,AN2,DL,ITER,IMAX,C2,IC,R,NR,IERR)
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION U0(M),Z(N),R(NR)
0005 EXTERNAL AK
0006 EXTERNAL PTICR0,PTICR1,
*PTILRA,PTILRB,PTILR1,PTILR2
C * working array mapping
C name length what
C IM: 1 number of active restrictions for PTILR1(PTILR2)
C A: N*M operator matrix
C C: N*N (C*X,X)
C D: N (D,X)
C CON: NN*N restrictions matrix
C B: NN restrictions vector
C MASK: NN mask of active restrictions
C AJ: NN*N matrix of active restrictions
C P: NN*N working array
C PI: N*N projector (working)
C GR: N working array
C W: MAX(NN,N,M) working arrays
C P1: N working array
C *
C * NR>=3*N*M+3*N*N+8*N+MAX(M,N+1)+3
C *
C * ICONT=0 - first call
C * ICONT=1 - next call for continuing
0007 ICONT=0
0008 101 CONTINUE
C * number of restrictions NN
0009 NN=N
0010 IF(IC.EQ.1)NN=N+1
C * get array starts
0011 NA=2
0012 NC=NA+N*M
0013 ND=NC+N*N
0014 NCON=ND+N
0015 NB=NCON+N*NN
0016 NMASK=NB+NN
0017 NAJ=NMASK+NN
0018 NP=NAJ+NN*N
0019 NPI=NP+NN*N
0020 NGR=NPI+N*N
0021 NW=NGR+N
0022 NP1=NW+MAX0(NN,M,N)
0023 NMAX=NP1+N
0024 IF(NMAX-1.GT.NR)GOTO 64
0025 IF(ICONT.EQ.1)GOTO 100
C * if ICONT doesn't equal 1 getting operator matrix
0026 CALL PTICR0(AK,R(NA),X1,X2,Y1,Y2,N,M)
C * transforming discrepancy functional to the form
C * (C*X,X)+(D,X)+E
0027 CALL PTILRA(R(NA),U0,R(NC),R(ND),N,M)
C * getting restrictions matrix
0028 CALL PTILRB(R(NCON),R(NB),N,NN,IC,C2,IE)
0029 IF(IE.NE.0)GOTO 65
C * minimizing discrepancy
0030 CALL PTILR1(M,MR,R(NMASK),R(NA),R(NAJ),
*R(NCON),R(NP),U0,R(NPI),R(NB),R(NGR),
*R(NW),R(NP1),.TRUE.,IMAX,ITER,R(NC),Z,
*R(ND),N,NN,AN2,DL*(M-1.)/(Y2-Y1),
*0.0D0,IERR)
0031 102 CONTINUE
0032 AN2=AN2*(Y2-Y1)/(M-1.)
0033 GOTO 999
0034 ENTRY PTILRE
C * entry to continue
0035 ICONT=1
0036 GOTO 101
0037 100 CONTINUE
C * call continue program
0038 CALL PTILR2(M,MR,R(NMASK),R(NA),R(NAJ),
*R(NCON),R(NP),U0,R(NPI),R(NB),R(NGR),
*R(NW),R(NP1),.TRUE.,IMAX,ITER,R(NC),Z,
*R(ND),N,NN,AN2,DL*(M-1.)/(Y2-Y1),
*0.0D0,IERR)
0039 GOTO 102
0040 64 CONTINUE
C * working array is too short
0041 IERR=64
0042 GOTO 999
0043 65 CONTINUE
C * invalid correctness set
0044 IERR=65
0045 GOTO 999
0046 999 CONTINUE
0047 RETURN
0048 END


0001 SUBROUTINE PTILR5(A,AJ,B1,RMASK,N,M,NN,
*GR,PI,P,P1,C,X,D,K,IED,DGR)
C * minimization of quadratic functional
C * on the set defined by active restrictions AJ
C * numbers of active restrictions in RMASK,
C * start point - X, end point - X,
C * projector - PI.
C * IED - return code:
C * 0-new face achieved,
C * 1-minimum found
C * 2-gradient norm is too small
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION AJ(NN,N),B1(NN),RMASK(NN),
*GR(N),PI(N,N),P(N),P1(N),C(N,N),
*X(N),A(NN,N),D(N)
0005 EXTERNAL PTICR1,PTICR2,
*PTICR3,PTICR6,PTILR7
C * K - iteration count on the subspace
0006 K=0
0007 AN1=1.
C * on K=0 descent direction is 0.
0008 CALL PTICR2(P,0.0D0,N)
C * loop of the conjugate gradient method
C * check minimum by number of iterations
0009 11 IF(K.EQ.N-M) GOTO 14
C * get gradient, if its norm less then DGR then exit
0010 CALL PTILR7(X,C,N,GR,D)
0011 CALL PTICR6(GR,GR,N,R)
0012 IF(R.LT.DGR) GOTO 17
C * P1 - antigradient projection
0013 CALL PTICR3(PI,GR,P1,N,N)
0014 CALL PTICR1(P1,P1,P1,N,-2.0D0)
C * AN - projection norm
C * AN2- conjugation coeff.
0015 CALL PTICR6(P1,P1,N,AN)
0016 AN2=AN/AN1
C * P - new conjugate (descent) direction
0017 CALL PTICR1(P1,P,P,N,AN2)
C * R:=(P,GR)
0018 CALL PTICR6(P,GR,N,R)
C * R1:=(CP,P)
0019 R1=0.0
0020 DO 4 I=1,N
0021 R2=0.0
0022 DO 5 J=1,N
0023 R2=R2+C(I,J)*P(J)
0024 5 CONTINUE
0025 R1=R1+P(I)*R2
0026 4 CONTINUE
0027 IF(R1.EQ.0.0) GOTO 14
C * optimal step =-0.5R/R1
0028 ALPHA=-R/R1/2.
C * finding the nearest plane
0029 ALM=ALPHA+1.0
0030 NP=0
0031 701 DO 6 I=1,NN
C * don't testing active constrains
0032 IF(RMASK(I).EQ.1) GOTO 6
0033 R=0.0
0034 R1=0.0
0035 DO 7 J=1,N
0036 R=R+A(I,J)*X(J)
0037 R1=R1+A(I,J)*P(J)
0038 7 CONTINUE
C * R1 > 0
0039 IF(R1.LE.0.0) GOTO 6
C * ALP - distance to plane
0040 ALP=(B1(I)-R)/R1
C * finding minimal one
0041 IF(ALM.LE.ALP) GOTO 6
0042 ALM=ALP
C * NP - restriction number
0043 NP=I
0044 6 CONTINUE
c ALM=ALM*0.99999
C * if ALPHA > ALM - new constraint
0045 IF(ALPHA.GE.ALM) GOTO 12
C * if ALPHA < 0 - minimum is found,
C * else - to the next step
0046 IF(ALPHA.LE.0.0) GOTO 14
0047 CALL PTICR1(X,P,X,N,ALPHA)
0048 K=K+1
0049 AN1=AN
0050 GOTO 11
0051 12 ALPHA=ALM
C * going to the new point on the boundary
0052 CALL PTICR1(X,P,X,N,ALPHA)
C * including new constraint into the restriction mask
0053 RMASK(NP)=1
0054 M=M+1
0055 IED=0
0056 GOTO 16
0057 14 IED=1
0058 GOTO 16
0059 17 CONTINUE
0060 IED=2
0061 16 CONTINUE
0062 RETURN
0063 END


0001 SUBROUTINE PTILR7(X,C,N,GR,D)
C * calculating of the gradient of
C * the quadratic functional
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION C(N,N),X(N),GR(N),D(N)
C * GR:=2LX+D
0005 DO 1 I=1,N
0006 R=0.0
0007 DO 2 J=1,N
0008 R=R+(C(I,J)+C(J,I))*X(J)
0009 2 CONTINUE
0010 GR(I)=R+D(I)
0011 1 CONTINUE
0012 RETURN
0013 END


0001 SUBROUTINE PTILR4(AJ,P,PI,M,N,NN)
C * getting projector
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION AJ(NN,N),P(NN,N),PI(1)
0005 EXTERNAL PTILR0,PTICR2
0006 IF(M.EQ.0) GOTO 1
C * PI:=AJ*AJ'
0007 DO 2 I=1,M
0008 DO 2 J=1,I
0009 R=0.0
0010 DO 3 K=1,N
0011 R=R+AJ(I,K)*AJ(J,K)
0012 3 CONTINUE
0013 L=M*(I-1)+J
0014 PI(L)=R
0015 L=M*(J-1)+I
0016 PI(L)=R
0017 2 CONTINUE
C * PI:=INV(PI)
0018 CALL PTILR0(PI,M)
C * P:=-PI*AJ
0019 DO 4 I=1,M
0020 DO 4 J=1,N
0021 R=0.0
0022 DO 5 K=1,M
0023 L=M*(I-1)+K
0024 R=R+PI(L)*AJ(K,J)
0025 5 CONTINUE
0026 P(I,J)=-R
0027 4 CONTINUE
C * PI:=AJ'*P
0028 DO 6 I=1,N
0029 DO 6 J=1,I
0030 R=0.0
0031 DO 7 K=1,M
0032 R=R+AJ(K,I)*P(K,J)
0033 7 CONTINUE
0034 L=N*(I-1)+J
0035 PI(L)=R
0036 L=N*(J-1)+I
0037 PI(L)=R
0038 6 CONTINUE
0039 GOTO 10
0040 1 CONTINUE
C * if there are no active constraints then PI:=0, P:=0
0041 CALL PTICR2(PI,0.0D0,N*N)
0042 CALL PTICR2(P,0.0D0,N*NN)
C * PI:=E+PI
0043 10 DO 9 I=1,N
0044 L=N*(I-1)+I
0045 PI(L)=PI(L)+1.0
0046 9 CONTINUE
0047 RETURN
0048 END


0001 SUBROUTINE PTILR3(AJ,A,N,NN,RMASK)
C * geting active consraints by the mask
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION AJ(NN,N),A(NN,N),RMASK(NN)
0005 K=0
0006 DO 1 I=1,NN
0007 IF(RMASK(I).EQ.0) GOTO 1
0008 K=K+1
0009 DO 2 J=1,N
0010 AJ(K,J)=A(I,J)
0011 2 CONTINUE
0012 1 CONTINUE
0013 RETURN
0014 END


0001 SUBROUTINE PTILR6(A,X,P,RMASK,
*GR,N,M,NN,C,D,IED)
C * excluding active constraints
C * IED - return code:
C * 0-nothing to exclude
C * 1-something is excluded
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION P(NN,N),GR(N),RMASK(NN)
0005 DIMENSION A(NN,N),X(N),C(N,N),D(N)
0006 EXTERNAL PTILR7
0007 CALL PTILR7(X,C,N,GR,D)
0008 K=0
0009 IED=0
0010 AL=1.0
0011 DO 1 I=1,NN
0012 IF(RMASK(I).EQ.0) GOTO 1
C * getting shadow parameter for active constraints
0013 K=K+1
0014 R=0.0
0015 DO 2 J=1,N
0016 R=R+P(K,J)*GR(J)
0017 2 CONTINUE
C * getting minimal shadow parameter
0018 IF (R.GE.AL) GOTO 1
0019 AL=R
0020 NNN=I
0021 1 CONTINUE
C * if minimal shadow parameter > 0, then exit
C * else excluding NNN-th constraint
0022 IF (AL.GE.0.00000) GOTO 3
0023 RMASK(NNN)=0
0024 IED=1
0025 M=M-1
0026 3 RETURN
0027 END


0001 SUBROUTINE PTILR1(MN,M,RMASK,AQ,AJ,A,P,
*U,PI,B,GR,W,P1,WORK,ICM,ICI,C,X,D,N,NN,
*R,DEL,DGR,IEND)
C * dispatcher
C * WORK=.TRUE.- working with the matrix A being set
C * ICI - large loops counter
C * ICM - maximum number of large loops
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION A(NN,N),RMASK(NN),AJ(NN,N),
*P(NN,N),PI(N,N),B(NN),GR(N),W(MN),
*P1(N),C(N,N),X(N),D(N),AQ(MN,N),U(MN)
0005 LOGICAL WORK
0006 EXTERNAL PTICI2,PTICR3,PTICR5,
*PTILR3,PTILR4,PTILR5,PTILR6
0007 CALL PTICR2(RMASK,0.0D0,NN)
C * M - number of active constraints
0008 M=0
0009 ENTRY PTILR2
C * starting large loops iterations
0010 ICI=0
0011 2 CONTINUE
0012 IF(ICI.EQ.ICM) GOTO 12
0013 ICI=ICI+1
0014 IF(.NOT.WORK) GOTO 101
C * calculating discrepancy - R1,
C * if it is less < DEL then exit
0015 CALL PTICR3(AQ,X,W,N,MN)
0016 CALL PTICR5(W,U,MN,R1)
0017 101 CONTINUE
0018 IF(R1.LE.DEL) GOTO 9
C * getting constraints matrix by the mask
0019 CALL PTILR3(AJ,A,N,NN,RMASK)
C * getting projector
0020 CALL PTILR4(AJ,P,PI,M,N,NN)
C * minimizing on the face
0021 CALL PTILR5(A,AJ,B,RMASK,N,M,NN,GR,
*PI,W,P1,C,X,D,K,IED,DGR)
0022 IF(IED.EQ.0) GOTO 2
C * if gradient norm less DGR then exit
0023 IF(IED.EQ.2) GOTO 11
C * trying to exclude active constraint
0024 CALL PTILR6(A,X,P,RMASK,GR,N,M,NN,C,D,IED)
C * if something was excluded then continue calculations
0025 IF(IED.EQ.1) GOTO 2
C * normal return
C exact mimimum is found
0026 IEND=0
0027 GOTO 10
0028 9 CONTINUE
C * normal return
C * discrepancy level achieved
0029 IEND=1
0030 GOTO 10
0031 12 CONTINUE
C * maximun large loop number achieved
0032 IEND=3
0033 GOTO 10
0034 11 CONTINUE
C * gradient norm is too small
0035 IEND=2
0036 10 CONTINUE
0037 IF(.NOT.WORK)GOTO 102
C * calculating discrepancy
0038 CALL PTICR3(AQ,X,W,N,MN)
0039 CALL PTICR5(W,U,MN,R)
0040 102 CONTINUE
0041 RETURN
0042 END


0001 SUBROUTINE PTILRA(AQ,U,C,D,N,MN)
C * transformation of discrepancy functional
C * to the form (C*X,X)+(D,X)+E
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION AQ(MN,N),U(MN),C(N,N),D(N)
C * D:=-2A'*U
0005 DO 6 I=1,N
0006 R=0.0
0007 DO 5 J=1,MN
0008 R=R+AQ(J,I)*U(J)
0009 5 CONTINUE
0010 D(I)=-2.0*R
0011 6 CONTINUE
C * L:=A'*A
0012 DO 7 I=1,N
0013 DO 7 J=1,I
0014 R=0.0
0015 DO 8 K=1,MN
0016 R=R+AQ(K,I)*AQ(K,J)
0017 8 CONTINUE
0018 C(I,J)=R
0019 C(J,I)=R
0020 7 CONTINUE
0021 RETURN
0022 END


0001 SUBROUTINE PTILRB(A,B,N,NN,ITASK,C,IERR)
C * forming constaints matrix
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION A(NN,N),B(NN)
0005 EXTERNAL PTICR2
C * invalid restriction type set
0006 IERR=99
0007 IF(ITASK.GT.3.OR.ITASK.LT.1) GOTO 777
C * for the monotoneous functions wa are to
C * test restriction constant to be non-zero
0008 IERR=98
0009 IF(ITASK.EQ.1.AND.C.EQ.0.0) GOTO 777
0010 IERR=0
C * send zero to the matrix
0011 CALL PTICR2(A,0.0D0,N*NN)
0012 L=N-1
0013 GOTO (1,2,3),ITASK
0014 1 CONTINUE
C * monotoneous functions
0015 DO 11 I=1,N
0016 A(I,I)=1.0
0017 A(I+1,I)=-1.0
0018 B(I+1)=0.0
0019 11 CONTINUE
0020 B(1)=C
0021 GOTO 777
0022 2 CONTINUE
C * convex functions
0023 DO 12 I=2,L
0024 A(I,I)=-2.0
0025 A(I,I-1)=1.0
0026 A(I,I+1)=1.0
0027 B(I)=0.0
0028 12 CONTINUE
0029 A(1,1)=-1.0
0030 A(N,N)=-1.0
0031 B(1)=0.0
0032 B(N)=0.0
0033 GOTO 777
0034 3 CONTINUE
C * monotoneous convex functions
0035 DO 13 I=2,L
0036 A(I,I)=-2.0
0037 A(I,I-1)=1.0
0038 A(I,I+1)=1.0
0039 B(I)=0.0
0040 13 CONTINUE
0041 A(1,1)=-1.0
0042 A(1,2)=1.0
0043 A(N,N)=-1.0
0044 B(1)=0.0
0045 B(N)=0.0
0046 777 CONTINUE
0047 RETURN
0048 END


0001 SUBROUTINE PTILR0(A,N)
C * getting inverse matrix
C * of symmetric positive definite matrix A(N,N)
C * N+1-th column of the matrix A ( A(N,N+1) )
C * is used as working array
C * only subdiagonal elements may be set
C * A(I,J), ГДE I>=J
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION A(N,1)
0005 DO 1 I=1,N
0006 I1=I+1
0007 DO 1 J=I,N
0008 J1=J+1
0009 X=A(J,I)
0010 I2=I-1
0011 IF(I2.LT.1)GOTO 2
0012 DO 3 K1=1,I2
0013 K=I-K1
0014 X=X-A(K,J1)*A(K,I1)
0015 3 CONTINUE
0016 2 CONTINUE
0017 IF(J.NE.I) GOTO 4
0018 Y=1./DSQRT(X)
0019 A(I,I1)=Y
0020 GOTO 1
0021 4 A(I,J1)=X*Y
0022 1 CONTINUE
0023 DO 5 I=1,N
0024 I3=I+1
0025 IF(I3.GT.N)GOTO 5
0026 DO 6 J=I3,N
0027 Z=0.
0028 J1=J+1
0029 J2=J-1
0030 DO 7 K1=I,J2
0031 K=J-1-K1+I
0032 7 Z=Z-A(K,J1)*A(I,K+1)
0033 A(I,J1)=Z*A(J,J1)
0034 6 CONTINUE
0035 5 CONTINUE
0036 DO 8 I=1,N
0037 DO 8 J=I,N
0038 Z=0.
0039 J1=N+1
0040 J4=J+1
0041 DO 9 K=J4,J1
0042 9 Z=Z+A(J,K)*A(I,K)
0043 A(I,J+1)=Z
0044 8 CONTINUE
0045 DO 11 I=1,N
0046 DO 11 J=I,N
0047 A(J,I)=A(I,J+1)
0048 A(I,J)=A(I,J+1)
0049 11 CONTINUE
0050 RETURN
0051 END


0001 SUBROUTINE PTISR(AK,U0,X1,X2,Y1,Y2,N,M,
*Z,AN2,ITER,DL,IMAX,IC,
*R,NR,IERR)
C * Program to solve integral equations
C * of the first kind
C * by conjugate gradient projection method
C * dispatcher
C * work array mapping
C * name length what
C * A: N*M operator matrix
C * H: N descent direction
C * G: N gradient
C * U: M operator value
C * U1: M working array
C * S: N working array
C *
C * NR=N*M+3N+2M
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION U0(M),Z(N),R(NR)
0005 EXTERNAL AK
0006 EXTERNAL PTICR0,PTICR1,
*PTISR1,PTISR2,PTISR3
0007 ICONT=0
C * ICONT - start/continue mode flag
C * ICONT=0 start mode
C * ICONT=1 continue mode
0008 100 CONTINUE
C * getting arrays starts
0009 NA=1
0010 NH=N*M+1
0011 NG=NH+N
0012 NU=NG+N
0013 NU1=NU+M
0014 NS=NU1+M
0015 NMAX=NS+N
0016 IF(NMAX-1.GT.NR)GOTO 64
0017 IF(ICONT.EQ.1)GOTO 101
C * getting operator matrix
0018 CALL PTICR0(AK,R(NA),X1,X2,Y1,Y2,N,M)
C * tranformation to the П-plus
0019 CALL PTISR2(R(NA),Z,N,M,IC,R(NS))
0020 CALL PTICR1(Z,Z,R(NH),N,0D0)
0021 101 CONTINUE
C * minimizing discrepancy
0022 CALL PTISR1(R(NA),R(NH),U0,N,M,ITER,
*DL*(M-1.)/(Y2-Y1),0D0,IMAX,AN2,0D0,0D0,
*Z,R(NU),R(NU1),R(NH),R(NG),R(NS),IERR)
0023 AN2=AN2*(Y2-Y1)/(M-1.)
C * return from П-plus to original coordinates
0024 CALL PTICR1(Z,Z,R(NH),N,0D0)
0025 CALL PTISR3(Z,N,IC,R(NS))
0026 GOTO 999
0027 ENTRY PTISRE
C * entry to continue calculations
0028 ICONT=1
0029 64 CONTINUE
C * working array is too short
0030 IERR=64
0031 GOTO 999
0032 999 CONTINUE
0033 RETURN
0034 END

0001 SUBROUTINE PTISR1(A,Z0,U0,N,M,ITER,DL2,
*ANGRD,IMAX,AN2,ALF,RO,
*Z,U,U1,H,G,PLUS,IERR)
C * minimization of Thikhonov's functional
C * by means of conjugate gradient method
C * in the first quadrant
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION U0(M),Z0(N),Z(N),U(M),
*A(M,N),H(N),U1(M),G(N),PLUS(N)
0005 EXTERNAL PTICI2,PTICR3,PTICR5,PTICR4
0006 EXTERNAL PTICR6,PTICR7,PTICR1
0007 EXTERNAL PTICR9,PTICR8
0008 ITERBN=0
0009 IEND=0
0010 ITER=0
0011 JCH=-1
0012 CALL PTICR1(Z0,Z0,Z,N,0D0)
0013 CALL PTICR3(A,Z,U,N,M)
0014 CALL PTICR5(U,U0,M,AN2)
0015 CALL PTICR9(AS2O,AN2,Z,N,ALF,RO)
0016 CALL PTICR4(G,U,U0,A,N,M)
0017 CALL PTICR8(G,Z,N,ALF,RO)
0018 CALL PTICR2(H,0D0,N)
C * starting iterations
0019 14 CONTINUE
0020 ITER=ITER+1
0021 13 ICH=0
C * ICH - set changed flag
C * ICH=0 - active constraints was not changed
0022 IF(JCH) 1,2,3
0023 1 CALL PTICR2(PLUS,1D0,N)
C * PLUS(I)=0 - I-th constraint is active
C * PLUS(I)=1 - I-th constraint isn't active
0024 DO 4 I=1,N
0025 IF(Z(I).GT.0..OR.G(I).LT..0)GOTO 5
0026 PLUS(I)=0
0027 Z(I)=0.
0028 5 CONTINUE
0029 4 CONTINUE
C * active constraints set is ready
0030 IDIM=0
0031 DO 15 I=1,N
0032 15 IDIM=IDIM+PLUS(I)+0.5
C * IDIM - current face dimension
0033 ICH=1
0034 GOTO 2
0035 3 CONTINUE
C * adding new active constraint
0036 PLUS(IPRIM)=0
0037 Z(IPRIM)=0
0038 6 ICH=1
0039 IDIM=IDIM-1
0040 2 CONTINUE
0041 CALL PTICR7(G,G,PLUS,N,ANGRI)
0042 ITERBN=ITERBN+1
C * ITERBN - face iterations count
0043 IF(ICH.EQ.1)ITERBN=1
0044 IF(ANGRI.GT.ANGRD.AND.ITERBN.NE.
.IDIM+1.AND.IEND.NE.1)GOTO 7
0045 IF(IEND.EQ.1.AND.IDIM.EQ.IDIMO)GOTO 99
C * if IEND=1 - test if exact minimum is reached
0046 IEND=1-IEND
0047 IDIMO=IDIM
0048 IF(IEND.EQ.0)GOTO 7
0049 JCH=-1
0050 GOTO 8
0051 7 CONTINUE
C * getting descent direction - H
0052 BT=0.
0053 IF(ICH.EQ.0)BT=ANGRI/ANGRIO
C * saving gradient projection norm
0054 ANGRIO=ANGRI
0055 DO 9 I=1,N
0056 9 H(I)=(BT*H(I)+G(I))*PLUS(I)
C * descent direction found
C * ****** Attention? mashine constant
0057 ALMAX=1.E18
0058 IPRIM=0
0059 DO 10 I=1,N
0060 IF(H(I).LE..0)GOTO 11
0061 AL=Z(I)/H(I)
0062 IF(ALMAX.LT.AL)GOTO 12
0063 ALMAX=AL
0064 IPRIM=I
0065 12 CONTINUE
0066 11 CONTINUE
0067 10 CONTINUE
C * maximum step ALMAX found
C * IPRIM - number of the new possible active constraint
0068 CALL PTICRO(A,Z,G,U1,H,ALMAX,AL,N,M,ALF,RO,IED)
0069 CALL PTICR1(Z,H,Z,N,-AL)
0070 CALL PTICR3(A,Z,U,N,M)
0071 CALL PTICR4(G,U,U0,A,N,M)
0072 CALL PTICR8(G,Z,N,ALF,RO)
0073 CALL PTICR7(G,G,PLUS,N,ANGRI)
0074 CALL PTICR5(U,U0,M,AN2)
0075 CALL PTICR9(AS2,AN2,Z,N,ALF,RO)
0076 JCH=0
0077 IF(IED.EQ.1)JCH=1
C * IED=1 - getting new constraint
0078 IF(IED.EQ.2)GOTO 22
0079 IF(AN2.LE.DL2.OR.ITER.GE.IMAX)GOTO 20
0080 IF(AS2.GE.AS2O)GOTO 21
0081 AS2O=AS2
0082 GOTO 14
0083 20 CONTINUE
C * exit by discrepancy level or iteration number
0084 IERR=1
0085 GOTO 999
0086 99 CONTINUE
C * the exact minimum is found
0087 IERR=0
0088 GOTO 999
0089 8 CONTINUE
C * return to the beginning to check exact minimum
0090 GOTO 13
0091 21 CONTINUE
C * functional isn't changed
0092 IERR=2
0093 GOTO 999
0094 22 CONTINUE
C * negative step found
0095 IERR=65
0096 GOTO 999
0097 999 CONTINUE
0098 RETURN
0099 END


0001 SUBROUTINE PTISR2(A,Z,N,M,IC,S)
C * transformation of the operator matrix A(M,N)
C * and initial approximation Z
C * to the first quadrant
C * S(N) - working array
C * IC - correcteness set
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION A(M,N),Z(N),S(N)
0005 EXTERNAL PTICR1,PTISR4
0006 IC1=IABS(IC)+1
0007 N1=N-1
0008 GOTO (800,801,802,803,804,805),IC1
0009 800 CONTINUE
C * non-negative functions
C * nothing to do
0010 GOTO 999
C * transformation of the initial point - Z
0011 801 CONTINUE
C * monotoneous functions, IC=1
0012 DO 101 I=2,N
0013 101 S(I-1)=Z(I-1)-Z(I)
0014 S(N)=Z(N)
0015 GOTO 799
0016 802 CONTINUE
C * descending, convex-up functions IC=2
0017 DO 102 I=2,N1
0018 102 S(I)=(I-N)*(Z(I-1)-2.*Z(I)+Z(I+1))
0019 S(N)=Z(N)
0020 S(1)=(Z(1)-Z(2))*(N-1.)
0021 GOTO 799
0022 803 CONTINUE
C * convex-up functions IC=3
0023 DO 103 I=2,N1
0024 103 S(I)=(Z(I-1)-2.*Z(I)+
+Z(I+1))*(N-I)*(I-1.)/(1.-N)
0025 S(N)=Z(N)
0026 S(1)=Z(1)
0027 GOTO 799
0028 804 CONTINUE
C * descending, convex-down functions IC=4
0029 DO 104 I=2,N1
0030 104 S(I-1)=(Z(I-1)-2.*Z(I)+Z(I+1))*(I-0.)
0031 S(N1)=(Z(N1)-Z(N))*(N-1.)
0032 S(N)=Z(N)
0033 GOTO 799
0034 805 CONTINUE
C * convex-down, positive at ends functions IC=5
0035 DO 105 I=2,N1
0036 105 S(I)=(Z(I-1)-2.*Z(I)+
+Z(I+1))*(N-I)*(I-1.)/(N-1.)
0037 S(N)=Z(N)
0038 S(1)=Z(1)
0039 GOTO 799
0040 799 CONTINUE
0041 CALL PTICR1(S,S,Z,N,0D0)
C * transformation of matrix A
0042 DO 1 K=1,M
C * filling array S
0043 CALL PTICR2(S,0D0,N)
0044 DO 2 J=1,N
0045 DO 3 I=1,N
C * calculating the element of transformation matrix
0046 T=PTISR4(I,J,IC,N)
0047 3 S(J)=S(J)+A(K,I)*T
0048 2 CONTINUE
C * K-th row of matrix A is converted.
0049 DO 4 J=1,N
0050 4 A(K,J)=S(J)
0051 1 CONTINUE
0052 999 CONTINUE
0053 RETURN
0054 END


0001 FUNCTION PTISR4(I,J,IC,N)
C * calculation of transformation matrix elements
C * from П-plus to compact, specified by parameter IC
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 IC1=IABS(IC)+1
0005 GOTO (800,801,802,803,804,805),IC1
0006 800 CONTINUE
C * non-negative
0007 H=0.
0008 IF(I.EQ.J)H=1.
0009 GOTO 998
0010 801 CONTINUE
C * descending functions - IC=1
0011 H=0.
0012 IF(I.LE.J)H=1.
0013 GOTO 998
0014 802 CONTINUE
C * descending, convex-up functions IC=2
0015 H=1.
0016 IF(I.GT.J)H=(N-I)/(N-J+0.)
0017 GOTO 998
0018 803 CONTINUE
C * convex-up IC=3
0019 IF(I.LE.J.AND.J.NE.1)H=(I-1.)/(J-1.)
0020 IF(I.GE.J.AND.J.NE.N)H=(N-I)/(N-J+0.)
0021 GOTO 998
0022 804 CONTINUE
C * descending, convex-down functions IC=4
0023 H=0.
0024 IF(I.LE.J.AND.J.NE.N)H=(J-I+1.)/(J+0.)
0025 IF(J.EQ.N)H=1.
0026 GOTO 998
0027 805 CONTINUE
C * convex-down, positive at ends functions IC=5
0028 IF(I.LE.J.AND.J.NE.1)H=(I-1.)/(J-1.)
0029 IF(I.GE.J.AND.J.NE.N)H=(N-I)/(N-J+0.)
0030 IF(J.NE.1.AND.J.NE.N)H=-H
0031 GOTO 998
0032 998 CONTINUE
0033 PTISR4=H
0034 RETURN
0035 END


0001 SUBROUTINE PTISR3(Z,N,IC,Y)
C * transformation from П-plus to the set specified by IC
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION Z(N),Y(N)
0005 EXTERNAL PTISR4,PTICR1
0006 DO 1 I=1,N
0007 S=0.
0008 DO 2 J=1,N
0009 2 S=S+Z(J)*PTISR4(I,J,IC,N)
0010 Y(I)=S
0011 1 CONTINUE
0012 CALL PTICR1(Y,Y,Z,N,0D0)
0013 RETURN
0014 END


0001 SUBROUTINE PTICR0(AK,A,X1,X2,Y1,Y2,N,M)
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
C * Gets system of linear equations
C * for finite-difference approximation of
C * integral equation of the first type
C * AK(X,Y) - integral equation kernel function
C * X1,X2 - integration limits
C * Y1,Y2 - Y-variable range
C * N - number of numerical integration points
C * M - number of points of Y-grid
C * A(M,N) - linear system matrix (output)
C * Integration by the second variable
C * of the function AK(X,Y)
0004 DIMENSION A(M,N)
0005 H=(X2-X1)/(N-1.)
0006 DO 1 I=1,N
0007 X=X1+(X2-X1)/(N-1.)*(I-1.)
0008 DO 2 J=1,M
0009 Y=Y1+(Y2-Y1)/(M-1.)*(J-1.)
0010 A(J,I)=AK(X,Y)*H
0011 IF(I.EQ.1.OR.I.EQ.N)A(J,I)=A(J,I)/2.
0012 2 CONTINUE
0013 1 CONTINUE
0014 RETURN
0015 END


0001 SUBROUTINE PTICRO(A,Z,G,U,H,ALM,AL,N,M,
*ALF,RO,IED)
C * one-dimensional minimization
C * of the smoosing functional
C * on half-line Z-AL*H AL>0
C * A(M,N) - linear operator matrix
C * Z(N) - start point
C * G(N) - gradient at the point Z
C * H(N) - descent direction
C * ALM - step restriction AL C * U(M) - working array
C * AL - step found
C * ALF - regularization parameter
C * RO - wight of difference derivative in norm W21
C * IED - return code
C * 0 - normal - minimum found,
C * 1 - minimum on the restriction found,
C * 2 - negative step found
C * 3 - zero divide
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION A(M,N),Z(N),G(N),U(M),H(N)
0005 EXTERNAL PTICR6,PTICR3
0006 H2=0.
0007 CALL PTICR6(G,H,N,GH)
0008 IF(GH.LT.0.)GOTO 2
0009 CALL PTICR3(A,H,U,N,M)
0010 CALL PTICR6(U,U,M,AH2)
0011 IF(ALF.EQ.0.)GOTO 3
C * regularization paremeter doesn't equal zero
C * find norm of the descent direction
0012 DO 4 I=2,N
0013 4 H2=H2+(H(I)-H(I-1))**2*RO+H(I)**2
0014 H2=H2+H(1)**2
0015 3 CONTINUE
0016 IF(AH2+H2*ALF.EQ.0.0)GOTO 998
0017 AL=GH/(AH2+H2*ALF)*.5
0018 IF(AL.GE.ALM)GOTO 1
0019 IED=0
0020 GOTO 999
0021 1 CONTINUE
C * minimum on the restriction (bound) found
0022 AL=ALM
0023 IED=1
0024 GOTO 999
0025 2 CONTINUE
C * negative step found, make it be zero
0026 AL=0.
0027 IED=2
0028 GOTO 999
0029 998 CONTINUE
C * zero divide, let step be zero
0030 IED=3
0031 AL=0.
0032 GOTO 999
0033 999 CONTINUE
0034 RETURN
0035 END


0001 SUBROUTINE PTICR1(B,C,A,N,P)
C * array moving
C * A(I):=B(I)+P*C(I), I=1,N
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION A(N),B(N),C(N)
0005 DO 1 I=1,N
0006 1 A(I)=B(I)+P*C(I)
0007 RETURN
0008 END


0001 SUBROUTINE PTICR2(A,R,N)
C * fill real array A of length N
C * by the value R
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION A(N)
0005 DO 1 I=1,N
0006 1 A(I)=R
0007 RETURN
0008 END


0001 SUBROUTINE PTICI2(A,R,N)
C * fill integer array A of length N
C * by the value R
0002 INTEGER*4 N,I,R,A
0003 DIMENSION A(N)
0004 DO 1 I=1,N
0005 1 A(I)=R
0006 RETURN
0007 END


0001 SUBROUTINE PTICR3(A,Z,U,N,M)
C * multiply matrix A(M,N) by vector Z(N)
C * U(M) - vector-result
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION U(M),A(M,N),Z(N)
0005 DO 1 I=1,M
0006 S=0.
0007 DO 2 J=1,N
0008 2 S=S+A(I,J)*Z(J)
0009 1 U(I)=S
0010 RETURN
0011 END


0001 SUBROUTINE PTICR4(G,U,U0,A,N,M)
C * calculating gradient of discrepancy (norm AZ-U0)
C * U=A*Z value of operator at the point Z
C * G(N) - result
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION G(N),U(M),U0(M),A(M,N)
0005 DO 1 J=1,N
0006 S=0.
0007 DO 2 I=1,M
0008 2 S=S+A(I,J)*(U(I)-U0(I))
0009 1 G(J)=S*2.
0010 RETURN
0011 END


0001 SUBROUTINE PTICR5(A,B,N,S)
C * descrepancy calculation
C * S - result
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION A(N),B(N)
0005 S=0.
0006 DO 1 I=1,N
0007 1 S=S+(A(I)-B(I))**2
0008 RETURN
0009 END


0001 SUBROUTINE PTICR6(A,B,N,S)
C * scalar product calculation of vectors A and B
C * S - result
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION A(N),B(N)
0005 S=0.
0006 DO 1 I=1,N
0007 1 S=S+A(I)*B(I)
0008 RETURN
0009 END


0001 SUBROUTINE PTICR7(A,B,C,N,S)
C * scalar product of vectors A and B of length N
C * with wights C(I)
C * S - result
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION A(N),B(N),C(N)
0005 S=0.
0006 DO 1 I=1,N
0007 1 S=S+A(I)*B(I)*C(I)
0008 RETURN
0009 END


0001 SUBROUTINE PTICR8(G,Z,N,ALF,RO)
C * calculation of stabilizer gradient
C * with wight RO
C * and adding it to vector G(N)
C * Z(N) - point to calculate
C * ALF - regularization parameter
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION Z(N),G(N)
0005 IF(ALF.EQ.0.)GOTO 999
0006 DO 1 I=1,N
0007 G(I)=G(I)+2.*ALF*Z(I)
0008 IF(I.NE.1.AND.I.NE.N)G(I)=G(I)+
*2.*ALF*(2.*Z(I)-Z(I-1)-Z(I+1))*RO
0009 IF(I.EQ.1)
+ G(I)=G(I)+2.*ALF*(Z(1)-Z(2))*RO
0010 IF(I.EQ.N)
+ G(I)=G(I)+2.*ALF*(Z(N)-Z(N-1))*RO
0011 1 CONTINUE
0012 999 CONTINUE
0013 RETURN
0014 END


0001 SUBROUTINE PTICR9(AS2,AN2,Z,N,ALF,RO)
C * calculate stabilizer and add it to discrepancy AN2
C * Z(N) - point
C * ALF - regularization parameter
C * RO - wight of difference derivative in norm W21
C * AS2 - result
0002 IMPLICIT REAL*8(A-H,O-Z)
0003 IMPLICIT INTEGER*4(I-N)
0004 DIMENSION Z(N)
0005 S=0.
0006 IF(ALF.EQ.0.)GOTO 999
0007 DO 4 I=2,N
0008 4 S=S+(Z(I)-Z(I-1))**2*RO+Z(I)**2
0009 S=S+Z(1)**2
0010 999 CONTINUE
0011 AS2=AN2+S*ALF
0012 RETURN
0013 END