CTITLEMCGRDX -- CONJUG. GRADIENT MINIMIZ. (3D STATICS + RNMO/GATE) 00010000 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR D.D. THOMPSON R&D 00020000 CA DESIGNER D.D. THOMPSON R&D 00030000 CA LANGUAGE VS FORTRAN 00040000 CA SYSTEM IBM 00050000 CA WRITTEN 2/15/77 00060000 C REVISED 12/15/77 MENDEKE - DELETED PRINT 00070000 C REVISED 12/21/77 D.D.REED - ADDED CODE TO INTERCEPT UNDERFLOWS IN00080000 C ITERATION LOOP. 00090000 C REVISED 01/18/78 D.D.THOMPSON - MODIFIED TO HANDLE 3-D AND 00100000 C SEPARATE RNMO/GATE. 00110000 C NAME CHANGED FROM MCGRDM. 00120000 C REVISED 11/21/80 POLAK - MOVED CHECK OF CZ <1.0D-10 00130000 C REVISED 02/04/81 POLAK - CORRECTED CALCULATION OF INDEX INTO 00140000 C XRNMO ARRAY. 00150000 C REVISED 06/18/81 POLAK - ADDED CALL TO MDQONX IN LOOP 20 BEFORE 00160000 C THE MAIN ITERATION LOOP. 00170000 C REVISED 11-08-82 COOPER - RECOMPILE PROGRAM TO MATCH LOAD. 00180000 C REVISED 01-15-84 REP - ADD CALLS TO USBFRX FOR BUFFERING TO DISK.00190000 C REVISED 08-14-85 REP. CHANGE NAME FROM MCGRD3. 00200000 C REVISED 11-09-88 ESN. INCORPORATE MEMORY PATH. 00210000 C REVISED 12-19-89 ESN. ADD IN CHECKPOINT PRINTS. 00220002 CA 00230000 CA CALL MCGRDX (IND, ICDP, KCDP, KCWIND, DSQ, WR, WS, X, XCDP, 00240000 CA XRNMO, Q, P, R, U, ZA, ZB, ZC, ZD, ZE, LD, LR, LS, 00250000 CA LWIND, NWIND, LINE, NLINE, NUMC, NUML, IPR) 00260000 CA 00270000 CA IN/OUT ARGUMENT TYPE DESCRIPTION 00280000 CA 00290000 CA IN IND I4 BIASED STACKING TABLES. MIN DIMENSION (2,NT) 00300000 CA WHERE NT = LWIND(NWIND) 00310000 CA IND(1,I) = BIASED RECEIVER INDEX FOR TRACE I 00320000 CA IND(2,I) = BIASED SHOT INDEX FOR TRACE I 00330000 CA BIASED ARE DETERMINED SO THAT THE LOWEST 00340000 CA RECEIVER INDEX TO APPEAR IS > OR = 1 AND 00350000 CA LOWEST SHOT INDEX TO APPEAR IS > OR = LR+1. 00360000 CA IN ICDP I4 ARRAY OF INDICES CORRESPONDING TO THE LAST 00370000 CA TRACE IN EACH CDP. LENGTH <= LD*NLINE*NWIND. 00380000 CA NOTE AS WITH LINE INDICES ARE RELATIVE TO THE00390000 CA START OF EACH WINDOW. 00400000 CA IN KCDP I4 ARRAY OF STARTING CDP NUMBERS FOR EACH LINE. 00410000 CA (SMALLEST MUST BE 1) LENGTH = NLINE*NWIND 00420000 CA IN KCWIND I4 ARRAY OF STARTING ADDRESSES IN ICDP CORRES- 00430000 CA PONDING TO THE BEGINNING OF EACH WINDOW. 00440000 CA (KCWIND(1) = 1) LENGTH = NWIND 00450000 CA IN DSQ R4 ARRAY OF OFFSET DISTANCES SQUARED. LENGTH =NT00460000 CA (IT IS RECOMMENDED THAT THE DISTANCES BE 00470000 CA NORMALIZED TO MAX. OF 1 BEFORE SQUARING.) 00480000 CA IN WR R4 RECEIVER AUXILLIARY WEIGHT. 00490000 CA IN WS R4 SHOT AUXILLIARY WEIGHT. 00500000 CA OUT X R8 RECEIVER AND SHOT STATICS. 00510000 CA X(1)--X(LR) = RECEIVER STATIC FACTORS 00520000 CA X(LR+1)--X(LS)= SHOT STATIC FACTORS 00530000 CA OUT XCDP R8 NEGATIVE OF THE PILOT STATIC. 00540000 CA ELEMENTS: 00550000 CA 1 -- LD FOR LINE=1, WINDOW=1 00560000 CA LD+1 -- 2*LD FOR LINE=2, WINDOW=1 00570000 CA ... 00580000 CA (NLINE-1)*LD+1 -- NLINE*LD FOR LINE=NLINE, 00590000 CA WINDOW=1 00600000 CA OUT XRNMO R8 RNMO FOR EACH CDP AT OFFSET = 1. 00610000 CA (ARRANGED SAME AS ZCDP) (IF DSQ IS NORMALIZED00620000 CA TO 1, THEN THESE TERMS WILL BE THE ACTUAL 00630000 CA RNMO OF A TRACE HAVING MAX. OFFSET DISTANCE 00640000 CA IF IT WERE TO APPEAR IN EACH CDP.) 00650000 CA IN Q R8 ARRAY OF MEASURED LAGS FOR EACH TRACE. LAGS 00660000 CA NEED ONLY BE RELATIVE WITHIN EACH GATHER. 00670000 CA REJECTED LAGS INDICATED BY Q > 1.0E6. Q IS 00680000 CA NOT PRESERVED UPON RETURN. 00690000 CA IN P R8 WORK ARRAY OF MIN. LENGTH LS 00700000 CA IN R R8 WORK ARRAY OF MIN. LENGTH LS 00710000 CA IN U R8 WORK ARRAY OF MIN. LENGTH LS 00720000 CA IN ZA R8 WORK ARRAY OF MIN. LENGTH NLINE * LD 00730000 CA IN ZB R8 WORK ARRAY OF MIN. LENGTH NLINE *LD 00740000 CA IN ZC R8 WORK ARRAY OF MIN. LENGTH NLINE*LD 00750000 CA IN ZD R8 WORK ARRAY OF MIN. LENGTH NLINE*LD 00760000 CA IN ZE R8 WORK ARRAY OF MIN. LENGTH NLINE* LD 00770000 CA IN LD I4 > OR = THE LAST BIASED CDP INDEX TO APPEAR. 00780000 CA IN LR I4 > OR = THE LAST BIASED REC. INDEX TO APPEAR. 00790000 CA IN LS I4 > OR = THE LAST BIASED SHOT INDEX TO APPEAR. 00800000 CA IN LWIND I4 LAST INDEX IN Y FOR EACH WINDOW. LENGTH = 00810000 CA NWIND 00820000 CA IN NWIND I4 NUMBER OF WINDOWS. 00830000 CA IN LINE I4 LAST ELEMENT FOR EACH LINE. ENTRIES FOR EACH00840000 CA WINDOW ARE RELATIVE TO START OF WINDOW. 00850000 CA LENGTH = NLINE*NWIND 00860000 CA IN NLINE I4 NUMBER OF LINES. 00870000 CA IN NUMC I4 # OF CDP IN EACH SIDE OF CENTER TO AVERAGE 00880000 CA RNMO. (NUMC = 0 IS OK) 00890000 CA IN NUML I4 # OF LINES ON EACH SIDE OF CENTER TO AVERAGE 00900000 CA RNMO. (NUML = 0 IS OK) 00910000 CA IN IPR I4 UNIT NUMBER TO PRINT 00920000 CA 00930000 CA 00940000 CA 00950000 CA THIS ROUTINE UNSCRAMBLES THE STATIC COMPONENTS FOR CDP DATA. 00960000 CA 00970000 SUBROUTINE MCGRDX (IND, ICDP, KCDP, KCWIND, DSQ, WR, WS, X, XCDP, 00980000 * XRNMO, Q, P, R, U, ZA, ZB, ZC, ZD, ZE, LD, LR, 00990000 * LS, LWIND, NWIND, LINE, NLINE, NUMC, NUML, IPR)01000000 C 01010000 IMPLICIT REAL*8 (A-H,O-Z) 01020000 EXTERNAL SAUFER 01030000 C 01040000 COMMON COM (1) 01050000 INTEGER COM 01060000 REAL XCOM (1) 01070000 REAL*8 ZCOM (1) 01080000 EQUIVALENCE (COM(1),XCOM(1),ZCOM(1)) 01090000 C 01100000 C COMMON FOR BUFFERING INFORMATION 01110000 C 01120000 COMMON /BFINFO/ BYND(15), BIND(15), BOFF(15), BQND(15) 01130000 C 01140000 INTEGER BYND, BIND, BOFF, BQND, B1, B2 01150000 C 01160000 C 01170000 C REAL ARRAYS IN PARAMETER LIST 01180000 C 01190000 DIMENSION P (1) 01200000 DIMENSION Q (1) 01210000 DIMENSION R (1) 01220000 DIMENSION U (1) 01230000 DIMENSION X (1) 01240000 DIMENSION XCDP (1) 01250000 DIMENSION XRNMO (1) 01260000 DIMENSION ZA (1) 01270000 DIMENSION ZB (1) 01280000 DIMENSION ZC (1) 01290000 DIMENSION ZD (1) 01300000 DIMENSION ZE (1) 01310000 REAL*4 DSQ (1) 01320000 C 01330000 C INTEGER ARRAYS IN PARAMETER LIST 01340000 C 01350000 INTEGER ICDP (1) 01360000 INTEGER IND (2,1) 01370000 INTEGER KCDP (1) 01380000 INTEGER KCWIND (1) 01390000 INTEGER LINE (1) 01400000 INTEGER LWIND (1) 01410000 C 01420000 C REAL VARIABLES IN PARAMETER LIST 01430000 C 01440000 REAL WR 01450000 REAL WS 01460000 C 01470000 C REAL VARIABLES -- LOCAL 01480000 C 01490000 REAL*16 D 01500000 REAL*16 S 01510000 REAL*16 Z 01520000 C 01530000 LOGICAL*1 ERRORF 01540000 COMMON /UFERR/ ERRORF 01550000 C 01560000 ERRORF = .FALSE. 01570000 C 01580000 NT = LWIND(NWIND) 01590000 WR2 = WR * WR 01600000 WS2 = WS * WS 01610000 C 01620000 WRITE (IPR, 9020) LS 01630000 C 01640000 LDP = 1 01650000 C 01660000 C CLEAR OUTPUT (X) AND LOAD STATIC VALUES INTO SCRATCH 01670000 C ARRAY (Q). 01680000 C 01690000 DO 10 I = LDP, LS 01700000 10 X(I) = 0.0 01710000 C 01720000 LA = 1 01730000 C 01740000 DO 20 I = 1, NWIND 01750000 JJ = NLINE * (I-1) + 1 01760000 J = (I-1) * LD * NLINE + 1 01770000 CALL MNOQDX (Q, LA-1, ICDP(KCWIND(I)), KCDP(JJ), LINE(JJ), LD, 01780000 * NLINE, DSQ, XCDP(J), XRNMO(J), ZA, ZB, ZC, NUMC, NUML) 01790000 C 01800000 CALL MDQONX (Q, LA-1, ICDP(KCWIND(I)), KCDP(JJ), LINE(JJ), LD, 01810000 * NLINE, DSQ, ZA, ZB, ZC, ZD, ZE, NUMC, NUML) 01820000 20 LA = LWIND(I) + 1 01830000 C 01840000 C CALL TO MTRANS REPLACED BY IN-LINE CODE. 01850000 C 01860000 DO 30 I = LDP, LS 01870000 30 P(I) = 0.0 01880000 C 01890000 IF (BQND(12) .EQ. 0) THEN 01900000 DO 40 I = 1, NT 01910000 CALL USBFRX (Q, I, B1, 0, BQND) 01920000 T = Q(B1) 01930000 IF(T .GT. 1.0E6) GO TO 40 01940000 C 01950000 C CALL TO MSCINM REPLACED BY IN-LINE CODE. 01960000 C 01970000 CALL USBFRX (IND, I, B2, 0, BIND) 01980000 KR = IND(1,B2) 01990000 KS = IND(2,B2) 02000000 P(KR) = P(KR) + T 02010000 P(KS) = P(KS) + T 02020000 40 CONTINUE 02030000 ELSE 02040000 B1 = (BQND(12)+1) / 2 02050000 DO 45 I = 1, NT 02060000 T = ZCOM(B1+I-1) 02070000 IF(T .GT. 1.0E6) GO TO 45 02080000 C 02090000 C CALL TO MSCINM REPLACED BY IN-LINE CODE. 02100000 C 02110000 KR = COM(BIND(12)+2*I-2) 02120000 KS = COM(BIND(12)+2*I-1) 02130000 P(KR) = P(KR) + T 02140000 P(KS) = P(KS) + T 02150000 45 CONTINUE 02160000 ENDIF 02170000 C 02180000 C COPY INITIAL CONJUGATE DIRECTION (P) INTO RESIDUAL VECTOR (R). 02190000 C 02200000 DO 50 I2 = LDP, LS 02210000 50 R(I2) = P(I2) 02220000 C 02230000 C SET # ITERATIONS (NPARP) AT MORE THAN MIN. FOR ACCURACY 02240000 C IMPROVEMENT. 02250000 C 02260000 NPARP = LS + LS / 10 02270000 NCHECK = NPARP / 20 02280000 WRITE (IPR,9100) NPARP,NCHECK 02290000 C 02300000 C ******** BEGIN MAIN CONJUGATE GRADIENT ITERATION LOOP **********02310000 C 02320000 C 02330000 CALL ERRSAV (208, TAB208) 02340000 CALL ERRSET (208, 2, -1, 1, SAUFER) 02350000 C 02360000 DO 160 I3 = 1, NPARP 02370000 I3X = I3 02380000 S = 0.0 02390000 C 02400000 DO 60 I4 = LDP, LS 02410000 60 S = S + P(I4) * R(I4) 02420000 C 02430000 C CS = P'R 02440000 C 02450000 CS = S 02460000 C 02470000 C U = X'B'TBX COMPUTED IN 5 SUBROUTINE CALLS BELOW. 02480000 C 02490000 C CALL TO MODELM REPLACED BY IN-LINE CODE. 02500000 C 02510000 IF (BQND(12) .EQ. 0) THEN 02520000 DO 70 I = 1, NT 02530000 CALL USBFRX (Q, I, B1, 1, BQND) 02540000 IF(Q(B1) .GT. 1.0E6) GO TO 70 02550000 CALL USBFRX (IND, I, B2, 0, BIND) 02560000 Q(B1) = P(IND(1,B2)) + P(IND(2,B2)) 02570000 70 CONTINUE 02580000 ELSE 02590000 B1 = (BQND(12)+1) / 2 02600000 DO 75 I = 1, NT 02610000 IF(ZCOM(B1+I-1) .GT. 1.0E6) GO TO 75 02620000 ZCOM(B1+I-1) = P(COM(BIND(12)+2*I-2)) + 02630000 * P(COM(BIND(12)+2*I-1)) 02640000 75 CONTINUE 02650000 ENDIF 02660000 C 02670000 LA = 1 02680000 C 02690000 DO 80 I = 1, NWIND 02700000 JJ = NLINE * (I-1) + 1 02710000 CALL MNOQDX (Q, LA-1, ICDP(KCWIND(I)), KCDP(JJ), LINE(JJ), LD,02720000 * NLINE, DSQ, ZA, ZB, ZC, ZD, ZE, NUMC, NUML) 02730000 C 02740000 CALL MDQONX (Q, LA-1, ICDP(KCWIND(I)), KCDP(JJ), LINE(JJ), LD,02750000 * NLINE, DSQ, ZA, ZB, ZC, ZD, ZE, NUMC, NUML) 02760000 80 LA = LWIND(I) + 1 02770000 C 02780000 C CALL TO MTRANS REPLACED BY IN-LINE CODE. 02790000 C 02800000 KA = LDP 02810000 C 02820000 DO 90 I = KA, LR 02830000 90 U(I) = WR2 * P(I) 02840000 C 02850000 KA = LR + 1 02860000 C 02870000 DO 100 I = KA, LS 02880000 100 U(I) = WS2 * P(I) 02890000 C 02900000 C MSCINM REPLACED BY IN-LINE CODE. 02910000 C 02920000 IF (BQND(12) .EQ. 0) THEN 02930000 DO 110 I = 1, NT 02940000 CALL USBFRX (IND, I, B2, 0, BIND) 02950000 CALL USBFRX (Q , I, B1, 0, BQND) 02960000 KR = IND(1,B2) 02970000 KS = IND(2,B2) 02980000 T = Q(B1) 02990000 IF(T .GT. 1.0E6) GO TO 110 03000000 U(KR) = U(KR) + T 03010000 U(KS) = U(KS) + T 03020000 110 CONTINUE 03030000 ELSE 03040000 B1 = (BQND(12)+1) / 2 03050000 DO 115 I = 1, NT 03060000 KR = COM(BIND(12)+2*I-2) 03070000 KS = COM(BIND(12)+2*I-1) 03080000 T = ZCOM(B1+I-1) 03090000 IF(T .GT. 1.0E6) GO TO 115 03100000 U(KR) = U(KR) + T 03110000 U(KS) = U(KS) + T 03120000 115 CONTINUE 03130000 ENDIF 03140000 C 03150000 C ...END OF COMPUTATION OF U. 03160000 C 03170000 Z = 0.0 03180000 C 03190000 DO 120 I49 = LDP, LS 03200000 120 Z = Z + P(I49) * U(I49) 03210000 C 03220000 C CZ = P'U, SOZ = STEP SIZE TO TAKE IN CONJUGATE DIRECTION(P). 03230000 C (MODIFY BOTH (X) AND RESIDUAL(R) FOR THIS STEP.) 03240000 C 03250000 CZ = Z 03260000 IF(DABS(CZ) .LT. 1.0D-10) GO TO 200 03270000 SOZ = CS / CZ 03280000 XINC = 0.0 03290000 C 03300000 DO 130 I6 = LDP, LS 03310000 X(I6) = X(I6) + SOZ * P(I6) 03320000 XINC = DABS(SOZ * P(I6)) + XINC 03330000 130 CONTINUE 03340000 C 03350000 AVINC = XINC / LS 03360000 IF (AVINC .LT. 1.0D-10) GO TO 200 03370000 D = 0.0 03380000 C 03390000 DO 140 I7 = LDP, LS 03400000 R(I7) = R(I7) - SOZ * U(I7) 03410000 D = D + R(I7) * U(I7) 03420000 IF (ERRORF) GO TO 210 03430000 140 CONTINUE 03440000 C 03450000 C CD = U'(NEW R), DOZ = AMOUNT OF OLD CONJUGATE DIRECTION (P) 03460000 C TO ADD TO NEW RESIDUAL TO PRODUCE NEW CONJUGATE DIRECTION. 03470000 C 03480000 CD = D 03490000 IF(DABS(CD) .LT. 1.0D-10) GO TO 200 03500000 IF(DABS(CS) .LT. 1.0D-10) GO TO 200 03510000 DOZ = -CD / CZ 03520000 C 03530000 DO 150 I8 = LDP, LS 03540000 150 P(I8) = R(I8) + DOZ * P(I8) 03550000 C 03560000 IF ((I3X/NCHECK)*NCHECK .EQ. I3X) THEN 03570001 WRITE (IPR,9110) I3X,CZ,AVINC,CD,CS 03580000 ENDIF 03590000 C 03600000 160 CONTINUE 03610000 C 03620000 C **END MAIN CONJUGATE GRADIENT ITERATION LOOP** 03630000 C 03640000 200 WRITE (IPR, 9000) I3X 03650000 GO TO 220 03660000 C 03670000 210 WRITE (IPR, 9010) I3X 03680000 C 03690000 220 CALL ERRSTR (208, TAB208) 03700000 ERRORF = .FALSE. 03710000 C 03720000 C CALL TO MODELM REPLACED BY IN-LINE CODE. 03730000 C 03740000 IF (BQND(12) .EQ. 0) THEN 03750000 DO 230 I = 1, NT 03760000 CALL USBFRX (IND, I, B2, 0, BIND) 03770000 CALL USBFRX (Q , I, B1, 1, BQND) 03780000 IF(Q(B1) .GT. 1.0E6) GO TO 230 03790000 Q(B1) = X(IND(1,B2)) + X(IND(2,B2)) 03800000 230 CONTINUE 03810000 ELSE 03820000 B1 = (BQND(12)+1) / 2 03830000 DO 235 I = 1, NT 03840000 IF(ZCOM(B1+I-1) .GT. 1.0E6) GO TO 235 03850000 ZCOM(B1+I-1) = X(COM(BIND(12)+2*I-2)) + 03860000 * X(COM(BIND(12)+2*I-1)) 03870000 235 CONTINUE 03880000 ENDIF 03890000 C 03900000 LA = 1 03910000 C 03920000 DO 250 II = 1, NWIND 03930000 JJ = NLINE * (II - 1) + 1 03940000 CALL MNOQDX (Q, LA-1, ICDP(KCWIND(II)), KCDP(JJ), LINE(JJ), LD, 03950000 * NLINE, DSQ, ZA, ZB, ZC, ZD, ZE, NUMC, NUML) 03960000 J = (II-1) * LD * NLINE 03970000 C 03980000 DO 240 K = 1, NLINE 03990000 KK = (K-1) * LD 04000000 C 04010000 DO 240 I = 1, LD 04020000 I1 = I + KK 04030000 I2 = I1 + J 04040000 XCDP(I2) = XCDP(I2) - ZA(I1) 04050000 240 XRNMO(I2) = XRNMO(I2) - ZB(I1) 04060000 C 04070000 250 LA = LWIND(II) + 1 04080000 C 04090000 9000 FORMAT (/// ' STOP ON SMALL DELTA X AT ITERATION ', I6) 04100000 C 04110000 9010 FORMAT (/// ' STOP ON UF AT ITERATION ', I6) 04120000 C 04130000 9020 FORMAT (// ' NUMBER OF SHOTS + RECEIVERS = ', I8) 04140000 C 04150000 9100 FORMAT (/// ' TOTAL NUMBER OF ITERATIONS WILL BE ', I6,/, 04160000 * ' CHECKPOINTS WILL BE EVERY ', I6,' ITERATIONS '///) 04170000 C 04180000 9110 FORMAT (' AT ITERATION ',I5,' ERRORS ARE ',4G10.4) 04190000 C 04200000 RETURN 04210000 C 04220000 END 04230000