C 00000030 C PROGRAM 'PX' COMPUTES A WEIGHT MATRIX (PX) OF (PHI,LAMDA) GIVEN A 00000040 C VARIANCE-COVARIANCE MATRIX OF (X,Y,Z) 00000050 C 00000060 C INPUT: CARD 1; JOB TITLE (FORMAT 20A4) 00000070 C 00000080 C CARD 2; NUMBER OF STATIONS 00000090 C INPUT UNIT NUMBER FOR SIGMA(X,Y,Z) 00000100 C FORMAT OF SIGMA(X,Y,Z) 0=UPPER TRIANGULAR PORTION ON00000110 C 1=FULL MATRIX 00000120 C OUTPUT FORMAT OF SIGMA(PHI,LAMDA) 00000130 C 0=LINE PRINTER ONLY 00000140 C .NE.0=AUXILLARY DEVICE (7=PUNCHED CARDS) 00000150 C FORMAT FOR AUXILLARY OUTPUT 0=UPPER TRIANG. PORTION 00000160 C 1=FULL MATRIX 00000170 C FORMAT(5I3) 00000180 C 00000190 C CARD 3; NAME OF ELLIPSOID 00000200 C SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID 00000210 C SEMI-MINOR AXIS OF REFERENCE ELLIPSOID 00000220 C TRANSLATION COMPONENTS (DX,DY,DZ) 00000230 C SCALE FOR COORDINATES IN PPM (EG. FOR 1 PPM INPUT 00000240 C IS 1.0) 00000250 C FORMAT(5A4,2F10.0,3F8.0,F6.0) 00000260 C 00000270 C CARD TYPE 4; SIGMA(X,Y,Z) IN METRES **2 ROW BY ROW FORMAT(8F10.0) 00000280 C 00000290 C CARD TYPE 5; STATION NUMBER,STATION NAME,X,Y,Z (GEODETIC OR 00000300 C GEOCENTRIC) 00000310 C FORMAT (I9,1X,2A8,4X,3F15.0) 00000320 C 00000330 C NOTES:-CARD TYPE 3 IS REPEATED UNTIL ENTIRE MATRIX IS INPUT 00000340 C -TO INPUT GEODETIC COORDINATES SET DX=DY=DZ=0.0 00000350 C -CARD TYPE 5 IS REPEATED ONCE FOR EACH STATION 00000360 C -PROGRAM WILL TAKE A MAXIMUM OF 25 STATIONS--FOR MORE STATIONS 00000370 C INCREASE DIMENSIONS AND PARAMETER 'IRD' TO 3*NUMBER OF STATION00000380 C 00000390 C 00000400 C 00000410 IMPLICIT REAL*8(A-H,O-Z) 00000420 REAL*8 NAME 00000430 REAL*4 SPHI,SLAM 00000440 DIMENSION X(75),DLONG(25),DLAT(25),H(25),C(75,75),C1(3,3),SIG(75, 00000450 175),CSIG(75,75),SIGX(75,75),CO(75,75),TITLE(20),NUMB(25),NAME(25, 00000460 22),ELNAME(5),RO(3) 00000470 IRD=75 00000480 C 00000490 C READ CARD 1 00000500 C 00000510 READ(5,1000)TITLE 00000520 C 00000530 C READ CARD 2 00000540 C 00000550 READ(5,1010)NST,IU,IFMT,IOUT,IOUTF 00000560 NNN=3*NST 00000570 N2=2*NST 00000580 IF(IFMT.NE.0)GO TO 30 00000590 C 00000600 C READ SIGMA(X,Y,Z) IN UPPER TRIANGULAR FORM AND FILL IN LOWER HALF 00000610 C (CARD TYPE 3) 00000620 C 00000630 DO 10 I=1,NNN 00000640 10 READ(IU,1020)(SIGX(I,J),J=I,NNN) 00000650 DO 20 I=1,NNN 00000660 K=I+1 00000670 IF(K.GT.NNN)GO TO 50 00000680 DO 20 J=K,NNN 00000690 20 SIGX(J,I)=SIGX(I,J) 00000700 GO TO 50 00000710 C 00000720 C READ FULL SIGMA(X,Y,Z) (CARD TYPE 3) 00000730 C 00000740 30 DO 40 I=1,NNN 00000750 40 READ(IU,1020)(SIGX(J,I),J=1,NNN) 00000760 C 00000770 C READ CARD 4 00000780 C 00000790 50 READ(5,1025)ELNAME,A,B,RO(1),RO(2),RO(3),SCALE 00000800 C 00000810 C READ CARD TYPE 5 SCALE AND TRANSLATE COORDINATES 00000820 C 00000830 K=1 00000840 DO 60 I=1,NST 00000850 READ(5,1030)NUMB(I),NAME(I,1),NAME(I,2),X(K),X(K+1),X(K+2) 00000860 M=1 00000870 L=K+3 00000880 DO 55 J=K,L 00000890 X(J)=((DABS(X(J))+DABS(X(J))*SCALE/1.0D6)*X(J)/DABS(X(J)))-RO(M) 00000900 55 M=M+1 00000910 60 K=L 00000920 C 00000930 C COMPUTE GEODETIC (CURVILINEAR) COORDINATES AND TRANSFORMATION MATRIX00000940 C 00000950 DO 70 I=1,NNN 00000960 DO 70 J=1,NNN 00000970 70 C(J,I)=0.0D0 00000980 DO 80 I=1,NST 00000990 II=3*I-2 00001000 IJ=3*I-1 00001010 IK=3*I 00001020 CALL TRACOT(A,B,X(II),X(IJ),X(IK),DLAT(I),DLONG(I),H(I),C1) 00001030 DO 80 K=1,3 00001040 DO 80 KK=1,3 00001050 80 C(3*I-3+K,3*I-3+KK)=C1(K,KK) 00001060 90 CONTINUE 00001070 C 00001080 C COMPUTE SIG(PHI,LAMBDA,H) 00001090 C 00001100 CALL MATMUL(C,IRD,IRD,SIGX,IRD,IRD,CSIG,IRD,IRD,NNN,NNN,NNN,1) 00001110 DO 100 I=1,NNN 00001120 DO 100 J=1,NNN 00001130 100 CO(J,I)=C(I,J) 00001140 CALL MATMUL(CSIG,IRD,IRD,CO,IRD,IRD,SIG,IRD,IRD,NNN,NNN,NNN,1) 00001150 C 00001160 C PRINT RESULTS 00001170 C 00001180 PRINT 1040,TITLE 00001190 PRINT 1045,ELNAME,A,B,RO(1),RO(2),RO(3),SCALE 00001200 NN=1 00001210 DO 130 I=1,NST 00001220 CALL RADARC(DLAT(I),K,L,SPHI) 00001230 CALL RADARC(DLONG(I),M,N,SLAM) 00001240 PRINT 1050,NUMB(I),NAME(I,1),NAME(I,2),X(NN),X(NN+1),X(NN+2),K,L, 00001250 1SPHI,M,N,SLAM,H(I) 00001260 130 NN=NN+3 00001270 C 00001280 C PRINT INPUT SIGMA(X,Y,Z) AND ITS CORRELATION COEFFICIENT MATRIX 00001290 C 00001300 PRINT 1060 00001310 CALL MATRIT(SIGX,NNN,NNN,IRD) 00001320 CALL CORCOF(SIGX,NNN,CO,IRD) 00001330 PRINT 1070 00001340 CALL MATRIT(CO,NNN,NNN,IRD) 00001350 C 00001360 C PRINT SIGMA(PHI,LAMDA,H) AND ITS CORRELATION COEFFICIENT MATRIX 00001370 C 00001380 PRINT 1080 00001390 CALL MATRIT(SIG,NNN,NNN,IRD) 00001400 CALL CORCOF(SIG,NNN,CO,IRD) 00001410 PRINT 1090 00001420 CALL MATRIT(CO,NNN,NNN,IRD) 00001430 C 00001440 C ELIMINATE PORTION OF SIGMA(PHI,LAMDA,H) PERTAINING TO H AND PRINT 00001450 C 00001460 CALL SQUASH(SIG,IRD,NNN) 00001470 PRINT 1100 00001480 CALL MATRIT(SIG,N2,N2,IRD) 00001490 C 00001500 C INVERT SIGMA(PHI,LAMDA) AND PRINT 00001510 C 00001520 CALL CHOLD(SIG,IRD,N2,DETA,&999) 00001530 PRINT 1110 00001540 CALL MATRIT(SIG,N2,N2,IRD) 00001550 IF(IOUT.EQ.0)GO TO 999 00001560 IF(IOUTF.NE.0)GO TO 150 00001570 IF(IOUT.EQ.6)PRINT 1110 00001580 C 00001590 C WRITE PX(PHI,LAMDA) ON AUXILLARY DEVICE IF IOUT.NE.0 00001600 C UPPER TRIANGULAR PORTION ONLY 00001610 C 00001620 DO 140 I=1,N2 00001630 140 WRITE(IOUT ,1120)(SIG(I,J),J=I,N2) 00001640 GO TO 999 00001650 C 00001660 C WRITE PX(PHI,LAMDA) ON AUXILLARY DEVICE FULL MATRIX 00001670 C 00001680 150 IF(IOUT.EQ.6)PRINT 1110 00001690 DO 160 I=1,N2 00001700 160 WRITE(IOUT,1120)(SIG(I,J),J=1,N2) 00001710 999 STOP 00001720 1000 FORMAT(20A4) 00001730 1010 FORMAT(5I3) 00001740 1020 FORMAT(8F10.0) 00001750 1025 FORMAT(5A4,2F10.0,3F8.0,F6.0) 00001760 1030 FORMAT(I9,1X,2A8,4X,3F15.0) 00001770 1040 FORMAT('1',25X,20A4) 00001780 1045 FORMAT('-',10X,5A4,3X,'A=',F10.2,3X,'B=',F10.2,3X,'DX=',F8.3,3X, 00001790 1'DY=',F8.3,3X,'DZ=',F8.3,3X,'SCALE=',F6.2,3X,'PPM' 00001800 2 ///3X,'STATION STATION',13X,'X',16X,'Y',00001810 316X,'Z',13X,'LATITUDE',11X,'LONGITUDE ELL HEIGHT'/4X,'NUMBER', 00001820 49X,'NAME',3(14X,'(M)'),47X,'(M)'//) 00001830 1050 FORMAT('0',1X,I9,2X,2A8,3F17.3,3X,I3,I4,F8.3,I7,I4,F8.3,F10.3) 00001840 1060 FORMAT(////30X,'VARIANCE COVARIANCE MATRIX (X,Y,Z) IN METRES **2')00001850 1070 FORMAT(////30X,'CORRELATION COEFFICIENT MATRIX (X,Y,Z)') 00001860 1080 FORMAT(////25X,'VARIANCE COVARIANCE MATRIX (PHI,LAMDA,H) IN SEC**200001870 1 AND METRES**2') 00001880 1090 FORMAT(////25X,'CORRELATION COEFFICIENT MATRIX OF (PHI,LAMDA,H)') 00001890 1100 FORMAT(////26X,'VARIANCE COVARIANCE MATRIX OF (PHI,LAMDA) IN SEC**00001900 12') 00001910 1110 FORMAT(////35X,'WEIGHT MATRIX PX OF (PHI,LAMDA)') 00001920 1120 FORMAT(8(F10.3)) 00001930 END 00001940 C 00001950 C SUBROUTINE 'RADARC' CONVERTS RADIANS TO DEGREES MINUTES AND 00001960 C SECONDS. FOR NEGATIVE ANGLES ONLY THE LEFTMOST NONZERO VALUE IS 00001970 C NEGATIVE (EGS. -50,15,30.5 ; 0,-35,30.0 ; 0,0,-50.5) 00001980 C 00001990 C NOTE: THE 0.0005 VALUE IS TO GUARD AGAINST ROUNDOFF 00002000 C 00002010 C INPUT: A = RADIAN VALUE OF ANGLE (REAL*8) 00002020 C 00002030 C OUTPUT: I = DEGREES (INTEGER) 00002040 C J = MINUTES (INTEGER) 00002050 C S = SECONDS (REAL*4) 00002060 C 00002070 SUBROUTINE RADARC(A,I,J,S) 00002080 DOUBLE PRECISION A,SEC,AD,AJ,RHO, SIGN 00002090 DATA RHO/206264.8062470963D0/ 00002100 C 00002110 C CHECK SIGN OF 'A' -- SET SIGN=-1 IF NEGATIVE AND CONVERT 'A' TO 00002120 C POSITIVE VALUE 00002130 C 00002140 SIGN=1.0D0 00002150 IF(A.LT.0.0)SIGN=-1.0D0 00002160 IF(SIGN.LT.0.0)A=-A 00002170 C 00002180 C CONVERT 'A' TO ARCSECONDS 00002190 C 00002200 SEC=A*RHO+0.0005D0 00002210 C 00002220 C FIND INTEGER DEGREES 00002230 C 00002240 I=SEC/3600.0D0 00002250 AD=I 00002260 C 00002270 C FIND INTEGER MINUTES 00002280 C 00002290 J=SEC/60.0D0-AD*60.0D0 00002300 AJ=J 00002310 C 00002320 C FIND REAL*4 SECONDS 00002330 C 00002340 S=SEC-AD*3600.0D0-AJ*60.0D0-0.0005D0 00002350 C 00002360 C SET LEFTMOST VALUE NEGATIVE IF SIGN=-1 00002370 C 00002380 IF(I.NE.0)GO TO 20 00002390 IF(J.EQ.0)GO TO 10 00002400 J=J*SIGN 00002410 GO TO 30 00002420 10 S=S*SIGN 00002430 GO TO 30 00002440 20 I=I*SIGN 00002450 C 00002460 C CONVERT 'A' BACK TO NEGATIVE IF SIGN=-1 00002470 C 00002480 30 IF(SIGN.LT.0.0)A=-A 00002490 RETURN 00002500 END 00002510 SUBROUTINE SQUASH(A,IRDA,N) 00002520 C 00002530 C SUBROUTINE 'SQUASH' TAKES A DOUBLE PRECISION MATRIX 'A' AND 00002540 C ELIMINATES EVERY THIRD ROW AND COLUMN - INTENDED USE IS TO CONVERT 00002550 C VARIANCE-COVARIANCE MATRIX OF (PHI,LAMDA,H) ((OR PX MATRIX)) TO 00002560 C VARIANCE COVARIANCE (OR PX)) MATRIX OF (PHI,LAMDA) ONLY 00002570 C INPUT: MATRIX 'A' TO BE REDUCED 00002580 C IRDA IS ROW DIMENSION OF 'A' 00002590 C N IS SIZE OF 'A' IN ANY PARTICULAR PROBLEM 00002600 C OUTPUT: MATRIX 'A' THE REDUCED MATRIX ((IN PLACE OF ORIGINAL)) 00002610 REAL*8 A(IRDA,N) 00002620 NM=N-N/3 00002630 NN=N-1 00002640 C 00002650 C ELIMINATE EVERY THIRD COLUMN 00002660 C 00002670 DO 30 I=1,N 00002680 K=1 00002690 L=4 00002700 LL=5 00002710 10 DO 20 J=L,LL 00002720 20 A(I,J-K)=A(I,J) 00002730 K=K+1 00002740 L=L+3 00002750 LL=LL+3 00002760 IF(LL.LT.N)GO TO 10 00002770 30 CONTINUE 00002780 C 00002790 C ELIMINATE EVERY THIRD ROW 00002800 C 00002810 DO 60 J=1,NM 00002820 K=1 00002830 L=4 00002840 LL=5 00002850 40 DO 50 I=L,LL 00002860 50 A(I-K,J)=A(I,J) 00002870 K=K+1 00002880 L=L+3 00002890 LL=LL+3 00002900 IF(LL.LT.N )GO TO 40 00002910 60 CONTINUE 00002920 RETURN 00002930 END 00002940 SUBROUTINE TRACOT(A,B,X,Y,Z,PHI,DAMBDA,H,D) 00002950 C***********************************************************************00002960 C* TRANSFORMATION OF THE TRIPLET (X,Y,Z) TO (PHI,LAMBDA,H), AND THE *00002970 C* FORMATION OF THE TRANSFORMATION MATRIX D *00002980 C***********************************************************************00002990 IMPLICIT REAL*8 (A-H,M-Z) 00003000 REAL*8 D(3,3) 00003010 ESQ=(A**2-B**2)/A**2 00003020 P=DSQRT(X**2+Y**2) 00003030 H=0.D0 00003040 PHI1=DATAN((Z/P)*(1.D0/(1.D0-ESQ))) 00003050 1 N=A/DSQRT(1.D0-ESQ*(DSIN(PHI1)**2)) 00003060 H=(P/DCOS(PHI1))-N 00003070 PHI=DATAN((Z/P)*(1.D0/(1.D0-ESQ*(N/(N+H))))) 00003080 IF((PHI-PHI1).LE.1.D-10) GO TO 2 00003090 PHI1=PHI 00003100 GO TO 1 00003110 2 H=(P/DCOS(PHI))-N 00003120 DAMBDA=DATAN2(Y,X) 00003130 C 00003140 C COMPUTE TRANSFORMATION MATRIX D 00003150 C 00003160 DS=DSIN(PHI) 00003170 DC=DCOS(PHI) 00003180 SL=DSIN(DAMBDA) 00003190 CL=DCOS(DAMBDA) 00003200 N=A/DSQRT(1.D0-ESQ*(DS**2)) 00003210 M=N*((1.D0-ESQ)/(1.D0-ESQ*(DS**2))) 00003220 MM=M+H 00003230 NN=(N+H)*DC 00003240 RHO=206264.8062D0 00003250 D(1,1)=-DS*CL/MM*RHO 00003260 D(1,2)=-DS*SL/MM*RHO 00003270 D(1,3)=DC/MM*RHO 00003280 D(2,1)=-SL/NN*RHO 00003290 D(2,2)=CL/NN*RHO 00003300 D(2,3)=0.D0 00003310 D(3,1)=DC*CL 00003320 D(3,2)=DC*SL 00003330 D(3,3)=DS 00003340 4 RETURN 00003350 END 00003360 SUBROUTINE CORCOF (SIG,N,COEF,NR) 00003370 C***********************************************************************00003380 C* CORCOF COMPUTES THE CORRELATION COEFFICIENT MATRIX OF THE VARIANCE *00003390 C* COVARIANCE MATRIX SIG, WHOSE DIMENSIONS ARE N , AND NR IS THE *00003400 C* MAXIMUM NUMBER OF ROWS IN SIG ( OR COEF ) *00003410 C***********************************************************************00003420 IMPLICIT REAL*8 (A-H,O-Z) 00003430 REAL*8 SIG(NR,N),COEF(NR,N) 00003440 I=1 00003450 J=1 00003460 1 COEF(I,J)=SIG(I,J)/(DSQRT(SIG(I,I))*DSQRT(SIG(J,J))) 00003470 IF(J.NE.I) COEF(J,I)=COEF(I,J) 00003480 J=J+1 00003490 IF(J.LE.N) GO TO 1 00003500 I=I+1 00003510 IF(I.GT.N) GO TO 2 00003520 J=I 00003530 GO TO 1 00003540 2 CONTINUE 00003550 RETURN 00003560 END 00003570 SUBROUTINE CHOLD(A,IRDA,NA,DETA,*) 00003580 C 00003590 C MATRIX INVERSION USING CHOLESKI DECOMPOSITION 00003600 C 00003610 C INPUT ARGUMENTS 00003620 C A = ARRAY CONTAINING POSITIVE DEFINITE SYMMETRIC INPUT MATRIX 00003630 C IRDA = ROW DIMENSION OF ARRAY CONTAINING INPUT MATRIX 00003640 C NA = SIZE OF INPUT MATRIX 00003650 C OUTPUT ARGUMENTS 00003660 C A = CONTAINS INVERSE OF INPUT MATRIX (INPUT DESTROYED) 00003670 C DETA = DETERMINANT OF INPUT MATRIX 00003680 C * = ERROR RETURN (TAKEN IF NA .LT. 1 OR IF DETA .LT. SING) 00003690 C 00003700 DOUBLE PRECISION A,DETA,SUM,SQRT,DSQRT,ABS,DABS,SING 00003710 DIMENSION A(IRDA,NA) 00003720 SQRT(SUM) = DSQRT(SUM) 00003730 ABS(DETA) = DABS(DETA) 00003740 DATA SING/1D-50/ 00003750 C 00003760 C CHOLESKI DECOMPOSITION OF INPUT MATRIX INTO TRIANGULAR MATRIX 00003770 C 00003780 IF(NA .LT. 1) GO TO 18 00003790 DETA = A(1,1) 00003800 A(1,1) = SQRT(A(1,1)) 00003810 IF(NA .EQ. 1) GO TO 6 00003820 DO 1 I = 2,NA 00003830 1 A(I,1) = A(I,1) / A(1,1) 00003840 DO 5 J = 2,NA 00003850 SUM = 0. 00003860 J1 = J - 1 00003870 DO 2 K = 1,J1 00003880 2 SUM = SUM + A(J,K) ** 2 00003890 DETA = DETA * (A(J,J) - SUM) 00003900 A(J,J) = SQRT(A(J,J) - SUM) 00003910 IF(J .EQ. NA) GO TO 5 00003920 J2 = J + 1 00003930 DO 4 I = J2,NA 00003940 SUM = 0. 00003950 DO 3 K = 1,J1 00003960 3 SUM = SUM + A(I,K) * A(J,K) 00003970 4 A(I,J) = (A(I,J) - SUM) / A(J,J) 00003980 5 CONTINUE 00003990 6 IF(ABS(DETA) .LT. SING) GO TO 16 00004000 C 00004010 C INVERSION OF LOWER TRIANGULAR MATRIX 00004020 C 00004030 DO 7 I = 1,NA 00004040 7 A(I,I) = 1. / A(I,I) 00004050 IF(NA .EQ. 1) GO TO 10 00004060 N1 = NA - 1 00004070 DO 9 J = 1,N1 00004080 J2 = J + 1 00004090 DO 9 I = J2,NA 00004100 SUM = 0. 00004110 I1 = I - 1 00004120 DO 8 K = J,I1 00004130 8 SUM = SUM + A(I,K) * A(K,J) 00004140 9 A(I,J) = - A(I,I) * SUM 00004150 C 00004160 C CONSTRUCTION OF INVERSE OF INPUT MATRIX 00004170 C 00004180 10 DO 15 J = 1,NA 00004190 IF(J .EQ. 1) GO TO 12 00004200 J1 = J - 1 00004210 DO 11 I = 1,J1 00004220 11 A(I,J) = A(J,I) 00004230 12 DO 14 I = J,NA 00004240 SUM = 0. 00004250 DO 13 K = I,NA 00004260 13 SUM = SUM + A(K,I) * A(K,J) 00004270 14 A(I,J) = SUM 00004280 15 CONTINUE 00004290 RETURN 00004300 16 WRITE(6,17) DETA 00004310 17 FORMAT(10X, 'SINGULAR MATRIX IN CHOLD. DET =',E20.5) 00004320 RETURN 1 00004330 18 WRITE(6,19) 00004340 19 FORMAT(10X,'MATRIX OF DIMENSION ZERO IN CHOLD') 00004350 RETURN 1 00004360 END 00004370 SUBROUTINE MATRIT(A,M,N,NZ) 00004380 C 00004390 C THIS SUBROUTINE PRINTS DOUBLE PRECISION MATRICES 5 COLUMNS AND 40 ROWS00004400 C PAGE 00004410 C A IS THE MATRIX TO BE PRINTED 00004420 C M IS THE ROW DIMENSION OF A 00004430 C N IS THE COLUMN DIMENSION OF A 00004440 C NZ IS THE MAXIMUM ROW DIMENSION OF A IN MAIN 00004450 C IN CASE OF MORE THAN 40 ROWS, THE MATRIX COLUMN CONTINUES IN THE NEXT00004460 C OF THE SAME PAGE 00004470 C A HEADING SHOULD BE PRINTED IN THE MAIN PROGRAM . THE MATRIX WILL STA00004480 C PRINTING ON THE SAME PAGE 00004490 C 00004500 IMPLICIT REAL*8(A-H,O-Z) 00004510 DIMENSION A(NZ,N),IR(2),NC1(5) 00004520 C 00004530 C INITIALIZE 'INT' TO 0 TO INDICATE THAT OUTPUT IS NOT TO START ON A NEW00004540 C 00004550 INT=0 00004560 NC=1 00004570 IF(M.GT.40) GO TO 1 00004580 C 00004590 C PROCEDURE TO WRITE MATRIX IF NO. OF ROWS < OR = 40 00004600 C 00004610 NEW=0 00004620 C 00004630 C PRINT COLUMN NUMBERS 00004640 C 00004650 6 NEW=NEW+5 00004660 IF(N.LT.NEW) NEW=N 00004670 IF (INT.EQ.0) WRITE(6,100)(I,I=NC,NEW) 00004680 IF (INT.NE.0) WRITE(6,200)(I,I=NC,NEW) 00004690 C 00004700 C PRINT M ROWS OF MATRIX 00004710 C 00004720 DO 4 I=1,M 00004730 4 WRITE(6,300) I,(A(I,K),K=NC,NEW) 00004740 IF (N.EQ.NEW) RETURN 00004750 INT=1 00004760 NC=NC+5 00004770 GO TO 6 00004780 C 00004790 C PROCEDURE TO PRINT MATRIX IF NO. OF ROWS IS > 40 00004800 C DETERMINE THE NUMBER OF PAGE COLUMNS EACH MATRIX COLUMN IS TO OCCUPY (00004810 C 00004820 1 NCL=M/40 00004830 IF(MOD(M,40).NE.0) NCL=NCL+1 00004840 MN=1 00004850 NPER=0 00004860 LEFT=0 00004870 25 IF(NCL-5) 83,84,85 00004880 83 NPER=5-NCL 00004890 84 MCL=NCL 00004900 GO TO 86 00004910 85 LEFT=NCL-5 00004920 MCL=5 00004930 86 IR(1)=1 00004940 24 IR(2)=1 00004950 ICL=MCL 00004960 NREM=NPER 00004970 JI=1 00004980 C 00004990 C DETERMINE AND PRINT COLUMN NUMBERS 00005000 C 00005010 12 DO 8 I=JI,ICL 00005020 8 NC1(I)=NC 00005030 21 IF (NREM.EQ.0.OR.NC.EQ.N) GO TO 9 00005040 NC=NC+1 00005050 JI=ICL+1 00005060 IF(NREM-NCL)10,10,11 00005070 10 ICL=NREM+ICL 00005080 LEFT=NCL-NREM 00005090 NREM=0 00005100 GO TO 12 00005110 11 ICL=NCL+ICL 00005120 NREM=NREM-NCL 00005130 GO TO 12 00005140 9 JI=ICL 00005150 IF(INT.EQ.0) WRITE(6,600)(NC1(I),I=1,JI) 00005160 IF(INT.NE.0) WRITE(6,700)(NC1(I),I=1,JI) 00005170 35 ICL=MCL 00005180 22 NREM=NPER 00005190 IF(MN.EQ.0) GO TO 81 00005200 MN=0 00005210 GO TO 82 00005220 81 IF(MOD((IR(1)-1),40).EQ.0.AND.MOD((IR(2)-1),40).EQ.0) GO TO 17 00005230 82 NR=IR(1) 00005240 IR(1)=IR(1)+1 00005250 JI=1 00005260 NM=0 00005270 C 00005280 C PRINT MATRIX VALUES 00005290 C 00005300 16 DO 7 I=JI,ICL 00005310 IF(NR.GT.M) GO TO 18 00005320 GO TO (30,31,32,33,34),I 00005330 30 WRITE(6,400) NR,A(NR,NC1(I)) 00005340 NM=1 00005350 GO TO 7 00005360 31 IF(NM.NE.0) WRITE(6,500) NR,A(NR,NC1(I)) 00005370 IF(NM.EQ.0) WRITE(6,410) NR,A(NR,NC1(I)) 00005380 GO TO 7 00005390 32 WRITE(6,510) NR,A(NR,NC1(I)) 00005400 GO TO 7 00005410 33 WRITE(6,520) NR,A(NR,NC1(I)) 00005420 GO TO 7 00005430 34 WRITE(6,530) NR,A(NR,NC1(I)) 00005440 7 NR=NR+40 00005450 GO TO 28 00005460 18 IF(NC1(ICL).EQ.N.AND.MOD(NR,40).EQ.0) RETURN 00005470 28 IF(NREM.EQ.0.OR.NC1(ICL).EQ.N) GO TO 35 00005480 JI=ICL+1 00005490 IF(NREM-NCL) 14,14,15 00005500 14 ICL=NREM+ICL 00005510 NREM=0 00005520 IF (NCL.GT.2) GO TO 26 00005530 NR=IR(2)-1 00005540 GO TO 16 00005550 26 NR=IR(2) 00005560 IR(2)=IR(2)+1 00005570 GO TO 16 00005580 15 ICL=NCL+ICL 00005590 NREM=NREM-NCL 00005600 NR=IR(2) 00005610 IR(2)=IR(2)+1 00005620 GO TO 16 00005630 17 INT=1 00005640 NPER=0 00005650 IF(LEFT.EQ.0) GO TO 20 00005660 MCL=LEFT 00005670 IR(1)=NR-39 00005680 MN=1 00005690 IF(MCL.LE.5) GO TO 19 00005700 MCL=5 00005710 LEFT=LEFT-5 00005720 GO TO 24 00005730 19 NPER=5-LEFT 00005740 GO TO 24 00005750 20 NC=NC+1 00005760 MN=1 00005770 GO TO 25 00005780 100 FORMAT('0',14X,I3,4(18X,I3)) 00005790 200 FORMAT('1',14X,I3,4(18X,I3)) 00005800 300 FORMAT(' ',I3,5(3X,D19.12)) 00005810 400 FORMAT(' ',I3,1X,D19.12) 00005820 410 FORMAT(' ',23X,I3,1X,D19.12) 00005830 500 FORMAT('+',23X,I3,1X,D19.12) 00005840 510 FORMAT('+',46X,I3,1X,D19.12) 00005850 520 FORMAT('+',69X,I3,1X,D19.12) 00005860 530 FORMAT('+',92X,I3,1X,D19.12) 00005870 600 FORMAT('0',11X,I3,4(20X,I3)) 00005880 700 FORMAT('1',11X,I3,4(20X,I3)) 00005890 END 00005900 SUBROUTINE MATMUL(A,IRDA,ICDA,B,IRDB,ICDB,C,IRDC,ICDC,M,N,K,ICODE)00005910 C 00005920 C SUBROUTINE MATMUL MULTIPLIES MATRICES COLUMN BY COLUMN IN ONE OF THRE00005930 C WAYS: A*B=C (ICODE=1) 00005940 C : A*B=A (ICODE=2) 00005950 C : A*B=B (ICODE=3) 00005960 C FOR REAL*8 MATRICES UP TO 1000 * 1000 00005970 C 00005980 C IN ALL CASES IT CHECKS THAT THE SIZE OF A MATRIX IS NOT GREATER THAN 00005990 C ITS DIMENSION AND IF THE INDICATED MULTIPLICATION IS POSSIBLE 00006000 C 00006010 C INPUT: A,IRDA,ICDA MATRIX A AND ITS DECLARED DIMENSIONS 00006020 C B,IRDB,ICDB MATRIX B AND ITS DECLARED DIMENSIONS 00006030 C M*N IS SIZE OF A 00006040 C N*K IS SIZE OF B 00006050 C ICODE =TYPE OF MULTIPLICATION DESIRED (SEE ABOVE) 00006060 C 00006070 C OUTPUT: M*K MATRIX C 00006080 C 00006090 C NOTE: IT IS IMPORTANT THAT THE COLUMN DIMENSIONS AS WELL AS THE ROW00006100 C DIMENSIONS BE SPECIFIED SINCE THEY ARE REQUIRED IN CHECKING 00006110 C IF THE INDICATED MULTIPLICATION IS POSSIBLE 00006120 C NOTE: IF ICODE .NE. 1 INPUT MATRICES ARE DESTROYED 00006130 C 00006140 C EXAMPLE CALL STATEMENTS 00006150 C ICODE=1 CALL MATMUL(A,10,15,B,10,25,C,25,25,10,10,15,1) 00006160 C ICODE=2 CALL MATMUL(A,10,15,B,10,25,A,10,15,10,10,10,2) 00006170 C ICODE=3 CALL MATMUL(A,10,15,B,10,25,B,10,25,10,10,15,3) 00006180 C 00006190 C C. CHAMBERLAIN FEB./1976 00006200 C 00006210 DOUBLE PRECISION A(IRDA,ICDA),B(IRDB,ICDB),C(IRDC,ICDC),SUM, 00006220 1WORK(1000) 00006230 C 00006240 C CHECK THAT SIZES ARE NOT GREATER THAN DIMENSIONS 00006250 C 00006260 IF(M.GT.IRDA.OR.N.GT.ICDA)GO TO 997 00006270 IF(N.GT.IRDB.OR.K.GT.ICDB)GO TO 997 00006280 IF(M.GT.IRDC.OR.K.GT.ICDC)GO TO 997 00006290 GO TO (10,40,80),ICODE 00006300 C 00006310 C CODING FOR A*B=C (ICODE=1) 00006320 C 00006330 10 DO 20 I=1,K 00006340 DO 20 J=1,M 00006350 20 C(J,I)=0.0D0 00006360 DO 30 I=1,K 00006370 DO 30 L=1,N 00006380 DO 30 J=1,M 00006390 30 C(J,I)=C(J,I)+A(J,L)*B(L,I) 00006400 RETURN 00006410 C 00006420 C CODING FOR A*B=A (ICODE=2) 00006430 C 00006440 40 IF(K.GT.ICDA)GO TO 998 00006450 DO 70 I=1,M 00006460 DO 60 L=1,K 00006470 SUM=0.0D0 00006480 DO 50 J=1,N 00006490 50 SUM=SUM+A(I,J)*B(J,L) 00006500 60 WORK(L)=SUM 00006510 DO 70 J=1,K 00006520 70 C(I,J)=WORK(J) 00006530 RETURN 00006540 C 00006550 C CODING FOR A*B=B (ICODE=3) 00006560 C 00006570 80 IF(M.GT.IRDB)GO TO 999 00006580 DO 110 I=1,K 00006590 DO 100 L=1,M 00006600 SUM=0.0D0 00006610 DO 90 J=1,N 00006620 90 SUM=SUM+A(L,J)*B(J,I) 00006630 100 WORK(L)=SUM 00006640 DO 110 J=1,M 00006650 110 C(J,I)=WORK(J) 00006660 RETURN 00006670 C 00006680 C WRITE ERROR MESSAGE AND STOP 00006690 C 00006700 997 PRINT 1000 00006710 STOP 00006720 998 PRINT 1010 00006730 STOP 00006740 999 PRINT 1020 00006750 STOP 00006760 1000 FORMAT('1',10X,'ERROR IN MATMUL ==> SIZE OF MATRIX > DIMENSION OF00006770 1 MATRIX') 00006780 1010 FORMAT('1',10X,'ERROR IN MATMUL ==> SIZE OF A*B > DIMENSION OF A 00006790 1(ICODE=2)') 00006800 1020 FORMAT('1',10X,'ERROR IN MATMUL ==> SIZE OF A*B > DIMENSION OF B 00006810 1(ICODE=3)') 00006820 END 00006830