CTITLEMPRDCT -- COMPUTES DDLOG SYSTEM MATRIX 00000010 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR D. D. THOMPSON 00000020 CA DESIGNER D. D. THOMPSON 00000030 CA LANGUAGE FORTRAN 77 00000040 CA SYSTEM IBM & CRAY 00000041 CA WRITTEN 1972 00000050 C REVISED MO-DA-YR 00000060 C REVISED 05-20-85 TWH. ADAPTED TO IBM & CRAY. 00000061 C 00000070 CA 00000080 CA 00000090 CA CALL MPRDCT (Z, KR, KW, KX, V, PR, KEY) 00000100 CA INPUT Z = REFLECTION COEFFICIENT (RC'S) AND R4 00000110 CA WAVELET ARRAY 00000120 CA INPUT KR = NUMBER OF RC'S IN THE TOP OF Z I4 00000130 CA INPUT KW = NUMBER OF WAVELET SAMPLES IN I4 00000140 CA BOTTOM OF Z 00000150 CA INPUT KX = NUMBER OF TRACE SAMPLES I4 00000160 CA INPUT V = DIAGONAL MATRIX WHICH NORMALIZES R4 00000170 CA MAIN DIAGONAL OF OUTPUT MATRIX 00000180 CA OUTPUT PR = OUTPUT MATRIX. UPPER TRIANGULAR R4 00000190 CA PART ONLY. STORED IN PACKED FORM 00000200 CA BY COLUMNS 00000210 CA INPUT KEY = INPUT FLAG I4 00000220 CA 1 - PREPARE PR FOR COMPUTING RC'S 00000230 CA ONLY 00000240 CA 2 - PREPARE PR FOR COMPUTING WAVELET 00000250 CA ONLY 00000260 CA 3 - PREPARE PR FOR COMPUTING BOTH 00000270 CA 00000280 CA 00000290 CA THIS ROUTINE COMPUTES THE POSITIVE DEFINITE 00000300 CA SYMMETRIC MATRIX WHICH MUST BE INVERTED FOR EACH 00000310 CA ITERATION OF THE DDLOG PROCESS. 00000320 CA 00000330 CA 00000340 C 00000350 C SUBROUTINES CALLED: NONE 00000360 C 00000370 C EJECT 00000380 C 00000390 SUBROUTINE MPRDCT (Z, KR, KW, KX, V, PR, KEY) 00000400 CAEND 00000410 C 00000420 DOUBLE PRECISION A 00000430 DIMENSION Z(1),PR(1),V(1) 00000440 IA=KW+1 00000450 IF((KW/2)*2.NE.KW)IA=KW 00000460 IQ=0 00000470 KRQ=0 00000480 IF(KEY.EQ.2) GO TO 70 00000490 C 00000500 DO 30 00000510 * M=1,KR 00000520 M2=2*M 00000530 ILM=M2-IA+1 00000540 IHM=KW+M2-IA 00000550 IF(ILM.LT.1) ILM=1 00000560 IF(IHM.GT.KX)IHM=KX 00000570 INDXM=IA-M2+KR 00000580 C 00000590 DO 20 00000600 * L=1,M 00000610 IQ=IQ+1 00000620 A=0.D0 00000630 L2=2*L 00000640 IL=L2-IA+1 00000650 IH=KW+L2-IA 00000660 IF(IL.LT.1) IL=1 00000670 IF(IH.GT.KX) IH=KX 00000680 IF(ILM.GT.IL) IL=ILM 00000690 IF(IHM.LT.IH) IH=IHM 00000700 IF(IH.LT.IL) GO TO 20 00000710 INDXL=IA-L2+KR-INDXM 00000720 IL=IL+INDXM 00000730 IH=IH+INDXM 00000740 C 00000750 DO 10 00000760 * I=IL,IH 00000770 C 00000780 10 A=A+DBLE(Z(I))*DBLE(Z(I+INDXL)) 00000790 C 00000800 20 PR(IQ)=A*DBLE(V(L))*DBLE(V(M)) 00000810 C 00000820 30 CONTINUE 00000830 C 00000840 IF(KEY.EQ.1) RETURN 00000850 IQQ=0 00000860 C 00000870 DO 60 00000880 * M=1,KW 00000890 IQ=IQ+IQQ 00000900 IQQ=IQQ+1 00000910 INDXMM=IA-M 00000920 ILM=1 00000930 INDXM=INDXMM/2 00000940 IF(INDXM*2.EQ.INDXMM) ILM=2 00000950 C 00000960 DO 50 00000970 * L=1,KR 00000980 IQ=IQ+1 00000990 A=0.D0 00001000 L2=2*L 00001010 IL=L2-IA+1 00001020 IH=KW+L2-IA 00001030 IF(IL.LT.1) IL=1 00001040 IF(IH.GT.KX)IH=KX 00001050 INDXL=IA-L2+KR 00001060 J=IL+ILM+1 00001070 IL=((J/2)*2)/J+IL 00001080 IF(IL.GT.IH) GO TO 50 00001090 INDXL=IL+INDXL 00001100 ILQ=IL+1 00001110 IH=INDXM+(IH+((ILQ/2)*2)/ILQ)/2 00001120 IL=ILQ/2+INDXM 00001130 C 00001140 DO 40 00001150 * I=IL,IH 00001160 A=A+DBLE(Z(I))*DBLE(Z(INDXL)) 00001170 C 00001180 40 INDXL=INDXL+2 00001190 C 00001200 50 PR(IQ)=A*DBLE( V(L))*DBLE(V(M+KR)) 00001210 C 00001220 60 CONTINUE 00001230 C 00001240 KRQ=KR 00001250 IQ=(KR*(KR+1))/2 00001260 C 00001270 70 DO 100 00001280 * M=1,KW 00001290 IQ=IQ+KRQ 00001300 INDXM=IA-M 00001310 C 00001320 DO 90 00001330 * L=1,M 00001340 IQ=IQ+1 00001350 A=0.D0 00001360 J=M+L 00001370 IF((J/2)*2.NE.J)GO TO 90 00001380 INDXL=IA-L 00001390 IC=INDXL/2 00001400 KXQ=KX/2 00001410 IF(IC*2.NE.INDXL) KXQ=(KX+1)/2 00001420 IL=IC+1 00001430 IH=KXQ+IC 00001440 ID=INDXM/2-IC 00001450 C 00001460 DO 80 00001470 * I=IL,IH 00001480 C 00001490 80 A=A+DBLE(Z(I))*DBLE(Z(I+ID)) 00001500 C 00001510 90 PR(IQ)=A*DBLE(V(L+KR))*DBLE(V(M+KR)) 00001520 C 00001530 100 CONTINUE 00001540 C 00001550 RETURN 00001560 END 00001570