C C* * 0020 C* MATRIX SUBROUTINE PACKAGE * 0030 C* * 0040 C* DAVID WELLS ( WITH HELP FROM CAPT. J.G. GIBBONS) * 0050 C* * 0060 C* UNIVERSITY OF NEW BRUNSWICK * 0070 C* * 0080 C* MARCH 30, 1971 * 0090 C* * 0100 C* * 0110 C* * 0120 C*********************************************************************** 0130 C 0140 C******************ROUTINES USING ONE MATRIX*************** 0150 C 0160 C MATINR(A,I,J,NA) READ A(I,J) BY ROWS 0170 C MATINC(A,I,J,NA) READ A(I,J) BY COLUMNS 0180 C MATOUT(A,I,J,NA) PRINT A(I,J) 0190 C MATZER(A,I,J,NA) SET A(I,J)=0 0200 C MIDENT(A,I,NA) SET SQUARE A(I,J)=IDENTITY 0210 C MATINV(A,S,I,NA) SET A(I,J)=INVERSE OF SQUARE A(I,J) 0220 C SET S=DETERMINANT OF A(I,J) 0230 C 0240 C******************ROUTINES USING TWO MATRICES************* 0250 C 0260 C MCOPY (A,B,I,J,NA,NB) SET B(I,J)=A(I,J) 0270 C MCOPYC(A,B,I,J,NA,NB) SET B(I,J)=-A(I,J) 0280 C TRNSPZ(A,B,I,J,NA,NB) SET B(I,J)=A(J,I) 0290 C SCAMLT(S,A,B,I,J,NA,NB) SET B(I,J)=S*A(I,J) 0300 C MEIGEN(A,B,I,NA,NB) SET B(I,J)=ORTHOGONAL MATRIX OF 0310 C EIGENVECTORS OF SYMMETRIC A(I,J) 0320 C SET A(I,J)=DIAGONAL MATRIX WHOSE 0330 C ELEMENTS ARE EIGENVALUES 0340 C 0350 C******************ROUTINES USING THREE MATRICES*********** 0360 C 0370 C MATADD(A,B,C,I,J,NA,NB,NC) SET C(I,J)=A(I,J)+B(I,J) 0380 C MATSUB(A,B,C,I,J,NA,NB,NC) SET C(I,J)=A(I,J)-B(I,J) 0390 C MATMLT(A,B,C,I,J,K,NA,NB,NC) SET C(I,J)=A(I,K)*B(K,J) 0400 C 0410 C******************NOTATION******************************** 0420 C 0430 C A,B,C, ARE MATRICES 0440 C NA,NB,NC ARE DIMENSIONED COLUMN LENGTHS OF A,B,C 0450 C I,J ARE THE ROW,COLUMN SIZE TO BE OPERATED ON 0460 C S IS A SCALAR 0470 C 0480 SUBROUTINE MATRIX(A,B,C) 0490 IMPLICIT REAL*8(A-H,O-Z) 0500 DIMENSION A(1),B(1),C(1) 0510 RETURN 0520 C 0530 C READ A(I,J) BY ROWS 0540 C 0550 ENTRY MATINR(A,IR,JC,NA) 0560 DO 2 I=1,IR 0570 DO 2 J=1,JC 0580 IA=I+(J-1)*NA 0590 2 A(IA)=FREAD(K) 0600 RETURN 0610 C 0620 C READ A(I,J) BY COLUMNS 0630 C 0640 ENTRY MATINC(A,IR,JC,NA) 0650 DO 3 J=1,JC 0660 DO 3 I=1,IR 0670 IA=I+(J-1)*NA 0680 3 A(IA)=FREAD(K) 0690 RETURN 0700 C 0710 C PRINT A(I,J) 0720 C 0730 ENTRY MATOUT(A,IR,JC,NA) 0740 PRINT 6 0750 K2=0 0760 4 K1=K2 0770 K3=K2+1 0780 K2=K2+5 0790 IF(JC.LE.K2) K2=JC 0800 PRINT 7,(K,K=K3,K2) 0810 DO 5 I=1,IR 0820 N1=K1*NA+I 0830 N2=(K2-1)*NA+I 0840 5 PRINT 8,I,(A(K),K=N1,N2,NA) 0850 IF(JC.GT.K2) GOTO 4 0860 PRINT 6 0870 6 FORMAT('0 LNCT= ',I5) 7 FORMAT(1H ,2X,5(15X,I5)) 0890 8 FORMAT(5H ROW ,I5,5D20.10) 0900 RETURN 0910 C 0920 C SET A(I,J)=0 0930 C 0940 ENTRY MATZER(A,IR,JC,NA) 0950 DO 9 I=1,IR 0960 DO 9 J=1,JC 0970 IA=I+(J-1)*NA 0980 9 A(IA)=0.D0 0990 RETURN 1000 C 1010 C SET A(I,J)=IDENTITY MATRIX 1020 C 1030 ENTRY MIDENT(A,IR,NA) 1040 DO 10 I=1,IR 1050 DO 10 J=1,JC 1060 IA=I+(J-1)*NA 1070 A(IA)=0.D0 1080 IF(I.EQ.J) A(IA)=1.D0 1090 10 CONTINUE 1100 RETURN 1110 C 1120 C SET B(I,J)=A(I,J) 1130 C 1140 ENTRY MCOPY (A,B,IR,JC,NA,NB) 1150 DO 11 I=1,IR 1160 DO 11 J=1,JC 1170 IA=I+(J-1)*NA 1180 IB=I+(J-1)*NB 1190 11 B(IB)=A(IA) 1200 RETURN 1210 C 1220 C SET B(I,J)= - A(I,J) 1230 C 1240 ENTRY MCOPYC(A,B,IR,JC,NA,NB) 1250 DO 12 I=1,IR 1260 DO 12 J=1,JC 1270 IA=I+(J-1)*NA 1280 IB=I+(J-1)*NB 1290 12 B(IB)= -A(IA) 1300 RETURN 1310 C 1320 C SET B(J,I)= A(I,J) 1330 C 1340 ENTRY TRNSPZ(A,B,IR,JC,NA,NB) 1350 DO 13 I=1,IR 1360 DO 13 J=1,JC 1370 IA=I+(J-1)*NA 1380 IB=J+(I-1)*NB 1390 13 B(IB)=A(IA) 1400 RETURN 1410 C 1420 C SET B(I,J)= S * A(I,J) 1430 C 1440 ENTRY SCAMLT(S,A,B,IR,JC,NA,NB) 1450 DO 14 I=1,IR 1460 DO 14 J=1,JC 1470 IA=I+(J-1)*NA 1480 IB=I+(J-1)*NB 1490 14 B(IB)=S*A(IA) 1500 RETURN 1510 C 1520 C SET C(I,J)=A(I,J) - B(I,J) 1530 C 1540 ENTRY MATSUB(A,B,C,IR,JC,NA,NB,NC) 1550 SIGN=-1.D0 1560 GOTO 15 1570 C 1580 C SET C(I,J)= A(I,J) + B(I,J) 1590 C 1600 ENTRY MATADD(A,B,C,IR,JC,NA,NB,NC) 1610 SIGN=1.D0 1620 15 DO 16 I=1,IR 1630 DO 16 J=1,JC 1640 IA=I+(J-1)*NA 1650 IB=I+(J-1)*NB 1660 IC=I+(J-1)*NC 1670 16 C(IC)=A(IA)+B(IB)*SIGN 1680 RETURN 1690 C 1700 C SET C(I,J) = A(I,K) * B(K,J) 1710 ENTRY MATMLT(A,B,C,IR,JC,KK,NA,NB,NC) 1720 DO 17 I=1,IR 1730 DO 17 J=1,JC 1740 IC=I+(J-1)*NC 1750 C(IC)=0.D0 1760 DO 17 K=1,KK 1770 IA=I+(K-1)*NA 1780 IB=K+(J-1)*NB 1790 17 C(IC)=C(IC)+A(IA)*B(IB) 1800 RETURN 1810 END 1820 C 1830 C SET A(I,J) = INVERSE OF A(I,J) 1840 C SET DET = DETERMINANT OF A(I,J) 1850 C 1860 SUBROUTINE MATINV(A,DET,NORDER,NDIM) 1870 IMPLICIT REAL*8(A-H,O-Z) 1880 DIMENSION A(NDIM,NORDER),Y(100),ICJ(100),IRR(100),JCC(100) 1890 IF(NDIM.LE.100) GOTO 20 1900 PRINT 10 1910 10 FORMAT(5X,'MATRIX LARGER THAN 100*100') 1920 RETURN 1930 20 DET=1.D0 1940 DO 205 K=1,NORDER 1950 PIVOT=0.D0 1960 C SET PIVOT=MAXIMUM ELEMENT LEFT 1970 DO 100 I=1,NORDER 1980 DO 100 J=1,NORDER 1990 IF(K.EQ.1) GOTO 50 2000 C DONT USE THE SAME PIVOT ROW OR COLUMN TWICE 2010 KA=K-1 2020 DO 30 IA=1,KA 2030 IF(I.EQ.IRR(IA)) GOTO 100 2040 IF(J.EQ.JCC(IA)) GOTO 100 2050 30 CONTINUE 2060 C PIVOT=LARGEST ELEMENT 2070 50 IF(DABS(A(I,J)).LE.PIVOT) GOTO 100 2080 PIVOT=DABS(A(I,J)) 2090 IRR(K)=I 2100 JCC(K)=J 2110 100 CONTINUE 2120 C RETURN IF PIVOT.LT.EPS 2130 IF(PIVOT.GT.0.1D-50) GOTO 120 2140 PRINT 110 2150 110 FORMAT('PIVOT ELEMENT TOO SMALL') 2160 DET=0.D0 2170 RETURN 2180 120 IR=IRR(K) 2190 JC=JCC(K) 2200 PIVOT=A(IR,JC) 2210 DET=DET*PIVOT 2220 C NORMALIZE PIVOTAL ROW (ITH ROW) 2230 DO 130 J=1,NORDER 2240 130 A(IR,J)=A(IR,J)/PIVOT 2250 A(IR,JC)=1./PIVOT 2260 IF(NORDER.EQ.1) RETURN 2270 C CLEAR ROW AND COLUMN 2280 DO 205 I=1,NORDER 2290 AJCK=A(I,JC) 2300 IF(I.EQ.IR) GOTO 205 2310 A(I,JC)=-AJCK/PIVOT 2320 DO 200 J=1,NORDER 2330 IF(J.EQ.JC) GOTO 200 2340 A(I,J)=A(I,J)-AJCK*A(IR,J) 2350 200 CONTINUE 2360 205 CONTINUE 2370 C COUNT INTERCHANGES NECESSARY TO REORDER 2380 DO 210 I=1,NORDER 2390 IR=IRR(I) 2400 JC=JCC(I) 2410 210 ICJ(IR)=JC 2420 NA=NORDER-1 2430 DO 220 I=1,NA 2440 IA=I+1 2450 DO 220 J=IA,NORDER 2460 IF(ICJ(J).GE.ICJ(I)) GOTO 220 2470 TEMP=ICJ(J) 2480 ICJ(J)=ICJ(I) 2490 ICJ(I)=TEMP 2500 DET=-DET 2510 220 CONTINUE 2520 C REARRANGE COLUMNS 2530 DO 240 J=1,NORDER 2540 DO 230 I=1,NORDER 2550 IR=IRR(I) 2560 JC=JCC(I) 2570 230 Y(JC)=A(IR,J) 2580 DO 240 I=1,NORDER 2590 240 A(I,J)=Y(I) 2600 C REARRANGE ROWS 2610 DO 260 I=1,NORDER 2620 DO 250 J=1,NORDER 2630 IR=IRR(J) 2640 JC=JCC(J) 2650 250 Y(IR)=A(I,JC) 2660 DO 260 J=1,NORDER 2670 260 A(I,J)=Y(J) 2680 RETURN 2690 END 2700 $GO SUBROUTINE MEIGEN(A,R,N,NA,NR) 2840 C 2710 C EIGENVALUES AND EIGENVECTORS BY JACOBI METHOD 2720 C SEE CARNAHAN PAGES 250-260 2730 C 2740 C INPUT REAL SYMMETRIC MATRIX B(N,N) 2750 C OUTPUT EIGENVECTORS AS COLUMNS OF ORTHOGONAL MATRIX R(N,N) 2760 C 2770 C ALGORITHM ITERATIVELY REDUCES INPUT MATRIX A TO NEAR DIAGONAL FORM 2780 C BY SUCCESSSIVELY ANNIHILATING OFF-DIAGONAL ELEMENTS 2790 C USING A SEQUENCE OF ROTATIONS. THE MATRIX R IS THE FINAL 2800 C PRODUCT ROTATION MATRIX, AND OUTPUT MATRIX A CONTAINS THE 2810 C EIGENVALUES AS DIAGONAL ELEMENTS. 2820 C 2830 IMPLICIT REAL*8(A-H,O-Z) 2850 DIMENSION A(NA,N),R(NR,N),Y(100) 2860 C SPECIFY ITMAX. MAXIMUM NUMBER OF ITERATIONS 2870 C EPS1. ROTATION ANGLE TEST(ROTATE BY PI/4 IF 2880 C EQUAL DIAGONAL ELEMENTS) 2890 C EPS2. IGNORE OFF-DIAGONAL ELEMENTS .LT.EPS2 2900 C EPS3. BREAKOUT TEST (SUM OF SQUARES OF DIAGONAL 2910 C ELEMENTS UNCHANGED BY LAST ITERATION) 2920 C 2930 ITMAX=50 2940 EPS1=0.1D-9 2950 EPS2=0.1D-14 EPS3=0.1D-4 2970 NM1=N-1 2980 C SET R=IDENTITY AND UPPER TRIANGULATE A 2990 DO 16 I=1,N 3000 DO 16 J=1,N 3010 R(I,J)=0. 3020 IF(I.EQ.J) R(I,J)=1. 3030 IF(I.GT.J) A(I,J)=0. 3040 16 CONTINUE 3050 C COMPUTE SIGMA1 AND S 3060 SIGMA1=0. 3070 OFFDSQ=0. 3080 DO 17 I=1,N 3090 SIGMA1=SIGMA1+A(I,I)*A(I,I) 3100 IP1=I+1 3110 IF(I.GE.N) GOTO 18 3120 DO 17 J=IP1,N 3130 17 OFFDSQ=OFFDSQ+A(I,J)*A(I,J) 3140 18 S=2.0*OFFDSQ+SIGMA1 3150 C 3160 C BEGIN JACOBI ITERATION 3170 DO 29 ITER=1,ITMAX 3180 DO 27 I=1,NM1 3190 IP1=I+1 3200 DO 27 J=IP1,N 3210 Q=DABS(A(I,I)-A(J,J)) 3220 C 3230 C COMPUTE SINE AND COSINE OF ROTATION ANGLE 3240 IF(Q.LE.EPS1) GOTO 19 3250 IF(DABS(A(I,J)).LE.EPS2) GOTO 27 3260 P=2.0*A(I,J)*Q/(A(I,I)-A(J,J)) 3270 SPQ=DSQRT(P*P+Q*Q) 3280 CSA=DSQRT((1.0+Q/SPQ)/2.0) 3290 SNA=P/(2.0*CSA*SPQ) 3300 GOTO 20 3310 19 CSA=1.0/DSQRT(2.0D0) 3320 SNA=CSA 3330 20 CONTINUE 3340 C 3350 C UPDATE COLUMNS I AND J OF R - EQUIVALENT TO 3360 C MULTIPLICATION BY THE ANNIHILATION MATRIX 3370 DO 21 K=1,N 3380 HOLDKI=R(K,I) 3390 R(K,I)=HOLDKI*CSA+R(K,J)*SNA 3400 21 R(K,J)=HOLDKI*SNA-R(K,J)*CSA 3410 C 3420 C COMPUTE NEW ELEMENTS OF A IN ROWS I AND J 3430 DO 24 K=1,N 3440 IF(K.GT.J) GOTO 23 3450 Y(K)=A(I,K) 3460 A(I,K)=CSA*Y(K)+SNA*A(K,J) 3470 IF(K.NE.J) GOTO 22 3480 A(J,K)=SNA*Y(K)-CSA*A(J,K) 3490 22 GOTO 24 3500 23 HOLDIK=A(I,K) 3510 A(I,K)=CSA*HOLDIK+SNA*A(J,K) 3520 A(J,K)=SNA*HOLDIK-CSA*A(J,K) 3530 24 CONTINUE 3540 C 3550 C COMPUTE NEW ELEMENTS OF A IN COLUMNS I AND J 3560 Y(J)=SNA*Y(I)-CSA*Y(J) 3570 C 3580 C WHEN K IS LARGER THAN I 3590 DO 26 K=1,J 3600 IF(K.LE.I) GOTO 25 3610 A(K,J)=SNA*Y(K)-CSA*A(K,J) 3620 GOTO 26 3630 25 HOLDKI=A(K,I) 3640 A(K,I)=CSA*HOLDKI+SNA*A(K,J) 3650 A(K,J)=SNA*HOLDKI-CSA*A(K,J) 3660 26 CONTINUE 3670 27 A(I,J)=0. 3680 C 3690 C FIND SIGMA2 FOR TRANSFORMED A AND TEST FOR CONVERGENCE 3700 SIGMA2=0. 3710 DO 28 I=1,N 3720 28 SIGMA2=SIGMA2+A(I,I)*A(I,I) 3730 IF(1.0-SIGMA1/SIGMA2.LT.EPS3) RETURN 3740 29 SIGMA1=SIGMA2 3750 WRITE(6,203) ITER,S,SIGMA1,SIGMA2 3760 203 FORMAT('0 EIGEN NOT CONVERGED'/' ITER=',I5,/' S=',F10.5,/ 1 ' SIGMA1=',F10.5,/'SIGMA2=',F10.5) RETURN 3790 END 3800 C*********************************************************************** 0690 C* * 0700 C* ******* * 0710 C* *FREAD* * 0720 C* ******* * 0730 C* * 0740 C* FUNCTION FOR READING FROM CARDS FORMAT FREE * 0750 C* * 0760 C* R.L.PARKER (SEE BEDFORD INSTITUTE COMPUTER * 0770 C* NOTE 1969-8-C FOR DETAILS) * 0780 C* * 0790 C* MODIFICATIONS BY D.E.WELLS UNB JAN 1971 * 0800 C* * 0810 C* *********** * 0820 C* *ARGUMENTS* * 0830 C* *********** * 0840 C* * 0850 C* INPUT IN =0 NO PRINT * 0860 C* =1 PRINT MESSAGE AND BAD NUMBERS ONLY * 0870 C* =2 PRINT ALL * 0880 C* =3 REREAD LAST NUMBER AND PRINT * 0890 C* * 0900 C* OUTPUT IOUT=0 IF SUCCESSFUL READ * 0910 C* = FIRST BAD CHARACTER IF ERROR HIT * 0920 C* * 0930 C* COLUMN COUNTER ICOL AVAILABLE AS ARGUMENT * 0940 C* USUALLY LEFT UNTOUCHED BY USER * 0950 C* HOWEVER FOR FLEXIBILITY IT CAN BE SET WHEN FUNCTION CALLED * 0960 C* ICOL=81 WILL SKIP TO END OF CURRENT CARD * 0970 C* ICOL.LT.81 WILL SKIP BACK OR FORWARD TO THAT COLUMN * 0980 C* * 0990 C* ************************ * 1000 C* *POSSIBLE MODIFICATIONS* * 1010 C* ************************ * 1020 C* * 1030 C* IF 'IN' NOT REQUIRED - * 1040 C* DROP FROM PARAMETER LISTS IN FUNCTION STATEMENT AND IN ALL * 1050 C* CALLING STATEMENTS * 1060 C* INITIALIZE TO DESIRED VALUE ( 0 1 2 OR 3) USING EITHER * 1070 C* DATA STATEMENT ( E.G. 'DATA IN/1/ ' ) OR * 1080 C* REPLACEMENT STATEMENT (E.G. 'IN=2 ' ) * 1090 C* * 1100 C* IF ' IOUT' NOT REQUIRED - * 1110 C* DROP FROM PARAMETER LISTS IN FUNCTION STATEMENT AND IN ALL * 1120 C* CALLING STATEMENTS * 1130 C* * 1140 C* IF 'ICOL' NOT REQUIRED - * 1150 C* DROP FROM PARAMETER LISTS IN FUNCTION STATEMENT AND IN ALL * 1160 C* CALLING STATEMENTS * 1170 C* INITIALIZE TO 81 USING DATA STATEMENT ( DATA ICOL/81/ ) * 1180 C* * 1190 C* NB - FORTRAN FUNCTIONS MUST HAVE AT LEAST ONE ARGUMENT * 1200 C* * 1210 C* IF DOUBLE PRECISION DATA EXPECTED (MORE THAN SEVEN DIGITS) - * 1220 C* ADD TO BOTH FUNCTION AND TO CALLING PROGRAM THE STATEMENT * 1230 C* ' DOUBLE PRECISION FREAD ' * 1240 C* * 1250 C*********************************************************************** 1260 FUNCTION FREAD(IOUT) DOUBLE PRECISION FREAD 1280 DOUBLE PRECISION X,XX,B 1290 DIMENSION IA(80),IDIGIT(17) 1300 DATA IDIGIT/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9, 1310 1 1H ,1H, ,1H.,1HE,1HD,1H-,1H+/ 1320 DATA IE/3/ 1330 DATA IN/0/,ICOL/81/ C 1340 C ...TEST FOR REREAD MODE (IN=3) 1350 C 1360 IF(IN.GE.3) GOTO(400,410,860),IE 1370 IOUT=0 1380 C 1390 C ...TEST FOR EMPTY CARD BUFFER. READ NEW CARD IF EMPTY 1400 C 1410 IF(ICOL.GT.80) GOTO 201 1420 C 1430 C ...IE=1 UNTIL FIRST 'E' OR 'D' FOUND. THEN 1440 C IE=NUMBER OF 'E' OR 'D' CHARACTERS FOUND + 1 1450 C IC1=SUBSCRIPT OF FIRST CHARACTER IN FREAD CALL 1460 C B=POWER OF 10 BY WHICH FRACTIONAL PART MUST BE DIVIDED 1470 C 1480 100 IE=1 1490 IC1=ICOL 1500 B=1.0 1510 C 1520 C ...RETURN HERE TO PICK UP EXPONENT 1530 C IDEC=1 UNTIL FIRST '.' FOUND. THEN 1540 C IDEC=NUMBER OF '.' + NUMBER OF 'E' OR 'D' FOUND + 2 1550 C X=ABSOLUTE VALUE OF NUMBER EVALUATED (FREAD) 1560 C SIGN=SIGN OF NUMBER EVALUATED (FREAD) 1570 C NC=LENGTH OF LEGAL CHARACTER LIST (IDIGIT). EITHER 15 OR 17 1580 C 1590 110 IDEC=1 1600 X=0.0 1610 SIGN=+1.0 1620 NC=17 1630 C 1640 C ...IGNORE LEADING SPACES AND COMMAS 1650 C 1660 DO 150 IC=ICOL,80 1670 IF(IA(IC).NE.IDIGIT(11).AND.IA(IC).NE.IDIGIT(12)) GOTO 202 1680 150 CONTINUE 1690 C 1700 C ...END OF CARD. ERROR IF EXPONENTIAL PART EXPECTED 1710 C 1720 ICOL=ICOL-1 1730 IF(IE.GT.1) GOTO 800 1740 C 1750 C ...NO ERROR. READ ANOTHER CARD 1760 C 1770 201 READ 2000,(IA(J),J=1,80) 1780 2000 FORMAT(80A1) 1790 ICOL=1 1800 GO TO 100 1810 C 1820 C ...NO MORE LEADING SPACES OR COMMAS 1830 C 1840 C START TESTING CHARACTER BY CHARACTER UNTIL TERMINATING 1850 C SPACE OR COMMA FOUND 1860 C 1870 202 DO 290 ICOL=IC,80 1880 JA=IA(ICOL) 1890 C 1900 C ...COMPARE EACH CHARACTER TO LIST OF LEGAL CHARACTERS (IDIGIT) 1910 C 1920 DO 205 I=1,NC 1930 IF(JA.EQ.IDIGIT(I)) GOTO 206 1940 205 CONTINUE 1950 C ...ERROR - CHARACTER NOT ON IDIGIT LIST 1960 C 1970 GO TO 800 1980 C ...CHARACTER ON LIST. TEST WHETHER NUMERAL 1990 C 2000 206 J=I-9 2010 IF(J.LE.0) GOTO 25 2020 GOTO(25,300,300,210,310,310,220,290),J 2030 C 2040 C ...CHARACTER = '.' ERROR IF PREVIOUS '.' OR 'E' OR 'D' FOUND 2050 C 2060 210 IDEC=IDEC+IE 2070 IF(IDEC.GT.2) GOTO 800 2080 GOTO 290 2090 C 2100 C ...CHARACTER = '-' RESET 'SIGN' AND SHORTEN IDIGIT LIST 2110 C 2120 220 SIGN=-1.0 2130 GO TO 290 2140 C 2150 C ...CHARACTER IS NUMERAL 2160 C 2170 25 DIGIT=I-1 2180 GOTO(251,252,800),IDEC 2190 C 2200 C ...ASSEMBLE INTEGRAL PART OF NUMBER 2210 C 2220 251 X=10.0D0*X+DIGIT 2230 GO TO 290 2240 C 2250 C ...ASSEMBLE FRACTIONAL PART OF NUMBER 2260 C 2270 252 B=B/10.0D0 2280 X=X+B*DIGIT 2290 C 2300 C ...CHARACTER = '+' SHORTEN IDIGIT LIST 2310 C 2320 290 NC=15 2330 C 2340 C ...NO TERMINATING SPACE OR COMMA BEFORE END OF CARD FOUND 2350 C 2360 ICOL=81 2370 C 2380 C ...CHARACTER IS TERMINATING SPACE OR COMMA. TEST FOR E OR F FORMAT 2390 C 2400 300 GO TO(400,410),IE 2410 C 2420 C ...TERMINATING CHARACTER IS 'E' OR 'D'. SET XX=MANTISSA 2430 C 2440 310 IE=IE+1 2450 XX=X*SIGN 2460 C 2470 C ...ERROR IF NO MANTISSA 2480 C 2490 IF(ICOL.EQ.IC) GOTO 800 2500 ICOL=ICOL+3-IE 2510 GO TO (110,110,800),IE 2520 C 2530 C ...TERMINATION OF F FORMAT 2540 C 2550 400 FREAD=X*SIGN 2560 GO TO 500 2570 C 2580 C ...TERMINATION OF E FORMAT 2590 C 2600 410 IX=X*SIGN 2610 FREAD=XX*10.0D0**IX 2620 C 2630 C ...PRINT IF IN=2 OR 3 2640 C 2650 500 IC2=ICOL-1 2660 IF (IN .GT.1) PRINT 5000,(IA(J),J=IC1,IC2) 2670 5000 FORMAT (10X,80A1) 2680 RETURN 2690 C 2700 C ...ERROR HANDLING 2710 C 2720 800 IOUT=IA(ICOL) 2730 C 2740 C ...SEARCH FOR NEXT TERMINATING BLANK OR COMMA 2750 C 2760 DO 810 I=ICOL,80 2770 IF(IA(I).EQ.IDIGIT(11).OR.IA(I).EQ.IDIGIT(12)) GO TO 815 2780 810 CONTINUE 2790 I=81 2800 815 ICOL=I 2810 C 2820 C ...PRINT ERROR MESSAGE IF N=2 OR 3 2830 C 2840 IF(IN.EQ.0) GOTO 850 2850 IC2=ICOL-1 2860 PRINT 820,IOUT,(IA(J),J=IC1,IC2) 2870 820 FORMAT(25H MALFORMED NUMBER DUE TO ,A1,4H IN ,80A1) 2880 C 2890 C ...EVALUATE FREAD AND RETURN 2900 C 2910 850 GO TO(400,410,400),IE 2920 C 2930 C ...ILLEGAL REREAD ERROR 2940 C 2950 860 PRINT 870 2960 870 FORMAT(25H ILLEGAL REREAD ATTEMPTED ) 2970 FREAD=0. 2980 RETURN 2990 END