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

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