C ALGORITHM 594, COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VO1983, P. 125-130. C PROGRAM MAIN(INPUT,OUTPUT,TAPE20 = INPUT,TAPE5 = OUTPUT) MAN 10 C MAN 20 C MAN 30 C MAN 40 C A PROGRAM FOR RELATIVE ERROR ANALYSIS MAN 50 C DATE JANUARY 6, 1981 UNIVERSITY OF ILLINOIS MAN 60 C CRAY RESEARCH, INC. MAN 70 C BELL LABORATORIES - NAPERVILLE MAN 80 C SANDIA NATIONAL LABORATORIES - ALBUQUERQUE MAN 90 C MAN 100 C AUTHORS J. L. LARSON, M. E. PASTERNAK, J. A. WISNIEWSKI MAN 110 C MAN 120 C MAN 130 C MAN 140 C MAN 150 C THIS IS THE MAIN PROGRAM WHICH COMPRISES THE RELATIVE MAN 160 C ERROR ANALYSIS PACKAGE. THE MAIN PROGRAM READS IN THE DATA AND MAN 170 C ALLOCATES STORAGE. IF NOT ENOUGH WORK SPACE HAS BEEN ALLOCATED, MAN 180 C THE PROGRAM PRINTS AN ERROR MESSAGE INDICATING THE AMOUNT OF MAN 190 C STORAGE NEEDED AND THEN TERMINATES. MAN 200 C MAN 210 C MAN 220 C IN ORDER TO RUN THE PACKAGE, THE USER SUPPLIES --- MAN 230 C MAN 240 C 1. AS INPUT DATA, THE OUTPUT FROM MILLER'S MINICOMPILER, WHICH MAN 250 C SPECIFIES THE PROGRAM TREE OF THE ALGORITHM BEING ANALYZED MAN 260 C (SEE REFERENCE 9). THE 7TH TO 80TH COLUMNS OF THE MAN 270 C FIRST LINE OF THIS INPUT ARE PRINTED AS A TITLE FOR THE MAN 280 C ANALYSIS. THE 4TH TO 80TH COLUMNS OF THE LAST LINE OF MAN 290 C THE INPUT ARE USED TO IDENTIFY THE DATA MAN 300 C SET BEING ANALYZED, AND ARE ALSO PRINTED. MAN 310 C MAN 320 C 2. THE DATA ASSOCIATED WITH A NUMBER OF DIFFERENT (POSSIBLY MAN 330 C ONLY ONE) ROUNDOFF ERROR ANALYSES ARE NEEDED. THESE DATA MAN 340 C ITEMS WERE DEFINED BY THE INPUT STATEMENTS OF THE MINI- MAN 350 C COMPILER PROGRAM IMPLEMENTING THE ALGORITHM. THE DATA MAN 360 C ITEMS ASSOCIATED WITH EACH PROBLEM ARE CONCATENATED, ONE MAN 370 C BEHIND THE OTHER, TO FORM THE OVERALL DATA STREAM. EACH MAN 380 C LINE OF THE DATA STREAM CONTAINS ONE DATA ITEM. EACH MAN 390 C ITEM IS READ USING FORTRAN FORMAT(E22.10). THIS MAN 400 C DATA STREAM IS TERMINATED BY A SENTINEL LINE HAVING A MAN 410 C NEGATIVE INTEGER IN ITS FIRST 3 COLUMNS. THIS LINE IS MAN 420 C USED TO MARK END-OF-DATA. MAN 430 C MAN 440 C THE INPUT ASSOCIATED WITH A SINGLE ERROR ANALYSIS HAS MAN 450 C THE FOLLOWING FORM --- MAN 460 C MAN 470 C - MAN 480 C ( PROGRAM TREE (OUTPUT FROM MINICOMPILER) MAN 490 C - MAN 500 C - MAN 510 C ( LIST OF DATA ITEMS, 1 PER LINE MAN 520 C - MAN 530 C - MAN 540 C ( -N --- SENTINEL LINE MAN 550 C - MAN 560 C MAN 570 C MAN 580 C MULTIPLE DATA SETS FOR A GIVEN MINICOMPILER OUTPUT MAY MAN 590 C ALSO BE ANALYZED WITHOUT INCLUDING A COPY OF THE PROGRAM MAN 600 C TREE FOR EACH PROBLEM. HERE, ONLY ONE COPY OF THE MAN 610 C PROGRAM TREE IS NEEDED FOR ALL OF THE DATA SETS FOR A GIVEN MAN 620 C PROBLEM. TO USE THIS OPTION, PLACE THE DATA SETS FOR A MAN 630 C SINGLE PROGRAM TREE, ONE AFTER THE OTHER, IMMEDIATELY MAN 640 C AFTER THE PROGRAM TREE, WITH AN EXTRA LINE PLACED MAN 650 C BETWEEN DATA SETS AS A DIVIDER. THE 4TH TO 80TH COLUMNS MAN 660 C OF EACH DIVIDER LINE ARE PRINTED AS A TITLE FOR THE DATA MAN 670 C SET FOLLOWING IT. THIS DIVIDER LINE MUST CONTAIN IN THE MAN 680 C FIRST 2 COLUMNS, FOR CHECKING PURPOSES, THE NUMBER OF MAN 690 C DATA ITEMS FOR A SINGLE ANALYSIS. IF THIS NUMBER IS NOT MAN 700 C CORRECT, AN ERROR MESSAGE IS PRINTED AND THE ENTIRE MAN 710 C ANALYSIS IS TERMINATED. THE LAST SET OF DATA FOR A GIVEN MAN 720 C PROGRAM TREE IS FOLLOWED BY A SEPARATOR LINE CONTAINING A MAN 730 C NEGATIVE INTEGER IN ITS FIRST 3 COLUMNS. MAN 740 C MAN 750 C THE FORM OF THE INPUT FILE FOR MULTIPLE DATA SETS IS MAN 760 C THE FOLLOWING. MAN 770 C MAN 780 C - MAN 790 C ( PROGRAM TREE MAN 800 C - MAN 810 C - MAN 820 C ( DATA SET 1 MAN 830 C - MAN 840 C ( M --- DIVIDER LINE MAN 850 C - MAN 860 C - MAN 870 C ( DATA SET 2 MAN 880 C - MAN 890 C - MAN 900 C ( -N --- SEPARATOR LINE MAN 910 C - MAN 920 C - MAN 930 C ( -N --- SENTINEL LINE MAN 940 C MAN 950 C MAN 960 C MULTIPLE ALGORITHMS MAY BE ANALYZED IN ONE RUN. TO MAN 970 C INDICATE THAT THIS OPTION IS IN USE, CHANGE THE FIRST MAN 980 C LINE OF MILLER'S OUTPUT SO THAT IT CONTAINS A NON-ZERO MAN 990 C INTEGER IN THE 4TH TO SIXTH COLUMNS. THE MAN 1000 C INPUT FOR THE ALGORITHMS ARE CONCATENATED TOGETHER, WITH MAN 1010 C A SENTINEL LINE PLACED AT THE VERY END. MAN 1020 C MAN 1030 C FOR 2 ALGORITHMS, EACH CONTAINING A SINGLE SET OF DATA MAN 1040 C ITEMS, THE STRUCTURE OF THE INPUT FILE IS THE FOLLOWING. MAN 1050 C MAN 1060 C - MAN 1070 C ( PROGRAM TREE MAN 1080 C - MAN 1090 C - MAN 1100 C ( DATA ITEMS MAN 1110 C - MAN 1120 C - MAN 1130 C ( PROGRAM TREE MAN 1140 C - MAN 1150 C - MAN 1160 C ( DATA ITEMS MAN 1170 C - MAN 1180 C - MAN 1190 C ( -N --- SENTINEL LINE MAN 1200 C - MAN 1210 C MAN 1220 C MAN 1230 C MAN 1240 C MULTIPLE ALGORITHMS WITH MULTIPLE DATA SETS FOR ANY AND ALL MAN 1250 C ALGORITHMS MAY BE RUN IN THE SAME ANALYSIS BY COMBINING THE MAN 1260 C INPUT FOR EACH PROBLEM AS DESCRIBED ABOVE. MAN 1270 C MAN 1280 C THE FORM OF THE INPUT FILE TO ANALYZE 2 ALGORITHMS, EACH MAN 1290 C WITH TWO SETS OF DATA IS AS FOLLOWS. MAN 1300 C MAN 1310 C - MAN 1320 C ( PROGRAM TREE MAN 1330 C - MAN 1340 C - MAN 1350 C ( DATA SET 1 MAN 1360 C - MAN 1370 C - MAN 1380 C ( M --- DIVIDER LINE MAN 1390 C - MAN 1400 C - MAN 1410 C ( DATA SET 2 MAN 1420 C - MAN 1430 C - MAN 1440 C ( -N --- SEPARATOR LINE MAN 1450 C - MAN 1460 C - MAN 1470 C ( PROGRAM TREE MAN 1480 C - MAN 1490 C - MAN 1500 C ( DATA SET 1 MAN 1510 C - MAN 1520 C - MAN 1530 C ( M --- DIVIDER LINE MAN 1540 C - MAN 1550 C - MAN 1560 C ( DATA SET 2 MAN 1570 C - MAN 1580 C - MAN 1590 C ( -N --- SEPARATOR LINE MAN 1600 C - MAN 1610 C - MAN 1620 C ( -N --- SENTINEL LINE MAN 1630 C MAN 1640 C MAN 1650 C MAN 1660 C MAN 1670 C MAN 1680 C MAN 1690 C THE SOFTWARE PRINTS --- MAN 1700 C MAN 1710 C 1. A LISTING OF THE USER-SUPPLIED INFORMATION-- MAN 1720 C A. MINICOMPILER PROGRAM TREE MAN 1730 C B. CONSTANTS USED WITHIN THE PROBLEM MAN 1740 C C. DATA ITEMS MAN 1750 C MAN 1760 C MAN 1770 C 2. A MESSAGE AS TO WHICH METHOD, HEURISTIC OR EXACT IS MAN 1780 C USED (SEE REFERENCES 4 AND 7). MAN 1790 C MAN 1800 C 3. THE NUMERIC OUTPUT VALUES OF THE PROBLEM, THE MAN 1810 C CONDITION NUMBER OF THE PROBLEM, AND THE CONDITION NUMBER MAN 1820 C OF THE ALGORITHM. MAN 1830 C A RELATIVE ERROR BOUND ON EACH OUTPUT VALUE IS ALSO MAN 1840 C PRINTED. MAN 1850 C MAN 1860 C 4. A VALUE THETA, WHICH IS A MEASURE OF THE MAXIMUM MAN 1870 C DISTANCE BETWEEN THE COMPUTED SOLUTION AND THE TRUE MAN 1880 C SOLUTION TO THE GIVEN PROBLEM WITH SLIGHTLY DIFFERENT MAN 1890 C DATA IS ALSO PRINTED. MAN 1900 C IF THE HEURISTIC METHOD IS USED, BOUNDS ON THE VALUE MAN 1910 C OF THETA ARE PRINTED. MAN 1920 C MAN 1930 C 5. THE AMOUNT OF STORAGE REQUIRED BY THE SOFTWARE. MAN 1940 C MAN 1950 C MORE SPECIFIC DETAILS CONCERNING THE OUTPUT PRINTED BY MAN 1960 C THE SOFTWARE CAN BE FOUND IN REFERENCES 4, 5 AND 7. MAN 1970 C MAN 1980 C MAN 1990 C MAN 2000 C EXTERNALS MAN 2010 C MAN 2020 C SEVERAL EXTERNAL SUBROUTINES ARE NEEDED TO RUN THIS PACKAGE. MAN 2030 C MAN 2040 C THE FOLLOWING BLAS ROUTINES ARE NEEDED (SEE REFERENCE 8) - MAN 2050 C ISAMAX MAN 2060 C SASUM MAN 2070 C SAXPY MAN 2080 C SCOPY MAN 2090 C SDOT MAN 2100 C SNRM2 MAN 2110 C SSCAL MAN 2120 C SSWAP MAN 2130 C MAN 2140 C TWO LINPACK ROUTINES ARE REFERENCED BY THE CODE (SEE MAN 2150 C REFERENCE 3). THEY ARE MAN 2160 C SQRDC MAN 2170 C SQRSL. MAN 2180 C MAN 2190 C THE FOLLOWING EISPACK ROUTINES (REFERENCE 10) ARE CALLED MAN 2200 C QZHES MAN 2210 C QZIT MAN 2220 C QZVAL MAN 2230 C QZVEC MAN 2240 C MAN 2250 C MAN 2260 C MAN 2270 C TREATMENT OF ZERO-VALUED DATA ITEMS MAN 2280 C MAN 2290 C MAN 2300 C IN ORDER TO USE THIS PROGRAM WITH PROBLEMS IN WHICH ANY OPERAND MAN 2310 C MAY HAVE THE VALUE ZERO(EXCEPT WHEN DIVIDING), THE PROBLEM MAN 2320 C IS PERTURBED SLIGHTLY TO SIMULATE THE RELATIVE ERROR OPERATION. A MAN 2330 C CALCULATED VALUE OF MACHINE EPSILON IS USED IN THIS SIMULATION. MAN 2340 C MAN 2350 C IF AN OPERAND EQUALS ZERO, A PERTURBATION IS MADE TO THE OPERANDS MAN 2360 C IN ORDER TO SIMULATE THE ROUNDING ERROR PROCESS. THE FOLLOWING MAN 2370 C PERTURBATIONS ARE MADE FOR THE GIVEN CASES, WHERE MAN 2380 C MAN 2390 C ARGR = RIGHT OPERAND FOR THE OPERATION INVOLVED MAN 2400 C ARGL = LEFT OPERAND FOR THE OPERATION INVOLVED. MAN 2410 C MAN 2420 C MAN 2430 C OPERATION VALUE OF OPERAND VALUE USED MAN 2440 C MAN 2450 C + 0.0 MAX(EPSILON,EPSILON*MAX(ARGL,ARGR) MAN 2460 C MAN 2470 C - 0.0 MAX(EPSILON,EPSILON*MAX(ARGL,ARGR) MAN 2480 C MAN 2490 C * 0.0 0.0 MAN 2500 C MAN 2510 C / ARGL = 0.0 0.0 MAN 2520 C MAN 2530 C / ARGR = 0.0 **ERROR MESSAGE-PROGRAM STOPS MAN 2540 C MAN 2550 C SQUARE ROOT .LT. 0.0 **ERROR MESSAGE-PROGRAM STOPS MAN 2560 C MAN 2570 C NEGATION 0.0 0.0 MAN 2580 C MAN 2590 C MAN 2600 C MAN 2610 C MAN 2620 C MAN 2630 C REFERENCES --- MAN 2640 C MAN 2650 C 1. BARRODALE, I. AND ROBERTS, F. D. K. AN IMPROVED ALGORITHM FOR MAN 2660 C DISCRETE L1 LINEAR APPROXIMATION. SIAM J. NUMER. ANAL. MAN 2670 C 10(1973), 839-848. MAN 2680 C 2. BARRODALE, I. AND ROBERTS, F. D. K. ALGORITHM 478. SOLUTION OF MAN 2690 C AN OVERDETERMINED SYSTEM OF EQUATIONS IN THE L1 NORM COMMUN. MAN 2700 C ACM 17(1974) 319-320. MAN 2710 C 3. DONGARRA, J. J., BUNCH, J. R., MOLER, C. B. AND STEWART, G. W. MAN 2720 C LINPACK USER'S GUIDE, SIAM PRESS, PHILADELPHIA, 1979. MAN 2730 C 4. LARSON, J. L. METHODS FOR AUTOMATIC ERROR ANALYSIS OF NUMERICALMAN 2740 C ALGORITHMS, PH.D. THESIS, UNIVERSITY OF ILLINOIS, DEPARTMENT OFMAN 2750 C COMPUTER SCIENCE TECHNICAL REPORT NO. UIUCDCS-R-78-916, MAN 2760 C APRIL 1978. MAN 2770 C 5. LARSON, J. L., PASTERNAK, M. E. AND WISNIEWSKI, J. A. MAN 2780 C SOFTWARE FOR RELATIVE ERROR ANALYSIS(TO APPEAR). MAN 2790 C 6. LARSON, J. L. AND SAMEH, A. H. EFFICIENT CALCULATION OF THE MAN 2800 C EFFECTS OF ROUNDOFF ERRORS, ACM TRANS. MATH. SOFTWARE 4, MAN 2810 C 3(SEPT. 1978), 228-236. MAN 2820 C 7. LARSON, J. L. AND SAMEH, A. H. ALGORITHMS FOR ROUNDOFF ERROR MAN 2830 C ANALYSIS --- A RELATIVE APPROACH, COMPUTING, 24, 4(1980), MAN 2840 C 275-297. MAN 2850 C 8. LAWSON, C. L., HANSON, R. J., KINCAID, D. R. AND KROGH, F. T. MAN 2860 C BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE, ACM TRANS. MAN 2870 C MATH. SOFTWARE 5, 3(SEPT 1979), 308-323. MAN 2880 C 9. MILLER, W. AND SPOONER, D. SOFTWARE FOR ROUNDOFF ERROR ANALYSISMAN 2890 C ACM TRANS. MATH. SOFTWARE 4, 4(DEC. 1978), 369-387. MAN 2900 C 10. SMITH, B.T., ET AL., MATRIX EIGENSYSTEM ROUTINES - EISPACK MAN 2910 C GUIDE, SECOND EDITION, SPRINGER-VERLAG, NEW YORK, (1976). MAN 2920 C MAN 2930 C MAN 2940 C MAN 2950 C MAN 2960 C***********************************************************************MAN 2970 C***********************************************************************MAN 2980 C***********************************************************************MAN 2990 C MAN 3000 C MAN 3010 C MAN 3020 DIMENSION TITLE(74), NAME(77) MAN 3030 C MAN 3040 INTEGER TOP(5), TOPTMP(3,5), TOPSAV(5) MAN 3050 C ARRAY FOR PROBLEM TITLE MAN 3060 C MAN 3070 C INDICES INTO RWORK FOR REAL ARRAYS MAN 3080 C MAN 3090 INTEGER ABAR, BBAR, AL1, BL1, AATORR, BBT, CONS, DATA, EL1, MAN 3100 * OUTORQ, OUTA, OUTB, LPRD, XL1, V, VALUE, WORK, COL MAN 3110 C MAN 3120 C INDICES INTO IWORK FOR INTEGER ARRAYS MAN 3130 C MAN 3140 INTEGER SL1, OP, LOP, ROP, OUTNDX, ACTIVE MAN 3150 C MAN 3160 C MAN 3170 INTEGER RWRKFR, IWRKFR, RWRKSZ, IWRKSZ, ISTORE, RSTORE MAN 3180 LOGICAL ERROR, FIRST, FIRSTF, FIRSTL MAN 3190 C MAN 3200 C MAN 3210 C MAN 3220 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ MAN 3230 C MAN 3240 C INSTALLATION INSTRUCTIONS MAN 3250 C MAN 3260 C TO USE THIS PACKAGE AT ANY SITE, SEVERAL INTERNAL PARAMETERS MAY MAN 3270 C NEED TO BE CHANGED. LINES THAT MAY NEED TO BE MODIFIED ARE MAN 3280 C BRACKETED WITH COMMENTS OF THE FORM... MAN 3290 C C+ INSTALL MAN 3300 C MAN 3310 C MAN 3320 C MAN 3330 C MAN 3340 C WORK ARRAYS TO ALLOCATE STORAGE MAN 3350 C MAN 3360 C+ INSTALL MAN 3370 REAL RWORK(5000) MAN 3380 INTEGER IWORK(5000) MAN 3390 C+ INSTALL MAN 3400 C MAN 3410 C RWORK AND IWORK ARE THE WORK STORAGE ARRAYS FROM WHICH MAN 3420 C STORAGE IS ALLOCATED. MAN 3430 C MAN 3440 C+ INSTALL MAN 3450 INTEGER INUNIT, OUTUNT MAN 3460 COMMON /IOINFO/ INUNIT, OUTUNT MAN 3470 C+ INSTALL MAN 3480 C MAN 3490 C INUNIT AND OUTUNT ARE THE DEVICE NUMBERS OF THE INPUT MAN 3500 C AND OUTPUT DEVICES, RESPECTIVELY. MAN 3510 C MAN 3520 C+ INSTALL MAN 3530 DATA RWRKSZ, IWRKSZ /5000,5000/ MAN 3540 C+ INSTALL MAN 3550 C MAN 3560 C RWRKSZ AND IWRKSZ ARE THE SIZES OF THE REAL AND INTEGER MAN 3570 C STORAGE ARRAYS, RESPECTIVELY. MAN 3580 C MAN 3590 C MAN 3600 C TO DETERMINE THE CHOICE CRITERION BETWEEN THE TWO METHODS, A LEAST MAN 3610 C SQUARES FIT TO SEVERAL SETS OF TIMING DATA WAS DONE. THE WORK DONE MAN 3620 C FOR THE EXACT METHOD WAS FOUND TO BE PROPORTIONAL TO THE CUBE OF MAN 3630 C THE NUMBER OF OUTPUT ITEMS TIMES THE NUMBER OF FACES (SEE REFERENCE MAN 3640 C 4). THE VALUE OF SWITCH WAS CHOSEN SO THAT THE EXACT METHOD WOULD MAN 3650 C CONSUME NO MORE THAN 30 SECONDS CPU TIME ON A CDC CYBER 175. THIS MAN 3660 C PARAMETER SHOULD BE APPROPRIATELY SCALED TO REPRESENT THE USER'S MAN 3670 C COMPUTING ENVIRONMENT. MAN 3680 C MAN 3690 C+ INSTALL MAN 3700 DATA SWITCH /3.0E7/ MAN 3710 C+ INSTALL MAN 3720 C MAN 3730 C SWITCH DECIDES WHICH METHOD TO USE. MAN 3740 C SET CORRECTLY FOR THE CYBER 175. MAN 3750 C MAN 3760 C MAN 3770 C MAN 3780 DATA ZERO, HALF, ONE, TWO /0.0E0,0.5E0,1.0E0,2.0E0/ MAN 3790 C MAN 3800 DATA MAXCOL, MAXROW /5,53/ MAN 3810 C MAN 3820 C USED TO PRINT OUT PROGRAM TREE MAN 3830 C MAXCOL = MAXIMUM NUMBER OF COLUMNS PER LINE MAN 3840 C MAXROW = NUMBER OF LINES PER PAGE OF PRINTOUT MAN 3850 C MAN 3860 C+ INSTALL MAN 3870 INUNIT = 20 MAN 3880 OUTUNT = 5 MAN 3890 C+ INSTALL MAN 3900 C MAN 3910 C SET THE INPUT AND OUTPUT DEVICE NUMBERS MAN 3920 C MAN 3930 C MAN 3940 C MAN 3950 C MAN 3960 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++MAN 3970 C MAN 3980 C MAN 3990 C CALCULATE MACHINE EPSILON = THE SMALLEST POSITIVE NUMBER SUCH MAN 4000 C THAT, 1 + MACHINE EPSILON .NE. 1 MAN 4010 C MAN 4020 C MAN 4030 EPS = ONE MAN 4040 10 CALL SDIFST(ONE, ONE+EPS, RESULT) MAN 4050 IF (RESULT.EQ.ZERO) GO TO 20 MAN 4060 EPS = EPS/TWO MAN 4070 GO TO 10 MAN 4080 C MAN 4090 20 CONTINUE MAN 4100 EPS = EPS*TWO MAN 4110 C MAN 4120 C MAN 4130 C INPUT DATA AND ALLOCATE STORAGE MAN 4140 C MAN 4150 30 CONTINUE MAN 4160 C MAN 4170 C THESE MUST BE RESET EACH TIME THROUGH LOOP MAN 4180 C MAN 4190 ISTORE = 0 MAN 4200 RSTORE = 0 MAN 4210 IWRKFR = 1 MAN 4220 RWRKFR = 1 MAN 4230 C MAN 4240 ERROR = .FALSE. MAN 4250 FIRST = .TRUE. MAN 4260 FIRSTF = .TRUE. MAN 4270 FIRSTL = .TRUE. MAN 4280 C MAN 4290 C READ OPERATIONS MAN 4300 C MAN 4310 READ (INUNIT,99991) NOP, IFLAG, (TITLE(I),I=1,74) MAN 4320 C MAN 4330 C NOP = NUMBER OF OPERATIONS IN PROGRAM TREE MAN 4340 C MAN 4350 C MAN 4360 C IF NOP=0, SKIP, READ IN NEXT PARAMETER MAN 4370 C MAN 4380 IF (NOP.EQ.0) GO TO 190 MAN 4390 C MAN 4400 C IF NOP .GT. 0, CONTINUE, ELSE, END-OF-FILE MAN 4410 C MAN 4420 IF (NOP.GT.0) GO TO 40 MAN 4430 STOP MAN 4440 C MAN 4450 40 CONTINUE MAN 4460 C MAN 4470 C WRITE OUT PROBLEM TITLE MAN 4480 C MAN 4490 WRITE (OUTUNT,99988) (TITLE(I),I=1,74) MAN 4500 ISTORE = ISTORE + 3*NOP MAN 4510 RSTORE = RSTORE + 3*NOP MAN 4520 C MAN 4530 C CHECKS FOR OVERFLOW OF WORK ARRAY DUE TO NOP MAN 4540 C IF OVERFLOW, READ IN DUMMY VARIABLES MAN 4550 C MAN 4560 IF (.NOT.((ISTORE.GT.IWRKSZ) .OR. (RSTORE.GT.RWRKSZ))) GO TO 50 MAN 4570 C MAN 4580 C WORKSPACE EXCEEDED, READ OVER DATA AND SET ERROR MAN 4590 C MAN 4600 C MAN 4610 C READS IN DUMMY VARIABLES MAN 4620 C MAN 4630 ERROR = .TRUE. MAN 4640 READ (INUNIT,99984) (IW1,IW2,IW3,I=1,NOP) MAN 4650 GO TO 180 MAN 4660 C MAN 4670 C ALLOCATE STORAGE FOR OPERATIONS AND READ IN MAN 4680 C MAN 4690 50 LOP = IWRKFR MAN 4700 OP = IWRKFR + NOP MAN 4710 ROP = IWRKFR + 2*NOP MAN 4720 IWRKFR = IWRKFR + 3*NOP MAN 4730 DO 60 I=1,NOP MAN 4740 ILOP = LOP + I - 1 MAN 4750 IOP = OP + I - 1 MAN 4760 IROP = ROP + I - 1 MAN 4770 READ (INUNIT,99984) IWORK(ILOP), IWORK(IOP), IWORK(IROP) MAN 4780 60 CONTINUE MAN 4790 C MAN 4800 C MAN 4810 C MAN 4820 C FORMAT PROGRAM TREE FOR OUTPUT IN COLUMNS, AND PRINT IT MAN 4830 C MAN 4840 C MAN 4850 WRITE (OUTUNT,99993) MAN 4860 WRITE (OUTUNT,99994) MAN 4870 C MAN 4880 MRMC = MAXROW*MAXCOL MAN 4890 NPAGE = 1 MAN 4900 NOPSAV = NOP MAN 4910 NDONE = 0 MAN 4920 INC = 0 MAN 4930 C MAN 4940 70 CONTINUE MAN 4950 C MAN 4960 C HOW MANY COLUMNS NEEDED MAN 4970 C MAN 4980 IF (NPAGE.NE.1) GO TO 80 MAN 4990 NCOLS = NOP/MAXROW MAN 5000 IF (NCOLS*MAXROW.NE.NOP) NCOLS = NCOLS + 1 MAN 5010 IF (NCOLS.LE.MAXCOL) GO TO 80 MAN 5020 C MAN 5030 C FILL ENTIRE PAGE IF .GT.1 PAGE WILL BE NEEDED MAN 5040 C MAN 5050 L = MAXROW MAN 5060 IREM = 0 MAN 5070 NCOLS = MAXCOL MAN 5080 NPAGE = 2 MAN 5090 GO TO 90 MAN 5100 80 CONTINUE MAN 5110 C MAN 5120 C HOW MANY ELEMENTS IN LAST ROW MAN 5130 C MAN 5140 IREM = MOD(NOP,NCOLS) MAN 5150 IF (NOP.LT.MAXCOL) IREM = 0 MAN 5160 C MAN 5170 C HOW MANY ROWS MAN 5180 C MAN 5190 L = NOP/NCOLS MAN 5200 IF (IREM.NE.0) L = L + 1 MAN 5210 90 CONTINUE MAN 5220 LM1 = L - 1 MAN 5230 C MAN 5240 C TOP MARKS THE INDEX OF THE FIRST ENTRY IN EACH COLUMN MAN 5250 C MAN 5260 TOP(1) = INC + 1 MAN 5270 TOP(2) = TOP(1) + L MAN 5280 IF (NCOLS.LT.3) GO TO 110 MAN 5290 DO 100 I=3,NCOLS MAN 5300 TOP(I) = TOP(I-1) + L MAN 5310 IF ((IREM.LT.I-1) .AND. (IREM.NE.0)) TOP(I) = TOP(I) - 1 MAN 5320 100 CONTINUE MAN 5330 110 CONTINUE MAN 5340 IF (IREM.EQ.0) LM1 = L MAN 5350 IF (L.EQ.1) LM1 = 1 MAN 5360 C MAN 5370 C PRINT OUT FULL ROWS MAN 5380 C MAN 5390 DO 140 I=1,LM1 MAN 5400 DO 120 K=1,NCOLS MAN 5410 ILOP = TOP(K) - 1 + LOP MAN 5420 IOP = TOP(K) - 1 + OP MAN 5430 IROP = TOP(K) - 1 + ROP MAN 5440 TOPTMP(1,K) = IWORK(ILOP) MAN 5450 TOPTMP(2,K) = IWORK(IOP) MAN 5460 TOPTMP(3,K) = IWORK(IROP) MAN 5470 TOPSAV(K) = TOP(K) MAN 5480 120 CONTINUE MAN 5490 C MAN 5500 C WRITE A LINE OF THE PROGRAM TREE LISTING IN A PRETTY FORMAT. MAN 5510 C MAN 5520 CALL WRITEL(TOP, TOPTMP, NCOLS) MAN 5530 DO 130 K=1,NCOLS MAN 5540 TOP(K) = TOPSAV(K) + 1 MAN 5550 130 CONTINUE MAN 5560 140 CONTINUE MAN 5570 IF (IREM.EQ.0) GO TO 160 MAN 5580 DO 150 K=1,IREM MAN 5590 ILOP = TOP(K) - 1 + LOP MAN 5600 IOP = TOP(K) - 1 + OP MAN 5610 IROP = TOP(K) - 1 + ROP MAN 5620 TOPTMP(1,K) = IWORK(ILOP) MAN 5630 TOPTMP(2,K) = IWORK(IOP) MAN 5640 TOPTMP(3,K) = IWORK(IROP) MAN 5650 150 CONTINUE MAN 5660 C MAN 5670 C WRITE A LINE OF THE PROGRAM TREE LISTING IN A PRETTY FORMAT. MAN 5680 C MAN 5690 CALL WRITEL(TOP, TOPTMP, IREM) MAN 5700 160 CONTINUE MAN 5710 C MAN 5720 NDONE = NDONE + MRMC MAN 5730 IF (NDONE.GT.NOPSAV) GO TO 170 MAN 5740 C MAN 5750 C IF .GT.1 PAGE NEEDED, RESET VALUES AND PROCESS NEXT PAGE MAN 5760 C MAN 5770 NOP = NOP - MRMC MAN 5780 INC = INC + MRMC MAN 5790 WRITE (OUTUNT,99999) MAN 5800 GO TO 70 MAN 5810 C MAN 5820 170 CONTINUE MAN 5830 NOP = NOPSAV MAN 5840 C MAN 5850 C MAN 5860 C MAN 5870 VALUE = RWRKFR MAN 5880 LPRD = RWRKFR + NOP MAN 5890 RWRKFR = RWRKFR + 3*NOP MAN 5900 180 CONTINUE MAN 5910 C MAN 5920 C READ OUTPUTS MAN 5930 C MAN 5940 190 READ (INUNIT,99981) M MAN 5950 C MAN 5960 C M = NUMBER OF OUTPUTS MAN 5970 C MAN 5980 C MAN 5990 C IF M = 0, READ IN NEXT PARAMETER MAN 6000 C MAN 6010 IF (M.EQ.0) GO TO 230 MAN 6020 C MAN 6030 C IF M .GT. 0, CONTINUE, ELSE, PRINT ERROR MESSAGE AND STOP MAN 6040 C MAN 6050 IF (M.GT.0) GO TO 200 MAN 6060 C MAN 6070 WRITE (OUTUNT,99990) MAN 6080 STOP MAN 6090 C MAN 6100 200 ISTORE = ISTORE + M MAN 6110 RSTORE = RSTORE + M*NOP MAN 6120 C MAN 6130 C CHECKS FOR OVERFLOW DUE TO M MAN 6140 C IF OVERFLOW, READ IN DUMMY VARIABLES MAN 6150 C MAN 6160 IF (.NOT.((ISTORE.GT.IWRKSZ) .OR. (RSTORE.GT.RWRKSZ) .OR. ERROR)) MAN 6170 * GO TO 210 MAN 6180 ERROR = .TRUE. MAN 6190 READ (INUNIT,99986) (IW1,I=1,M) MAN 6200 GO TO 220 MAN 6210 210 CONTINUE MAN 6220 C MAN 6230 OUTNDX = IWRKFR MAN 6240 IWRKFR = IWRKFR + M MAN 6250 C MAN 6260 BBAR = RWRKFR MAN 6270 RWRKFR = RWRKFR + M*NOP MAN 6280 C MAN 6290 IEND = IWRKFR - 1 MAN 6300 READ (INUNIT,99986) (IWORK(I),I=OUTNDX,IEND) MAN 6310 220 CONTINUE MAN 6320 C MAN 6330 C READ CONSTANTS MAN 6340 C MAN 6350 230 READ (INUNIT,99981) NCONS MAN 6360 C MAN 6370 C NCONS = NUMBER OF CONSTANTS MAN 6380 C MAN 6390 C MAN 6400 C IF NCONS = 0, READ IN NEXT PARAMETER MAN 6410 C MAN 6420 CONS = RWRKFR MAN 6430 NCDIM = 1 MAN 6440 IF (NCONS.EQ.0) GO TO 270 MAN 6450 C MAN 6460 C IF NCONS .GT. 0, CONTINUE, ELSE, PRINT ERROR MESSAGE AND STOP MAN 6470 C MAN 6480 IF (NCONS.GT.0) GO TO 240 MAN 6490 WRITE (OUTUNT,99989) MAN 6500 STOP MAN 6510 C MAN 6520 240 CONTINUE MAN 6530 RSTORE = RSTORE + NCONS MAN 6540 C MAN 6550 C CHECKS FOR OVERFLOW DUE TO NCONS MAN 6560 C MAN 6570 IF ((RSTORE.LT.RWRKSZ) .AND. (.NOT.ERROR)) GO TO 250 MAN 6580 ERROR = .TRUE. MAN 6590 READ (INUNIT,99985) (W1,I=1,NCONS) MAN 6600 GO TO 270 MAN 6610 250 CONTINUE MAN 6620 C MAN 6630 NCDIM = NCONS MAN 6640 RWRKFR = RWRKFR + NCONS MAN 6650 IEND = RWRKFR - 1 MAN 6660 READ (INUNIT,99985) (RWORK(I),I=CONS,IEND) MAN 6670 WRITE (OUTUNT,99980) MAN 6680 DO 260 K=1,NCONS MAN 6690 I = CONS + K - 1 MAN 6700 WRITE (OUTUNT,99982) K, RWORK(I) MAN 6710 260 CONTINUE MAN 6720 270 CONTINUE MAN 6730 C MAN 6740 FIRST = .TRUE. MAN 6750 C MAN 6760 C READ DATA ITEMS MAN 6770 C MAN 6780 280 CONTINUE MAN 6790 C MAN 6800 READ (INUNIT,99992) N, (NAME(I),I=1,77) MAN 6810 C MAN 6820 C N = NUMBER OF DATA ITEMS MAN 6830 C MAN 6840 C PRINTOUT DATA SET TITLE MAN 6850 C MAN 6860 IF ((.NOT.FIRST) .OR. (NCONS.EQ.0)) WRITE (OUTUNT,99999) MAN 6870 C MAN 6880 WRITE (OUTUNT,99998) (NAME(I),I=1,77) MAN 6890 C MAN 6900 C IF N = 0, DONE READING PARAMETERS MAN 6910 C MAN 6920 IF (N.EQ.0) GO TO 350 MAN 6930 C MAN 6940 IF (N.GT.0) GO TO 290 MAN 6950 IF ((N.LT.0) .AND. (.NOT.FIRST)) GO TO 30 MAN 6960 WRITE (OUTUNT,99995) MAN 6970 STOP MAN 6980 C MAN 6990 290 CONTINUE MAN 7000 IF (FIRST) NFIRST = N MAN 7010 IF (N.EQ.NFIRST) GO TO 300 MAN 7020 WRITE (OUTUNT,99995) MAN 7030 STOP MAN 7040 C MAN 7050 300 CONTINUE MAN 7060 IF (.NOT.FIRST) GO TO 310 MAN 7070 C MAN 7080 C ACTIVE DIMENSIONED TO MAX (NOP, N+M) MAN 7090 C MAN 7100 MAX = AMAX0(NOP,N+M) MAN 7110 C MAN 7120 C MAN 7130 ISTORE = ISTORE + MAX MAN 7140 RSTORE = RSTORE + N + M*(N+M) MAN 7150 C MAN 7160 C MAN 7170 C CHECKS FOR OVERFLOW DUE TO MAX MAN 7180 C IF OVERFLOW, READ IN DUMMY VARIABLES MAN 7190 C MAN 7200 310 CONTINUE MAN 7210 IF ((RSTORE.LE.RWRKSZ) .AND. (ISTORE.LE.IWRKSZ) .AND. MAN 7220 * (.NOT.ERROR)) GO TO 320 MAN 7230 ERROR = .TRUE. MAN 7240 READ (INUNIT,99985) (W1,I=1,N) MAN 7250 GO TO 350 MAN 7260 320 CONTINUE MAN 7270 C MAN 7280 IF (.NOT.FIRST) GO TO 330 MAN 7290 C MAN 7300 ACTIVE = IWRKFR MAN 7310 C MAN 7320 C MAN 7330 IWRKFR = ACTIVE + MAX MAN 7340 C MAN 7350 DATA = RWRKFR MAN 7360 C MAN 7370 ABAR = DATA + N MAN 7380 C MAN 7390 RWRKFR = ABAR + M*(N+M) MAN 7400 C MAN 7410 C MAN 7420 330 CONTINUE MAN 7430 C MAN 7440 C MAN 7450 IEND = DATA + N - 1 MAN 7460 READ (INUNIT,99985) (RWORK(I),I=DATA,IEND) MAN 7470 WRITE (OUTUNT,99983) MAN 7480 DO 340 K=1,N MAN 7490 I = DATA + K - 1 MAN 7500 WRITE (OUTUNT,99979) K, RWORK(I) MAN 7510 340 CONTINUE MAN 7520 350 CONTINUE MAN 7530 C MAN 7540 C IF WORK SPACE WAS EXCEEDED, PRINT AN ERROR MESSAGE AND CONTINUE MAN 7550 C WITH THE NEXT INPUT SET. MAN 7560 C MAN 7570 IF (ERROR) GO TO 360 MAN 7580 C MAN 7590 CALL BOUNDS(RWORK(ABAR), RWORK(BBAR), RWORK(LPRD), RWORK(VALUE), MAN 7600 * RWORK(CONS), IWORK(LOP), IWORK(OP), IWORK(ROP), RWORK(DATA), MAN 7610 * EPS, M, N, NOP, IWORK(OUTNDX), IWORK(ACTIVE), K2PM, K1, N+M, MAN 7620 * MAX, COMB, NCDIM) MAN 7630 C MAN 7640 GO TO 370 MAN 7650 C MAN 7660 360 WRITE (OUTUNT,99977) MAN 7670 WRITE (OUTUNT,99987) ISTORE, RSTORE MAN 7680 GO TO 460 MAN 7690 C MAN 7700 370 CONTINUE MAN 7710 C MAN 7720 C ALLOCATE ADDITIONAL SPACE AS REQUIRED. MAN 7730 C MAN 7740 IF (.NOT.FIRST) GO TO 380 MAN 7750 RSTORE = RSTORE + M*(M+2) MAN 7760 ERROR = RSTORE.GT.RWRKSZ MAN 7770 OUTA = RWRKFR MAN 7780 OUTB = OUTA + M MAN 7790 AATORR = OUTB + M MAN 7800 RWRKFR = AATORR + M**2 MAN 7810 380 CONTINUE MAN 7820 C MAN 7830 C CHOOSE METHOD USED FOR THE ANALYSIS. MAN 7840 C MAN 7850 IF (FLOAT((M**3))*COMB.LE.SWITCH) GO TO 420 MAN 7860 C MAN 7870 C MESSAGE----LARSON CHOSEN FOR ANALYSIS MAN 7880 C MAN 7890 WRITE (OUTUNT,99997) MAN 7900 C MAN 7910 C ADDITIONAL ALLOCATION FOR LARSON MAN 7920 C MAN 7930 IF (.NOT.FIRSTL) GO TO 390 MAN 7940 C MAN 7950 ISTORE = ISTORE + M + N + 1 MAN 7960 RSTORE = RSTORE + M*(2*M+N+8) + 4*N + 9 MAN 7970 C MAN 7980 390 CONTINUE MAN 7990 IF ((ISTORE.GT.IWRKSZ) .OR. (RSTORE.GT.RWRKSZ)) ERROR = .TRUE. MAN 8000 IF (.NOT.ERROR) GO TO 400 MAN 8010 WRITE (OUTUNT,99978) MAN 8020 GO TO 450 MAN 8030 C MAN 8040 400 CONTINUE MAN 8050 C MAN 8060 IF (.NOT.FIRST) GO TO 410 MAN 8070 SL1 = IWRKFR MAN 8080 C MAN 8090 IWRKFR = SL1 + M + N + 1 MAN 8100 C MAN 8110 XL1 = RWRKFR MAN 8120 C MAN 8130 BBT = XL1 + M + 1 MAN 8140 C MAN 8150 BL1 = BBT + M**2 MAN 8160 C MAN 8170 EL1 = BL1 + M + N + 1 MAN 8180 C MAN 8190 AL1 = EL1 + M + N + 1 MAN 8200 C MAN 8210 RWRKFR = AL1 + (M+2)*(M+N+3) MAN 8220 C MAN 8230 410 CONTINUE MAN 8240 FIRST = .FALSE. MAN 8250 CALL LARSON(IWORK(OUTNDX), IWORK(SL1), RWORK(ABAR), RWORK(BBAR), MAN 8260 * RWORK(AL1), RWORK(BL1), RWORK(AATORR), RWORK(BBT), RWORK(EL1), MAN 8270 * RWORK(OUTA), RWORK(OUTB), RWORK(XL1), RWORK(VALUE), NOP, M, N+M, MAN 8280 * N+M+1, N+M+3, M+1, M+2, EPS, K2PM, K1, RWORK(DATA), N) MAN 8290 C MAN 8300 FIRSTL = .FALSE. MAN 8310 C MAN 8320 GO TO 450 MAN 8330 C MAN 8340 420 CONTINUE MAN 8350 C MAN 8360 C MAN 8370 C MESSAGE----FACE CHOSEN FOR ANALYSIS MAN 8380 C MAN 8390 C MAN 8400 WRITE (OUTUNT,99996) MAN 8410 C MAN 8420 C MAN 8430 C ADDITIONAL STORAGE ALLOCATION FOR FACE MAN 8440 C MAN 8450 C MAN 8460 IF (FIRST) ISTORE = ISTORE + 3*M MAN 8470 IF (FIRST) RSTORE = RSTORE + 3*M MAN 8480 C MAN 8490 ERROR = ISTORE.GT.IWRKSZ .OR. RSTORE.GT.RWRKSZ MAN 8500 IF (.NOT.ERROR) GO TO 430 MAN 8510 WRITE (OUTUNT,99978) MAN 8520 GO TO 450 MAN 8530 C MAN 8540 430 CONTINUE MAN 8550 C MAN 8560 C MAN 8570 IF (.NOT.FIRSTF) GO TO 440 MAN 8580 COL = IWRKFR MAN 8590 C MAN 8600 LIMIT = COL + M MAN 8610 C MAN 8620 JPVT = LIMIT + M MAN 8630 C MAN 8640 IWRKFR = JPVT + M MAN 8650 C MAN 8660 OUTORQ = RWRKFR MAN 8670 C MAN 8680 V = OUTORQ + M MAN 8690 C MAN 8700 WORK = V + M MAN 8710 C MAN 8720 RWRKFR = WORK + M MAN 8730 C MAN 8740 440 CONTINUE MAN 8750 FIRST = .FALSE. MAN 8760 NCOMB = COMB + HALF MAN 8770 CALL FACE(RWORK(ABAR), RWORK(BBAR), RWORK(OUTA), RWORK(OUTB), MAN 8780 * RWORK(OUTORQ), RWORK(AATORR), RWORK(V), RWORK(VALUE), MAN 8790 * RWORK(WORK), IWORK(COL), IWORK(JPVT), IWORK(LIMIT), MAN 8800 * IWORK(OUTNDX), NOP, M, N+M, EPS, K2PM, K1, NCOMB, RWORK(DATA), N)MAN 8810 C MAN 8820 FIRSTF = .FALSE. MAN 8830 C MAN 8840 450 CONTINUE MAN 8850 WRITE (OUTUNT,99987) ISTORE, RSTORE MAN 8860 460 IF (IFLAG.NE.0) GO TO 280 MAN 8870 GO TO 30 MAN 8880 C MAN 8890 C FORMAT STATEMENTS MAN 8900 C MAN 8910 99999 FORMAT (1H1) MAN 8920 99998 FORMAT (1H , 77A1) MAN 8930 99997 FORMAT (21H-THE HEURISTIC METHOD, 27H WILL PERFORM THE ANALYSIS.) MAN 8940 99996 FORMAT (22H-THE EXACT METHOD WILL, 22H PERFORM THE ANALYSIS.) MAN 8950 99995 FORMAT (31H-ERROR IN NUMBER OF DATA ITEMS./19H NOT ALL SETS OF DA,MAN 8960 * 12HTA AGREED ON, 25HTHE NUMBER OF DATA ITEMS./15H RECHECK INPUT.)MAN 8970 99994 FORMAT (/13H PROGRAM TREE) MAN 8980 99993 FORMAT (/49H V(I) IS THE VALUE OF THE ITH COMPUTATIONAL STEP./ MAN 8990 * 40H D(I) IS THE VALUE OF THE ITH DATA ITEM./17H C(I) IS THE VALU,MAN 9000 * 22HE OF THE ITH CONSTANT.) MAN 9010 99992 FORMAT (I2, 77A1) MAN 9020 99991 FORMAT (2I3, 74A1) MAN 9030 99990 FORMAT (1X, 34HERROR- NEGATIVE NUMBER OF OUTPUTS/12HRECHECK INPU,MAN 9040 * 22HT---PROGRAM TERMINATED) MAN 9050 99989 FORMAT (1X, 36HERROR- NEGATIVE NUMBER OF CONSTANTS/10HRECHECK IN,MAN 9060 * 24HPUT---PROGRAM TERMINATED) MAN 9070 99988 FORMAT (1H1, 74A1) MAN 9080 99987 FORMAT (19H-STORAGE REQUIRED -/9H INTEGER , I6/9H REAL , I6) MAN 9090 99986 FORMAT (I3) MAN 9100 99985 FORMAT (E22.10) MAN 9110 99984 FORMAT (I3, I2, I4) MAN 9120 99983 FORMAT (5H-DATA) MAN 9130 99982 FORMAT (3H C(, I2, 3H) =, 1PE15.7) MAN 9140 99981 FORMAT (I2) MAN 9150 99980 FORMAT (10H1CONSTANTS) MAN 9160 99979 FORMAT (3H D(, I2, 3H) =, 1PE15.7) MAN 9170 99978 FORMAT (6H- , 30(1H*)/25H0ERROR - STORAGE OVERFLOW) MAN 9180 99977 FORMAT (42H-NOT ENOUGH STORAGE FOR SUBROUTINE BOUNDS.) MAN 9190 END MAN 9200 SUBROUTINE WRITEL(TOP, TOPTMP, NCOLS) WRI 10 C C OUTPUTS A SINGLE LINE OF THE PROGRAM TREE LISTING IN A PRETTY C FORMAT. THIS IS A MULTIPLE COLUMN VERSION OF MILLER'S FINISH C ROUTINE. CONSEQUENTLY, IT CONSERVES ON PAPER, AND IS MUCH MORE C READABLE THAN A SINGLE COLUMN LISTING FOR LONG PROGRAM TREES. C THE CODE IS LIMITED TO A MAXIMUM OF 999 PROGRAM STEPS. C MODIFICATION IS REQUIRED FOR LONGER PROGRAMS. C C JOHN A. WISNIEWSKI SANDIA NATIONAL LABORATORIES C JANUARY 6,1982 ALBUQUERQUE, N.M. C C INTEGER TOP(NCOLS), TOPTMP(3,NCOLS), LINE(131), RPAREN, LPAREN, * NCOLS INTEGER STAR, PLUS, MINUS, SLASH, S, Q, R, T, IBLANK, V, D INTEGER DIGIT(10), EQUALS INTEGER INUNIT, OUTUNT C COMMON /IOINFO/ INUNIT, OUTUNT C DATA RPAREN, LPAREN, STAR, PLUS, MINUS, SLASH /1H),1H(,1H*,1H+, * 1H-,1H// DATA S, Q, R, T, IBLANK, V, D, EQUALS /1HS,1HQ,1HR,1HT,1H ,1HV, C * 1HD,1H=/ DATA C /1HC/ DATA DIGIT(1), DIGIT(2), DIGIT(3), DIGIT(4), DIGIT(5), DIGIT(6), * DIGIT(7), DIGIT(8), DIGIT(9), DIGIT(10) /1H0,1H1,1H2,1H3,1H4,1H5, * 1H6,1H7,1H8,1H9/ C DO 10 I=1,131 LINE(I) = IBLANK 10 CONTINUE C DO 90 I=1,NCOLS C C SET UP FIELD TO THE LEFT OF THE EQUALS SIGN. C IOFF = 26*(I-1) + 1 LINE(IOFF+1) = V LINE(IOFF+2) = LPAREN LINE(IOFF+6) = RPAREN LINE(IOFF+8) = EQUALS IF (TOP(I).LE.0) GO TO 20 J = MOD(TOP(I),10) + 1 LINE(IOFF+5) = DIGIT(J) TOP(I) = TOP(I)/10 IF (TOP(I).LE.0) GO TO 20 J = MOD(TOP(I),10) + 1 LINE(IOFF+4) = DIGIT(J) TOP(I) = TOP(I)/10 IF (TOP(I).LE.0) GO TO 20 J = MOD(TOP(I),10) + 1 LINE(IOFF+3) = DIGIT(J) 20 CONTINUE IF (TOPTMP(2,I).NE.5) GO TO 40 C C UNARY OPERATION SQUARE ROOT. C LINE(IOFF+10) = S LINE(IOFF+11) = Q LINE(IOFF+12) = R LINE(IOFF+13) = T LINE(IOFF+14) = LPAREN LINE(IOFF+15) = V IF (TOPTMP(1,I).LE.100) LINE(IOFF+15) = D IF (TOPTMP(1,I).LT.0) LINE(IOFF+15) = C IF (TOPTMP(1,I).GT.100) TOPTMP(1,I) = TOPTMP(1,I) - 100 TOPTMP(1,I) = IABS(TOPTMP(1,I)) LINE(IOFF+16) = LPAREN IF (TOPTMP(1,I).LE.0) GO TO 30 J = MOD(TOPTMP(1,I),10) + 1 LINE(IOFF+19) = DIGIT(J) TOPTMP(1,I) = TOPTMP(1,I)/10 IF (TOPTMP(1,I).LE.0) GO TO 30 J = MOD(TOPTMP(1,I),10) + 1 LINE(IOFF+18) = DIGIT(J) TOPTMP(1,I) = TOPTMP(1,I)/10 IF (TOPTMP(1,I).LE.0) GO TO 30 J = MOD(TOPTMP(1,I),10) + 1 LINE(IOFF+17) = DIGIT(J) 30 CONTINUE LINE(IOFF+20) = RPAREN LINE(IOFF+21) = RPAREN GO TO 90 40 CONTINUE IF (TOPTMP(2,I).NE.6) GO TO 60 C C UNARY OPERATION NEGATION. C LINE(IOFF+10) = MINUS LINE(IOFF+11) = V IF (TOPTMP(1,I).LE.100) LINE(IOFF+11) = D IF (TOPTMP(1,I).LT.0) LINE(IOFF+11) = C IF (TOPTMP(1,I).GT.100) TOPTMP(1,I) = TOPTMP(1,I) - 100 TOPTMP(1,I) = IABS(TOPTMP(1,I)) LINE(IOFF+12) = LPAREN IF (TOPTMP(1,I).LE.0) GO TO 50 J = MOD(TOPTMP(1,I),10) + 1 LINE(IOFF+15) = DIGIT(J) TOPTMP(1,I) = TOPTMP(1,I)/10 IF (TOPTMP(1,I).LE.0) GO TO 50 J = MOD(TOPTMP(1,I),10) + 1 LINE(IOFF+14) = DIGIT(J) TOPTMP(1,I) = TOPTMP(1,I)/10 IF (TOPTMP(1,I).LE.0) GO TO 50 J = MOD(TOPTMP(1,I),10) + 1 LINE(IOFF+13) = DIGIT(J) 50 CONTINUE LINE(IOFF+16) = RPAREN GO TO 90 60 CONTINUE C C BINARY OPERATIONS. C LINE(IOFF+10) = V IF (TOPTMP(1,I).LE.100) LINE(IOFF+10) = D IF (TOPTMP(1,I).LT.0) LINE(IOFF+10) = C IF (TOPTMP(1,I).GT.100) TOPTMP(1,I) = TOPTMP(1,I) - 100 TOPTMP(1,I) = IABS(TOPTMP(1,I)) LINE(IOFF+11) = LPAREN LINE(IOFF+15) = RPAREN IF (TOPTMP(1,I).LE.0) GO TO 70 J = MOD(TOPTMP(1,I),10) + 1 LINE(IOFF+14) = DIGIT(J) TOPTMP(1,I) = TOPTMP(1,I)/10 IF (TOPTMP(1,I).LE.0) GO TO 70 J = MOD(TOPTMP(1,I),10) + 1 LINE(IOFF+13) = DIGIT(J) TOPTMP(1,I) = TOPTMP(1,I)/10 IF (TOPTMP(1,I).LE.0) GO TO 70 J = MOD(TOPTMP(1,I),10) + 1 LINE(IOFF+13) = DIGIT(J) 70 CONTINUE IF (TOPTMP(2,I).EQ.1) LINE(IOFF+16) = PLUS IF (TOPTMP(2,I).EQ.2) LINE(IOFF+16) = MINUS IF (TOPTMP(2,I).EQ.3) LINE(IOFF+16) = STAR IF (TOPTMP(2,I).EQ.4) LINE(IOFF+16) = SLASH C LINE(IOFF+17) = V IF (TOPTMP(3,I).LE.100) LINE(IOFF+17) = D IF (TOPTMP(3,I).LT.0) LINE(IOFF+17) = C IF (TOPTMP(3,I).GT.100) TOPTMP(3,I) = TOPTMP(3,I) - 100 TOPTMP(3,I) = IABS(TOPTMP(3,I)) LINE(IOFF+18) = LPAREN LINE(IOFF+22) = RPAREN IF (TOPTMP(3,I).LE.0) GO TO 80 J = MOD(TOPTMP(3,I),10) + 1 LINE(IOFF+21) = DIGIT(J) TOPTMP(3,I) = TOPTMP(3,I)/10 IF (TOPTMP(3,I).LE.0) GO TO 80 J = MOD(TOPTMP(3,I),10) + 1 LINE(IOFF+20) = DIGIT(J) TOPTMP(3,I) = TOPTMP(3,I)/10 IF (TOPTMP(3,I).LE.0) GO TO 80 J = MOD(TOPTMP(3,I),10) + 1 LINE(IOFF+19) = DIGIT(J) 80 CONTINUE 90 CONTINUE C WRITE (OUTUNT,99999) (LINE(I),I=1,131) 99999 FORMAT (131A1) RETURN END SUBROUTINE BOUNDS(ABAR, BBAR, LPRD, VALUE, CONS, LOP, OP, ROP, BOU 10 * DATA, EPS, M, N, NOP, OUTNDX, ACTIVE, K2PM, K1, NPM, MAX, COMB, * NCDIM) C REAL ABAR(M,NPM), BBAR(M,NOP), CONS(NCDIM), DATA(N), LPRD(NOP,2), * VALUE(NOP) C INTEGER ACTIVE(MAX), OUTNDX(M), OP(NOP), ROP(NOP), LOP(NOP) C REAL ONE(1), ZERO(1) C C C NPM = N + M C NPMP1 = N + M + 1 C NPMP3 = N + M + 3 C MP1 = M + 1 C MP2 = M + 2 C C NOP = NUMBER OF OPERATIONS C NCDIM = NUMBER OF CONSTANTS C M = NUMBER OF OUTPUTS C N = NUMBER OF DATA C INTEGER OJ C INTEGER INUNT, OUTUNT COMMON /IOINFO/ INUNT, OUTUNT C DATA ONE(1), ZERO(1), HALF /1.0E0,0.0E0,0.5E0/ C C INITIALIZE ABAR=ADOT, BBAR C DO 10 I=1,M CALL SCOPY(NPM, ZERO(1), 0, ABAR(I,1), M) CALL SCOPY(NOP, ZERO(1), 0, BBAR(I,1), M) 10 CONTINUE C C C C C THIS LOOP COMPUTES THE NODE VALUES AND THE ARC WEIGHTS. C LPRD(K,1) IS THE WEIGHT ON THE LEFT ARC. C LPRD(K,2) IS THE WEIGHT ON THE RIGHT ARC. C DO 180 K=1,NOP ILOP = LOP(K) IOP = OP(K) IROP = ROP(K) IF (ILOP.LT.0) GO TO 30 IF (ILOP.LE.100) GO TO 20 ARGL = VALUE(ILOP-100) GO TO 40 20 ARGL = DATA(ILOP) GO TO 40 30 CONTINUE ILOP = -ILOP ARGL = CONS(ILOP) ILOP = -ILOP 40 IF (IROP.LT.0) GO TO 60 IF (IROP.EQ.0) GO TO 70 IF (IROP.LE.100) GO TO 50 ARGR = VALUE(IROP-100) GO TO 70 50 ARGR = DATA(IROP) GO TO 70 60 CONTINUE IROP = -IROP ARGR = CONS(IROP) IROP = -IROP C C WHICH OPERATION--- C 70 GO TO (80, 100, 120, 130, 140, 150), IOP C C ADDITION OPERATION C 80 VALUE(K) = ARGL + ARGR IF (VALUE(K).NE.ZERO(1)) GO TO 90 C C IF VALUE(K) = 0, PERTURB THE PROBLEM BY REPLACING IT WITH A VALUE C CLOSE TO ZERO IN ORDER TO DIVIDE BY VALUE(K) C VALUE(K) = EPS*AMAX1(ARGL,ARGR) IF ((ARGL.EQ.ZERO(1)) .AND. (ARGR.EQ.ZERO(1))) VALUE(K) = EPS 90 LPRD(K,1) = ARGL/VALUE(K) LPRD(K,2) = ARGR/VALUE(K) GO TO 180 C C SUBTRACTION OPERATION C 100 VALUE(K) = ARGL - ARGR IF (VALUE(K).NE.ZERO(1)) GO TO 110 C C IF VALUE(K) = 0, PERTURB THE PROBLEM BY REPLACING IT WITH A VALUE C CLOSE TO ZERO IN ORDER TO DIVIDE BY VALUE(K) C VALUE(K) = EPS*AMAX1(ARGL,ARGR) IF ((ARGL.EQ.ZERO(1)) .AND. (ARGR.EQ.ZERO(1))) VALUE(K) = EPS 110 LPRD(K,1) = ARGL/VALUE(K) LPRD(K,2) = -ARGR/VALUE(K) GO TO 180 C C MULTIPLICATION OPERATION C 120 VALUE(K) = ARGL*ARGR LPRD(K,1) = ONE(1) LPRD(K,2) = ONE(1) GO TO 180 C C DIVISION OPERATION C 130 IF (ARGR.EQ.ZERO(1)) GO TO 160 VALUE(K) = ARGL/ARGR LPRD(K,1) = ONE(1) LPRD(K,2) = -ONE(1) GO TO 180 C C SQUARE ROOT OPERATION C 140 IF (ARGL.LT.ZERO(1)) GO TO 170 VALUE(K) = SQRT(ARGL) LPRD(K,1) = HALF LPRD(K,2) = ZERO(1) GO TO 180 C C NEGATION OPERATION C 150 VALUE(K) = -ARGL LPRD(K,1) = ONE(1) LPRD(K,2) = ZERO(1) GO TO 180 160 WRITE (OUTUNT,99998) K STOP 170 WRITE (OUTUNT,99999) K STOP 180 CONTINUE C C THIS LOOP COMPUTES ABAR AND BBAR FROM THE ARC WEIGHTS. C THIS PARTICULAR IMPLEMENTATION CORRESPONDS TO (4.1) IN TOMS. C DO 210 J=1,M OJ = OUTNDX(J) BBAR(J,OJ) = ONE(1) DO 200 II=1,OJ I = OJ - II + 1 ILOP = LOP(I) IF (ILOP.LE.0) GO TO 190 IF (ILOP.LE.100) ABAR(J,ILOP) = ABAR(J,ILOP) + * LPRD(I,1)*BBAR(J,I) IF (ILOP.GT.100) BBAR(J,ILOP-100) = BBAR(J,ILOP-100) + * LPRD(I,1)*BBAR(J,I) 190 IROP = ROP(I) IF (IROP.LE.0) GO TO 200 IF (IROP.LE.100) ABAR(J,IROP) = ABAR(J,IROP) + * LPRD(I,2)*BBAR(J,I) IF (IROP.GT.100) BBAR(J,IROP-100) = BBAR(J,IROP-100) + * LPRD(I,2)*BBAR(J,I) 200 CONTINUE 210 CONTINUE C C ABAR BECOMES ADOT. C CALL SCOPY(M, ONE(1), 0, ABAR(1,N+1), M+1) C C CONDENSE ADOT AND BBAR. C CALL CONDNS(ABAR(1,1), M, M, NPM, K2PM, ACTIVE(1), EPS) CALL CONDNS(BBAR(1,1), M, M, NOP, K1, ACTIVE(1), EPS) C C MM1 = M - 1 COMB = RCOMB(K2PM,MM1) C RETURN C 99999 FORMAT (40H-PROGRAM TRIES TO TAKE SQUARE ROOT OF NE, 9HGATIVE NU, * 5HMBER./33H ERROR OCCURS AT OPERATION NUMBER, I3/12H RECHECK INP, * 3HUT.) 99998 FORMAT (33H-PROGRAM TRIES TO DIVIDE BY ZERO//17H ERROR OCCURS AT , * 16HOPERATION NUMBER, I3/15H RECHECK INPUT.) C END SUBROUTINE LARSON(OUTNDX, SL1, ABAR, BBAR, AL1, BL1, AAT, BBT, LAR 10 * EL1, OUTA, OUTB, XL1, VALUE, NOP, M, NPM, NPMP1, NPMP3, MP1, * MP2, EPS, K2PM, K1, DATA, NDATA) C C LARSON SUBROUTINE C C C PERFORMS THE ERROR ANALYSIS USING AN HEURISTIC METHOD C C THIS SUBROUTINE CONTAINS MOST OF THE CODE FROM THE LARSON C MAIN PROGRAM. THE LARSON PROGRAM HAS BEEN CONVERTED TO C A MAIN PROGRAM WHICH READS IN THE DATA AND ALLOCATES ALL C INPUT DEPENDENT ARRAY STORAGE IN WORK ARRAYS AND PASSES C THE ARRAYS TO THIS SUBROUTINE. C C C EXTERNALS--ISAMAX,SASUM,SAXPY,SCOPY,SDOT. C C INTEGER ARRAYS C INTEGER OUTNDX(M), SL1(NPMP1) C C C REAL ARRAYS C REAL ABAR(M,NPM), BBAR(M,NOP), AL1(NPMP3,MP2), BL1(NPMP1), * AAT(M,M), BBT(M,M), EL1(NPMP1), OUTA(M), OUTB(M), XL1(MP1), * VALUE(NOP), DATA(NDATA) C REAL LAMBDA, ZERO(1) C C INTEGER INUNT, OUTUNT COMMON /IOINFO/ INUNT, OUTUNT C C DATA PTONE, ONE, TWO, ZERO(1) /0.1E0,1.0E0,2.0E0,0.0E0/ C C IF (M.EQ.1) GO TO 90 C C COMPUTE ADOT*ADOT**T, AND BBAR*BBAR**T C DO 20 I=1,M DO 10 J=1,I BBT(I,J) = SDOT(K1,BBAR(I,1),M,BBAR(J,1),M) BBT(J,I) = BBT(I,J) AAT(J,I) = SDOT(K2PM,ABAR(I,1),M,ABAR(J,1),M) AAT(I,J) = AAT(J,I) 10 CONTINUE 20 CONTINUE C C USE THE QZ ALGORITHM TO COMPUTE THE EIGENVALUES AND C CORRESPONDING EIGENVECTORS OF THE NORMAL EQUATIONS FORMULATION C OF THE GENERALIZED SVD BBAR*X = SIGMA*ADOT*X. C CALL QZHES(M, M, BBT(1,1), AAT(1,1), .TRUE., AL1(1,1)) CALL QZIT(M, M, BBT(1,1), AAT(1,1), EPS, .TRUE., AL1(1,1), INFO) CALL QZVAL(M, M, BBT(1,1), AAT(1,1), EL1(1), XL1(1), BL1(1), * .TRUE., AL1(1,1)) IF (INFO.EQ.0) GO TO 30 WRITE (OUTUNT,99994) RETURN 30 CONTINUE CALL QZVEC(M, M, BBT(1,1), AAT(1,1), EL1(1), XL1(1), BL1(1), * AL1(1,1)) C C FIND THE LARGEST AND SMALLEST ABSOLUTE ELEMENTS OF QZ-BETA. C BMIN = ABS(BL1(1)) BMAX = BMIN DO 40 K=2,M BTEST = ABS(BL1(K)) IF (BTEST.LT.BMIN) BMIN = BTEST IF (BTEST.GT.BMAX) BMAX = BTEST 40 CONTINUE C C CHECK FOR A NUMERICALLY SINGULAR AAT MATRIX. IF SO, THEN THETA C IS INFINITE TO MACHINE PRECISION. C IF (PTONE*BMIN.LT.EPS*BMAX) GO TO 90 C C COMPUTE LARGEST EIGENVALUE. C LAMBDA = ZERO(1) DO 50 K=1,M RATIO = EL1(K)/BL1(K) IF (RATIO.LT.LAMBDA) GO TO 50 LAMBDA = RATIO I = K 50 CONTINUE CALL SCOPY(M, ZERO(1), 0, XL1(1), 1) C C CHOOSE GAMMA(I) SO THAT BBAR(I) HAS POSITIVE PROJECTION ON C EIGENVECTOR. C DO 60 J=1,K1 SUM = SDOT(M,AL1(1,I),1,BBAR(1,J),1) SIGN = ONE IF (SUM.EQ.ZERO(1)) WRITE (OUTUNT,99995) J IF (SUM.LT.ZERO(1)) SIGN = -ONE CALL SAXPY(M, SIGN, BBAR(1,J), 1, XL1(1), 1) 60 CONTINUE C C COMPUTE SCALING FACTOR FOR L1 LINEAR PROGRAM. C C C FIND THE INFINITY NORM OF MATRIX ABAR C ANRM = ZERO(1) DO 70 I=1,M SUM = SASUM(K2PM,ABAR(I,1),M) IF (SUM.GT.ANRM) ANRM = SUM 70 CONTINUE C I = ISAMAX(M,XL1(1),1) C C FIND THE INFINITY NORM OF VECTOR XL1. C OUTNRM = ABS(XL1(I)) C KL1 = TWO*ANRM/OUTNRM + ONE K2PMP1 = K2PM + 1 C C INITIALIZE LINEAR PROGRAM ARRAY AND RIGHT HAND SIDE. C DO 80 J=1,M CALL SCOPY(K2PM, ABAR(J,1), M, AL1(1,J), 1) AL1(K2PMP1,J) = -FLOAT(KL1)*XL1(J) 80 CONTINUE CALL SCOPY(K2PM, ZERO(1), 0, BL1(1), 1) BL1(K2PMP1) = ONE TOLER = SQRT(SQRT(EPS**3)) C C ASCHER VERSION OF BARRODALE AND ROBERTS L1 LINEAR PROGRAM. C CALL L1(K2PMP1, M, NPMP3, MP2, AL1(1,1), BL1(1), TOLER, XL1(1), * EL1(1), SL1(1)) C C RECOVER SOLUTION TO ORIGINAL PROBLEM. C TEMP = AL1(K2PM+2,M+1) IF (TEMP.NE.ZERO(1)) THETA = ONE/FLOAT(KL1)/TEMP 90 CONTINUE C C FORWARD ANALYSIS, CONDITION OF PROBLEM AND ALGORITHM. C DO 100 I=1,M OUTA(I) = SASUM(K2PM,ABAR(I,1),M) - ONE OUTB(I) = SASUM(K1,BBAR(I,1),M) 100 CONTINUE WRITE (OUTUNT,99997) DO 120 I=1,M C C SUBTRACT IDENTITY FROM ADOT TO GET ABAR FOR FORWARD C ANALYSIS. C J = OUTNDX(I) IF (J.LT.0) GO TO 110 TOTOUT = OUTA(I) + OUTB(I) WRITE (OUTUNT,99998) J, VALUE(J), OUTA(I), OUTB(I), TOTOUT GO TO 120 110 CONTINUE MINUSJ = -J WRITE (OUTUNT,99998) J, DATA(MINUSJ) 120 CONTINUE IF (M.EQ.1) GO TO 140 IF (PTONE*BMIN.LT.EPS*BMAX) TEMP = ZERO(1) C IF (TEMP.NE.ZERO(1)) GO TO 130 WRITE (OUTUNT,99999) RETURN 130 CONTINUE C WRITE (OUTUNT,99996) THETA GO TO 160 140 CONTINUE C C SPECIAL CASE WHEN ONLY ONE OUTPUT. C IF (OUTA(1).NE.(-ONE)) GO TO 150 WRITE (OUTUNT,99999) RETURN C 150 CONTINUE THETA = OUTB(1)/(OUTA(1)+ONE) WRITE (OUTUNT,99996) THETA 160 CONTINUE RLAMDA = SQRT(LAMBDA) BNDL = THETA BNDR = RLAMDA*SQRT(FLOAT(K1)) WRITE (OUTUNT,99993) BNDL, BNDR RETURN C 99999 FORMAT (18H-THETA IS INFINITE) 99998 FORMAT (1H , I6, 4X, 1PE20.7, 4X, 1PE20.7, 13X, 1PE20.7, 12X, * 1PE20.7) 99997 FORMAT (12H-OUTPUT NODE, 10X, 5HVALUE, 12X, 18HCONDITION OF PROBL, * 2HEM, 12X, 22HCONDITION OF ALGORITHM, 11X, 18HRELATIVE ERROR BOU, * 2HND) 99996 FORMAT (9H-THETA = , 1PE13.5) 99995 FORMAT (23H-CHOOSING +1 FOR GAMMA(, I3, 1H)) 99994 FORMAT (35H-ABNORMAL TERMINATION FROM EISPACK./15H CONTACT AUTHOR, * 26HS CONCERNING THIS PROBLEM.) 99993 FORMAT (25H-INTERVAL ABOUT THETA = (, 1PE15.5, 2H ,, 1PE15.5, * 2H )) END SUBROUTINE FACE(ABAR, BBAR, OUTA, OUTB, QRAUX, R, V, VALUE, WORK, FAC 10 * COL, JPVT, LIMIT, OUTNDX, NOP, M, NPM, EPS, K2PM, K1, NCOMB, * DATA, NDATA) C C C C PERFORMS THE ERROR ANALYSIS USING THE EXACT METHOD C C THIS SUBROUTINE IS THE CODE FROM THE MAIN PROGRAM OF THE C FACE SUBROUTINE PACKAGE EXCEPT FOR THE INPUT OF DATA VALUES. C THIS SUBROUTINE IS CALLED BY A MAIN PROGRAM WHICH INPUTS THE C DATA VALUES AND ALLOCATES ALL ARRAY STORAGE USED BY THIS C SUBROUTINE. C C C EXTERNALS NEEDED = C BLAS - SASUM C SCOPY C SDOT C LINPACK - SQRDC C SQRSL C C INTEGER NOP, M C REAL ABAR(M,NPM), BBAR(M,NOP), OUTA(M), OUTB(M), QRAUX(M), * R(M,M), V(M), VALUE(NOP), WORK(M), DATA(NDATA) C INTEGER COL(M), JPVT(M), LIMIT(M), OUTNDX(M) C REAL ZERO(1) C INTEGER INUNT, OUTUNT COMMON /IOINFO/ INUNT, OUTUNT C DATA ONE, ZERO(1), FIFTN /1.0E0,0.0E0,15.0E0/ C MM1 = M - 1 C C C NCOMB IS THE POSSIBLE NUMBER OF FACE PAIRS FOR THE POLYTOPE C OF THE PROBLEM. C THETA = ZERO(1) C C THE 9999 LOOP INVESTIGATES EACH POSSIBLE FACE FOR B-ANAYSIS. C C DO 80 ITER=1,NCOMB ITR = ITER C C SELECT ONE COMBINATION OF M-1 COLUMNS OF ADOT. C LDIM = MAX0(1,MM1) CALL COMBIN(COL, LIMIT, K2PM, MM1, ITR, LDIM) C C DO 10 J=1,MM1 JJ = COL(J) CALL SCOPY(M, ABAR(1,JJ), 1, R(1,J), 1) 10 CONTINUE C C PERFORM A QR DECOMPOSITION. C CALL SQRDC(R(1,1), M, M, MM1, QRAUX(1), JPVT(1), WORK(1), 1) C C CHECK IF THESE M-1 COLUMNS ARE LINEARLY DEPENDENT. C RMAX = ABS(R(1,1)) C RMIN = RMAX IF (MM1.GT.0) RMIN = ABS(R(MM1,MM1)) C C C TOL IS A CONSERVATIVE ESTIMATE FOR THE TOLERANCE GIVEN IN C WILKINSON, ALGEBRAIC EIGENVALUE PROBLEM , CH.3, SEC.45 C TOL = FIFTN*FLOAT(M)*EPS C C IF (RMIN.LE.(TOL*RMAX)) GO TO 80 C C FIND THE LAST COLUMN OF Q. THIS VECTOR, V, IS ORTHOGONAL TO C THE FACE. C CALL SCOPY(MM1, ZERO(1), 0, V(1), 1) V(M) = ONE C C NOW, FORM Q * V C CALL SQRSL(R(1,1), M, M, MM1, QRAUX(1), V(1), V(1), WORK(1), * WORK(1), WORK(1), WORK(1), 10000, INFO) IF (INFO.EQ.0) GO TO 20 WRITE (OUTUNT,99999) STOP C C C C GET RID OF NOISE. C 20 CONTINUE C DO 30 I=1,M IF (ABS(V(I)).LT.TOL) V(I) = ZERO(1) 30 CONTINUE OMEGA = ZERO(1) C C FIND OMEGA, THE PROJECTION OF THE FACE ON THE NORMAL, V. C COL(M) = K2PM + 1 I = 1 C DO 60 J=1,K2PM C C COLUMNS WHICH ARE ORTHOGONAL TO THE FACE HAVE ZERO C PROJECTIONS. SKIP. C 40 IF (J.EQ.COL(I)) GO TO 60 IF (J.LT.COL(I)) GO TO 50 I = I + 1 GO TO 40 50 OMEGA = OMEGA + ABS(SDOT(M,V(1),1,ABAR(1,J),1)) C 60 CONTINUE C IF (OMEGA.EQ.ZERO(1)) GO TO 90 C THETA1 = ZERO(1) C C FIND THETA1, THE LARGEST PROJECTION ON V OF A VERTEX OF THE C POLYTOPE OF THE ALGORITHM. C DO 70 J=1,K1 THETA1 = THETA1 + ABS(SDOT(M,V(1),1,BBAR(1,J),1)) 70 CONTINUE C C THETA = DISTANCE OF VERTEX(BBAR) / DISTANCE OF FACE. C THETA1 = THETA1/OMEGA C C KEEP TRACK OF LARGEST THETA THUS ENCOUNTERED. C IF (THETA1.GE.THETA) THETA = THETA1 C 80 CONTINUE C 90 CONTINUE C C C SUBTRACT IDENTITY FROM ADOT TO GET ABAR FOR FORWARD ANALYSIS. C DO 100 I=1,M C C FORWARD ANALYSIS, CONDITION OF THE PROBLEM AND ALGORITHM. C OUTA(I) = -ONE + SASUM(K2PM,ABAR(I,1),M) OUTB(I) = SASUM(K1,BBAR(I,1),M) 100 CONTINUE WRITE (OUTUNT,99996) C C FORWARD ERROR BOUND C DO 120 I=1,M J = OUTNDX(I) IF (J.LT.0) GO TO 110 TOTOUT = OUTA(I) + OUTB(I) WRITE (OUTUNT,99995) J, VALUE(J), OUTA(I), OUTB(I), TOTOUT GO TO 120 110 CONTINUE MINUSJ = -J WRITE (OUTUNT,99995) J, DATA(MINUSJ) 120 CONTINUE IF (OMEGA.NE.ZERO(1)) GO TO 130 WRITE (OUTUNT,99998) RETURN C 130 CONTINUE WRITE (OUTUNT,99997) THETA RETURN C 99999 FORMAT (30H-ERROR OCCURS IN CALL TO SQRSL/20H RECHECK DATA FOR ER, * 5HRORS.) 99998 FORMAT (18H-THETA IS INFINITE) 99997 FORMAT (9H-THETA = , 1PE13.5) 99996 FORMAT (12H-OUTPUT NODE, 10X, 5HVALUE, 12X, 18HCONDITION OF PROBL, * 2HEM, 12X, 22HCONDITION OF ALGORITHM, 11X, 18HRELATIVE ERROR BOU, * 2HND) 99995 FORMAT (1X, I6, 4X, 1PE20.7, 4X, 1PE20.7, 13X, 1PE20.7, 12X, * 1PE20.7) C END SUBROUTINE COMBIN(COL, LIMIT, K, L, INIT, LDIM) COM 10 C C THIS ROUTINE SUCCESSIVELY COMPUTES COL, A VECTOR OF COLUMN INDICE C REPRESENTING ONE OF THE COMINATIONS OF K ITEMS TAKEN L AT A TIME. C INTEGER COL(LDIM), LIMIT(LDIM) C IF (INIT.NE.1) GO TO 20 DO 10 I=1,L COL(I) = I LIMIT(I) = K - L + I 10 CONTINUE RETURN 20 DO 30 II=1,L I = L + 1 - II COL(I) = COL(I) + 1 IF (COL(I).LE.LIMIT(I)) GO TO 40 30 CONTINUE 40 IF (I.EQ.L) RETURN IP1 = I + 1 DO 50 JJ=IP1,L J = JJ - I COL(JJ) = COL(I) + J 50 CONTINUE RETURN C END SUBROUTINE CONDNS(A, MDIM, M, N, NCOL, ACTIVE, EPS) CON 10 C C ROUTINE TO CONDENSE THE M BY N MATRIX A TO AN M BY NCOL MATRIX. C THE MATRIX A HAS ROW DIMENSION MDIM IN THE CALLING PROGRAM. C THE VECTOR ACTIVE INDICATES WHICH COLUMNS ARE CANDIDATES FOR C CONDENSATION. 1 = YES, 0 = NO. THE CONDENSED MATRIX IS RETURNED C IN THE FIRST NCOL COLUMNS OF A. TWO COLUMNS, U AND V, ARE C CONDENSED IF INNERPRODUCT(U,V)/TWONORM(U)/TWONORM(V) EQUALS C ONE TO WITHIN MACHINE EPSILON, OR IF EITHER TWO-NORM IS MUCH C SMALLER THAN THE OTHER. C C CONDNS CALLS THE FOLLOWING BLAS--- C SAXPY C SCOPY C SDOT C SNRM2 C C DIMENSION A(MDIM,N) INTEGER ACTIVE(N) DATA ONE, ZERO /1.0E0,0.0E0/ DO 10 I=1,N ACTIVE(I) = 1 10 CONTINUE NCOL = 0 DO 60 ICOL=1,N IF (ACTIVE(ICOL).EQ.0) GO TO 60 NCOL = NCOL + 1 ICOLP1 = ICOL + 1 IF (ICOLP1.GT.N) GO TO 50 DO 40 IC=ICOLP1,N IF (ACTIVE(IC).EQ.0) GO TO 40 SUM = SDOT(M,A(1,ICOL),1,A(1,IC),1) C C MUST CHECK IF EITHER 2-NORM IS CLOSE TO ZERO, IE. MUCH SMALLER C THAN THE OTHER. IF ONE IS, THAT COLUMN IS CONDENSED C SNICOL = SNRM2(M,A(1,ICOL),1) SNIC = SNRM2(M,A(1,IC),1) IF (SNIC.LE.(EPS*SNICOL)) GO TO 20 IF (.NOT.(SNICOL.LE.(EPS*SNIC))) GO TO 30 CALL SCOPY(M, A(1,IC), 1, A(1,ICOL), 1) 20 CONTINUE ACTIVE(IC) = 0 GO TO 40 30 CONTINUE C COSINE = SUM/SNICOL/SNIC ERROR = ABS(COSINE) - ONE IF (ABS(ERROR).GT.EPS) GO TO 40 C C IF THE ANGLE BETWEEN COLUMNS A(1,IC) AND A(1,ICOL) IS ZERO, IE. C COSINE(ANGLE) IS CLOSE TO ONE, THEN CONDENSE THOSE TWO COLUMNS C ACTIVE(IC) = 0 SIGN = ONE IF (COSINE.LT.ZERO) SIGN = -ONE CALL SAXPY(M, SIGN, A(1,IC), 1, A(1,ICOL), 1) 40 CONTINUE 50 IF (ICOL.EQ.NCOL) GO TO 60 CALL SCOPY(M, A(1,ICOL), 1, A(1,NCOL), 1) 60 CONTINUE RETURN END REAL FUNCTION RCOMB(NTOP, NBOT) RCO 10 C C THIS FUNCTION COMPUTES THE NUMBER OF WAYS NTOP ITEMS CAN BE C TAKEN NBOT AT A TIME. C DATA ONE /1.0E0/ C IDIF = NTOP - NBOT COMB = ONE DO 10 I=1,IDIF COMB = COMB*FLOAT(NBOT+I)/FLOAT(I) 10 CONTINUE RCOMB = COMB RETURN END SUBROUTINE L1(M, N, MDIMA, N2, A, B, TOLER, X, E, S) L1 10 C C REFERENCES --- C C 1. BARRODALE, I. AND ROBERTS, F. D. R. AN IMPROVED ALGORITHM FOR C DISCRETE L1 LINEAR APPROXIMATION. SIAM J. NUMER. ANAL. C 10(1973), 839-848. C 2. BARRODALE, I. AND ROBERTS, F. D. R. ALGORITHM 478. SOLUTION C OF AN OVERDETERMINED SYSTEM OF EQUATIONS IN THE L1 NORM. C COMMUN. ACM 17(1974), 319-320. C C C THIS SUBROUTINE USES A MODIFICATION OF THE SIMPLEX METHOD C OF LINEAR PROGRAMMING TO CALCULATE AN L1 SOLUTION TO AN C OVER-DETERMINED SYSTEM OF LINEAR EQUATIONS. C DESCRIPTION OF PARAMETERS. C M NUMBER OF EQUATIONS. C N NUMBER OF UNKNOWNS (M.GE.N). C MDIMA ROW DIMENSION OF A IN CALLING PROGRAM C M2 SET EQUAL TO M+2 FOR ADJUSTABLE DIMENSIONS. C N2 SET EQUAL TO N+2 FOR ADJUSTABLE DIMENSIONS. C A TWO DIMENSIONAL REAL ARRAY OF SIZE (M2,N2). C ON ENTRY, THE COEFFICIENTS OF THE MATRIX MUST BE C STORED IN THE FIRST M ROWS AND N COLUMNS OF A. C THESE VALUES ARE DESTROYED BY THE SUBROUTINE. C B ONE DIMENSIONAL REAL ARRAY OF SIZE M. ON ENTRY, B C MUST CONTAIN THE RIGHT HAND SIDE OF THE EQUATIONS. C THESE VALUES ARE DESTROYED BY THE SUBROUTINE. C TOLER A SMALL POSITIVE TOLERANCE. EMPERICAL EVIDENCE C SUGGESTS TOLER=10**(-D*2/3) WHERE D REPRESENTS C THE NUMBER OF DECIMAL DIGITS OF ACCURACY AVAILABLE. C X ONE DIMENSIONAL REAL ARRAY OF SIZE N. ON EXIT, THIS C ARRAY CONTAINS A SOLUTION TO THE L1 PROBLEM. C E ONE DIMENSIONAL REAL ARRAY OF SIZE M. ON EXIT, THIS C ARRAY CONTAINS THE RESIDUALS IN THE EQUATIONS. C S INTEGER ARRAY OF SIZE M USED FOR WORKSPACE. C ON EXIT FROM THE SUBROUTINE, THE ARRAY A CONTAINS THE C FOLLOWING INFORMATION. C A(M+1,N+1) THE MINIMUM SUM OF THE ABSOLUTE VALUES OF C THE RESIDUALS. C A(M+1,N+2) THE RANK OF THE MATRIX OF COEFFICIENTS. C A(M+2,N+1) EXIT CODE WITH VALUES. C 0 - OPTIMAL SOLUTION WHICH IS PROBABLY NON-UNIQUE. C 1 - UNIQUE OPTIMAL SOLUTION. C 2 - CALCULATIONS TERMINATED PREMATURELY DUE TO C ROUNDING ERRORS. C A(M+2,N+2) NUMBER OF SIMPLEX ITERATIONS PERFORMED. C REAL MIN, MAX, A(MDIMA,N2), X(N), E(M), B(M) REAL ZERO(1), ONE(1) INTEGER OUT, S(M) LOGICAL STAGE, TEST C C BIG MUST BE SET EQUAL TO ANY VERY LARGE REAL CONSTANT. C ITS VALUE HERE IS APPROPRIATELY SMALL FOR MOST MACHINES C HAVING FLOATING POINT HARDWARE. C C L1 CALLS THE FOLLOWING BLAS--- C SAXPY C SCOPY C SDOT C SSCAL C SSWAP C C DATA BIG /1.E38/, ZERO(1) /0.0E0/, ONE(1) /1.0E0/, TWO /2.0E0/ C C INITIALIZATION. C M1 = M + 1 M2 = M + 2 N1 = N + 1 DO 10 J=1,N A(M2,J) = J 10 CONTINUE CALL SCOPY(N, ZERO(1), 0, X(1), 1) DO 20 I=1,M A(I,N2) = N + I A(I,N1) = B(I) IF (B(I).LT.ZERO(1)) CALL SSCAL(N2, -ONE(1), A(I,1), MDIMA) 20 CONTINUE CALL SCOPY(M, ZERO(1), 0, E(1), 1) C C COMPUTE THE MARGINAL COSTS. C DO 30 J=1,N1 A(M1,J) = SDOT(M,ONE(1),0,A(1,J),1) 30 CONTINUE C C STAGE I. C DETERMINE THE VECTOR TO ENTER THE BASIS. C STAGE = .TRUE. KOUNT = 0 KR = 1 KL = 1 40 MAX = -ONE(1) DO 50 J=KR,N IF (ABS(A(M2,J)).GT.FLOAT(N)) GO TO 50 D = ABS(A(M1,J)) IF (D.LE.MAX) GO TO 50 MAX = D IN = J 50 CONTINUE IF (A(M1,IN).GE.ZERO(1)) GO TO 60 CALL SSCAL(M2, -ONE(1), A(1,IN), 1) C C DETERMINE THE VECTOR TO LEAVE THE BASIS. C 60 K = 0 DO 70 I=KL,M D = A(I,IN) IF (D.LE.TOLER) GO TO 70 K = K + 1 B(K) = A(I,N1)/D S(K) = I TEST = .TRUE. 70 CONTINUE 80 IF (K.GT.0) GO TO 90 TEST = .FALSE. GO TO 110 90 MIN = BIG DO 100 I=1,K IF (B(I).GE.MIN) GO TO 100 J = I MIN = B(I) OUT = S(I) 100 CONTINUE B(J) = B(K) S(J) = S(K) K = K - 1 C C CHECK FOR LINEAR DEPENDENCE IN STAGE I. C 110 IF (TEST .OR. .NOT.STAGE) GO TO 120 CALL SSWAP(M2, A(1,KR), 1, A(1,IN), 1) KR = KR + 1 GO TO 160 120 IF (TEST) GO TO 130 A(M2,N1) = TWO GO TO 230 130 PIVOT = A(OUT,IN) IF (A(M1,IN)-PIVOT-PIVOT.LE.TOLER) GO TO 140 IF (M1.NE.OUT) CALL SAXPY(N1-KR+1, -TWO, A(OUT,KR), MDIMA, * A(M1,KR), MDIMA) CALL SSCAL(N1-KR+1, -ONE(1), A(OUT,KR), MDIMA) A(OUT,N2) = -A(OUT,N2) GO TO 80 C C PIVOT ON A(OUT,IN) C 140 ASAVE = A(OUT,IN) CALL SSCAL(N1-KR+1, ONE(1)/PIVOT, A(OUT,KR), MDIMA) A(OUT,IN) = ASAVE DO 150 J=KR,N1 IF (J.EQ.IN) GO TO 150 C C FORM A(*,J) = A(*,J) - A(OUT,J) * A(1,IN), EXCEPT FOR C FOR ELEMENT A(OUT,J) WHICH REMAINS UNCHANGED C TEMP = A(OUT,J) CALL SAXPY(M1, -A(OUT,J), A(1,IN), 1, A(1,J), 1) A(OUT,J) = TEMP 150 CONTINUE CALL SSCAL(M1, -ONE(1)/PIVOT, A(1,IN), 1) A(OUT,IN) = ONE(1)/PIVOT D = A(OUT,N2) A(OUT,N2) = A(M2,IN) A(M2,IN) = D KOUNT = KOUNT + 1 IF (.NOT.STAGE) GO TO 170 C C INTERCHANGE ROWS IN STAGE I. C KL = KL + 1 CALL SSWAP(N2-KR+1, A(OUT,KR), MDIMA, A(KOUNT,KR), MDIMA) 160 IF (KOUNT+KR.NE.N1) GO TO 40 C C STAGE II. C STAGE = .FALSE. C C DETERMINE THE VECTOR TO ENTER THE BASIS. C 170 MAX = -BIG DO 190 J=KR,N D = A(M1,J) IF (D.GE.ZERO(1)) GO TO 180 IF (D.GT.(-ONE(1))) GO TO 190 D = -D - TWO 180 IF (D.LE.MAX) GO TO 190 MAX = D IN = J 190 CONTINUE IF (MAX.LE.TOLER) GO TO 200 IF (A(M1,IN).GT.ZERO(1)) GO TO 60 CALL SSCAL(M2, -ONE(1), A(1,IN), 1) A(M1,IN) = A(M1,IN) - TWO GO TO 60 C C PREPARE OUTPUT. C 200 L = KL - 1 DO 210 I=1,L IF (A(I,N1).LT.ZERO(1)) CALL SSCAL(N2-KR+1, -ONE(1), A(I,KR), * MDIMA) 210 CONTINUE A(M2,N1) = ZERO(1) IF (KR.NE.1) GO TO 230 DO 220 J=1,N D = ABS(A(M1,J)) IF (D.LE.TOLER .OR. TWO-D.LE.TOLER) GO TO 230 220 CONTINUE A(M2,N1) = ONE(1) 230 DO 260 I=1,M K = A(I,N2) D = A(I,N1) IF (K.GT.0) GO TO 240 K = -K D = -D 240 IF (I.GE.KL) GO TO 250 X(K) = D GO TO 260 250 K = K - N E(K) = D 260 CONTINUE A(M2,N2) = KOUNT A(M1,N2) = N1 - KR A(M1,N1) = SDOT(M-KL+1,ONE(1),0,A(KL,N1),1) RETURN END SUBROUTINE SDIFST(A, B, C) SDI 10 C C C = A - B RETURN END SUBROUTINE SQRDC(X, LDX, N, P, QRAUX, JPVT, WORK, JOB) SQR 10 INTEGER LDX, N, P, JOB INTEGER JPVT(1) REAL X(LDX,1), QRAUX(1), WORK(1) C C SQRDC USES HOUSEHOLDER TRANSFORMATIONS TO COMPUTE THE QR C FACTORIZATION OF AN N BY P MATRIX X. COLUMN PIVOTING C BASED ON THE 2-NORMS OF THE REDUCED COLUMNS MAY BE C PERFORMED AT THE USERS OPTION. C C ON ENTRY C C X REAL(LDX,P), WHERE LDX .GE. N. C X CONTAINS THE MATRIX WHOSE DECOMPOSITION IS TO BE C COMPUTED. C C LDX INTEGER. C LDX IS THE LEADING DIMENSION OF THE ARRAY X. C C N INTEGER. C N IS THE NUMBER OF ROWS OF THE MATRIX X. C C P INTEGER. C P IS THE NUMBER OF COLUMNS OF THE MATRIX X. C C JPVT INTEGER(P). C JPVT CONTAINS INTEGERS THAT CONTROL THE SELECTION C OF THE PIVOT COLUMNS. THE K-TH COLUMN X(K) OF X C IS PLACED IN ONE OF THREE CLASSES ACCORDING TO THE C VALUE OF JPVT(K). C C IF JPVT(K) .GT. 0, THEN X(K) IS AN INITIAL C COLUMN. C C IF JPVT(K) .EQ. 0, THEN X(K) IS A FREE COLUMN. C C IF JPVT(K) .LT. 0, THEN X(K) IS A FINAL COLUMN. C C BEFORE THE DECOMPOSITION IS COMPUTED, INITIAL COLUMNS C ARE MOVED TO THE BEGINNING OF THE ARRAY X AND FINAL C COLUMNS TO THE END. BOTH INITIAL AND FINAL COLUMNS C ARE FROZEN IN PLACE DURING THE COMPUTATION AND ONLY C FREE COLUMNS ARE MOVED. AT THE K-TH STAGE OF THE C REDUCTION, IF X(K) IS OCCUPIED BY A FREE COLUMN C IT IS INTERCHANGED WITH THE FREE COLUMN OF LARGEST C REDUCED NORM. JPVT IS NOT REFERENCED IF C JOB .EQ. 0. C C WORK REAL(P). C WORK IS A WORK ARRAY. WORK IS NOT REFERENCED IF C JOB .EQ. 0. C C JOB INTEGER. C JOB IS AN INTEGER THAT INITIATES COLUMN PIVOTING. C IF JOB .EQ. 0, NO PIVOTING IS DONE. C IF JOB .NE. 0, PIVOTING IS DONE. C C ON RETURN C C X X CONTAINS IN ITS UPPER TRIANGLE THE UPPER C TRIANGULAR MATRIX R OF THE QR FACTORIZATION. C BELOW ITS DIAGONAL X CONTAINS INFORMATION FROM C WHICH THE ORTHOGONAL PART OF THE DECOMPOSITION C CAN BE RECOVERED. NOTE THAT IF PIVOTING HAS C BEEN REQUESTED, THE DECOMPOSITION IS NOT THAT C OF THE ORIGINAL MATRIX X BUT THAT OF X C WITH ITS COLUMNS PERMUTED AS DESCRIBED BY JPVT. C C QRAUX REAL(P). C QRAUX CONTAINS FURTHER INFORMATION REQUIRED TO RECOVER C THE ORTHOGONAL PART OF THE DECOMPOSITION. C C JPVT JPVT(K) CONTAINS THE INDEX OF THE COLUMN OF THE C ORIGINAL MATRIX THAT HAS BEEN INTERCHANGED INTO C THE K-TH COLUMN, IF PIVOTING WAS REQUESTED. C C LINPACK. THIS VERSION DATED 08/14/78 . C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C C SQRDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS. C C BLAS SAXPY,SDOT,SSCAL,SSWAP,SNRM2 C FORTRAN ABS,AMAX1,MIN0,SQRT C C INTERNAL VARIABLES C INTEGER J, JP, L, LP1, LUP, MAXJ, PL, PU REAL MAXNRM, SNRM2, TT REAL SDOT, NRMXL, T LOGICAL NEGJ, SWAPJ C C PL = 1 PU = 0 IF (JOB.EQ.0) GO TO 60 C C PIVOTING HAS BEEN REQUESTED. REARRANGE THE COLUMNS C ACCORDING TO JPVT. C DO 20 J=1,P SWAPJ = JPVT(J).GT.0 NEGJ = JPVT(J).LT.0 JPVT(J) = J IF (NEGJ) JPVT(J) = -J IF (.NOT.SWAPJ) GO TO 10 IF (J.NE.PL) CALL SSWAP(N, X(1,PL), 1, X(1,J), 1) JPVT(J) = JPVT(PL) JPVT(PL) = J PL = PL + 1 10 CONTINUE 20 CONTINUE PU = P DO 50 JJ=1,P J = P - JJ + 1 IF (JPVT(J).GE.0) GO TO 40 JPVT(J) = -JPVT(J) IF (J.EQ.PU) GO TO 30 CALL SSWAP(N, X(1,PU), 1, X(1,J), 1) JP = JPVT(PU) JPVT(PU) = JPVT(J) JPVT(J) = JP 30 CONTINUE PU = PU - 1 40 CONTINUE 50 CONTINUE 60 CONTINUE C C COMPUTE THE NORMS OF THE FREE COLUMNS. C IF (PU.LT.PL) GO TO 80 DO 70 J=PL,PU QRAUX(J) = SNRM2(N,X(1,J),1) WORK(J) = QRAUX(J) 70 CONTINUE 80 CONTINUE C C PERFORM THE HOUSEHOLDER REDUCTION OF X. C LUP = MIN0(N,P) DO 200 L=1,LUP IF (L.LT.PL .OR. L.GE.PU) GO TO 120 C C LOCATE THE COLUMN OF LARGEST NORM AND BRING IT C INTO THE PIVOT POSITION. C MAXNRM = 0.0E0 MAXJ = L DO 100 J=L,PU IF (QRAUX(J).LE.MAXNRM) GO TO 90 MAXNRM = QRAUX(J) MAXJ = J 90 CONTINUE 100 CONTINUE IF (MAXJ.EQ.L) GO TO 110 CALL SSWAP(N, X(1,L), 1, X(1,MAXJ), 1) QRAUX(MAXJ) = QRAUX(L) WORK(MAXJ) = WORK(L) JP = JPVT(MAXJ) JPVT(MAXJ) = JPVT(L) JPVT(L) = JP 110 CONTINUE 120 CONTINUE QRAUX(L) = 0.0E0 IF (L.EQ.N) GO TO 190 C C COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L. C NRMXL = SNRM2(N-L+1,X(L,L),1) IF (NRMXL.EQ.0.0E0) GO TO 180 IF (X(L,L).NE.0.0E0) NRMXL = SIGN(NRMXL,X(L,L)) CALL SSCAL(N-L+1, 1.0E0/NRMXL, X(L,L), 1) X(L,L) = 1.0E0 + X(L,L) C C APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS, C UPDATING THE NORMS. C LP1 = L + 1 IF (P.LT.LP1) GO TO 170 DO 160 J=LP1,P T = -SDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L) CALL SAXPY(N-L+1, T, X(L,L), 1, X(L,J), 1) IF (J.LT.PL .OR. J.GT.PU) GO TO 150 IF (QRAUX(J).EQ.0.0E0) GO TO 150 TT = 1.0E0 - (ABS(X(L,J))/QRAUX(J))**2 TT = AMAX1(TT,0.0E0) T = TT TT = 1.0E0 + 0.05E0*TT*(QRAUX(J)/WORK(J))**2 IF (TT.EQ.1.0E0) GO TO 130 QRAUX(J) = QRAUX(J)*SQRT(T) GO TO 140 130 CONTINUE QRAUX(J) = SNRM2(N-L,X(L+1,J),1) WORK(J) = QRAUX(J) 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE C C SAVE THE TRANSFORMATION. C QRAUX(L) = X(L,L) X(L,L) = -NRMXL 180 CONTINUE 190 CONTINUE 200 CONTINUE RETURN END SUBROUTINE SQRSL(X, LDX, N, K, QRAUX, Y, QY, QTY, B, RSD, XB, SQR 10 * JOB, INFO) INTEGER LDX, N, K, JOB, INFO REAL X(LDX,1), QRAUX(1), Y(1), QY(1), QTY(1), B(1), RSD(1), XB(1) C C SQRSL APPLIES THE OUTPUT OF SQRDC TO COMPUTE COORDINATE C TRANSFORMATIONS, PROJECTIONS, AND LEAST SQUARES SOLUTIONS. C FOR K .LE. MIN(N,P), LET XK BE THE MATRIX C C XK = (X(JPVT(1)),X(JPVT(2)), ... ,X(JPVT(K))) C C FORMED FROM COLUMNNS JPVT(1), ... ,JPVT(K) OF THE ORIGINAL C N X P MATRIX X THAT WAS INPUT TO SQRDC (IF NO PIVOTING WAS C DONE, XK CONSISTS OF THE FIRST K COLUMNS OF X IN THEIR C ORIGINAL ORDER). SQRDC PRODUCES A FACTORED ORTHOGONAL MATRIX Q C AND AN UPPER TRIANGULAR MATRIX R SUCH THAT C C XK = Q * (R) C (0) C C THIS INFORMATION IS CONTAINED IN CODED FORM IN THE ARRAYS C X AND QRAUX. C C ON ENTRY C C X REAL(LDX,P). C X CONTAINS THE OUTPUT OF SQRDC. C C LDX INTEGER. C LDX IS THE LEADING DIMENSION OF THE ARRAY X. C C N INTEGER. C N IS THE NUMBER OF ROWS OF THE MATRIX XK. IT MUST C HAVE THE SAME VALUE AS N IN SQRDC. C C K INTEGER. C K IS THE NUMBER OF COLUMNS OF THE MATRIX XK. K C MUST NNOT BE GREATER THAN MIN(N,P), WHERE P IS THE C SAME AS IN THE CALLING SEQUENCE TO SQRDC. C C QRAUX REAL(P). C QRAUX CONTAINS THE AUXILIARY OUTPUT FROM SQRDC. C C Y REAL(N) C Y CONTAINS AN N-VECTOR THAT IS TO BE MANIPULATED C BY SQRSL. C C JOB INTEGER. C JOB SPECIFIES WHAT IS TO BE COMPUTED. JOB HAS C THE DECIMAL EXPANSION ABCDE, WITH THE FOLLOWING C MEANING. C C IF A.NE.0, COMPUTE QY. C IF B,C,D, OR E .NE. 0, COMPUTE QTY. C IF C.NE.0, COMPUTE B. C IF D.NE.0, COMPUTE RSD. C IF E.NE.0, COMPUTE XB. C C NOTE THAT A REQUEST TO COMPUTE B, RSD, OR XB C AUTOMATICALLY TRIGGERS THE COMPUTATION OF QTY, FOR C WHICH AN ARRAY MUST BE PROVIDED IN THE CALLING C SEQUENCE. C C ON RETURN C C QY REAL(N). C QY CONNTAINS Q*Y, IF ITS COMPUTATION HAS BEEN C REQUESTED. C C QTY REAL(N). C QTY CONTAINS TRANS(Q)*Y, IF ITS COMPUTATION HAS C BEEN REQUESTED. HERE TRANS(Q) IS THE C TRANSPOSE OF THE MATRIX Q. C C B REAL(K) C B CONTAINS THE SOLUTION OF THE LEAST SQUARES PROBLEM C C MINIMIZE NORM2(Y - XK*B), C C IF ITS COMPUTATION HAS BEEN REQUESTED. (NOTE THAT C IF PIVOTING WAS REQUESTED IN SQRDC, THE J-TH C COMPONENT OF B WILL BE ASSOCIATED WITH COLUMN JPVT(J) C OF THE ORIGINAL MATRIX X THAT WAS INPUT INTO SQRDC.) C C RSD REAL(N). C RSD CONTAINS THE LEAST SQUARES RESIDUAL Y - XK*B, C IF ITS COMPUTATION HAS BEEN REQUESTED. RSD IS C ALSO THE ORTHOGONAL PROJECTION OF Y ONTO THE C ORTHOGONAL COMPLEMENT OF THE COLUMN SPACE OF XK. C C XB REAL(N). C XB CONTAINS THE LEAST SQUARES APPROXIMATION XK*B, C IF ITS COMPUTATION HAS BEEN REQUESTED. XB IS ALSO C THE ORTHOGONAL PROJECTION OF Y ONTO THE COLUMN SPACE C OF X. C C INFO INTEGER. C INFO IS ZERO UNLESS THE COMPUTATION OF B HAS C BEEN REQUESTED AND R IS EXACTLY SINGULAR. IN C THIS CASE, INFO IS THE INDEX OF THE FIRST ZERO C DIAGONAL ELEMENT OF R AND B IS LEFT UNALTERED. C C THE PARAMETERS QY, QTY, B, RSD, AND XB ARE NOT REFERENCED C IF THEIR COMPUTATION IS NOT REQUESTED AND IN THIS CASE C CAN BE REPLACED BY DUMMY VARIABLES IN THE CALLING PROGRAM. C TO SAVE STORAGE, THE USER MAY IN SOME CASES USE THE SAME C ARRAY FOR DIFFERENT PARAMETERS IN THE CALLING SEQUENCE. A C FREQUENTLY OCCURING EXAMPLE IS WHEN ONE WISHES TO COMPUTE C ANY OF B, RSD, OR XB AND DOES NOT NEED Y OR QTY. IN THIS C CASE ONE MAY IDENTIFY Y, QTY, AND ONE OF B, RSD, OR XB, WHILE C PROVIDING SEPARATE ARRAYS FOR ANYTHING ELSE THAT IS TO BE C COMPUTED. THUS THE CALLING SEQUENCE C C CALL SQRSL(X,LDX,N,K,QRAUX,Y,DUM,Y,B,Y,DUM,110,INFO) C C WILL RESULT IN THE COMPUTATION OF B AND RSD, WITH RSD C OVERWRITING Y. MORE GENERALLY, EACH ITEM IN THE FOLLOWING C LIST CONTAINS GROUPS OF PERMISSIBLE IDENTIFICATIONS FOR C A SINGLE CALLINNG SEQUENCE. C C 1. (Y,QTY,B) (RSD) (XB) (QY) C C 2. (Y,QTY,RSD) (B) (XB) (QY) C C 3. (Y,QTY,XB) (B) (RSD) (QY) C C 4. (Y,QY) (QTY,B) (RSD) (XB) C C 5. (Y,QY) (QTY,RSD) (B) (XB) C C 6. (Y,QY) (QTY,XB) (B) (RSD) C C IN ANY GROUP THE VALUE RETURNED IN THE ARRAY ALLOCATED TO C THE GROUP CORRESPONDS TO THE LAST MEMBER OF THE GROUP. C C LINPACK. THIS VERSION DATED 08/14/78 . C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C C SQRSL USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS. C C BLAS SAXPY,SCOPY,SDOT C FORTRAN ABS,MIN0,MOD C C INTERNAL VARIABLES C INTEGER I, J, JJ, JU, KP1 REAL SDOT, T, TEMP LOGICAL CB, CQY, CQTY, CR, CXB C C C SET INFO FLAG. C INFO = 0 C C DETERMINE WHAT IS TO BE COMPUTED. C CQY = (JOB/10000).NE.0 CQTY = (MOD(JOB,10000)).NE.0 CB = (MOD(JOB,1000)/100).NE.0 CR = (MOD(JOB,100)/10).NE.0 CXB = (MOD(JOB,10)).NE.0 JU = MIN0(K,N-1) C C SPECIAL ACTION WHEN N=1. C IF (JU.NE.0) GO TO 40 IF (CQY) QY(1) = Y(1) IF (CQTY) QTY(1) = Y(1) IF (CXB) XB(1) = Y(1) IF (.NOT.CB) GO TO 30 IF (X(1,1).NE.0.0E0) GO TO 10 INFO = 1 GO TO 20 10 CONTINUE B(1) = Y(1)/X(1,1) 20 CONTINUE 30 CONTINUE IF (CR) RSD(1) = 0.0E0 GO TO 250 40 CONTINUE C C SET UP TO COMPUTE QY OR QTY. C IF (CQY) CALL SCOPY(N, Y, 1, QY, 1) IF (CQTY) CALL SCOPY(N, Y, 1, QTY, 1) IF (.NOT.CQY) GO TO 70 C C COMPUTE QY. C DO 60 JJ=1,JU J = JU - JJ + 1 IF (QRAUX(J).EQ.0.0E0) GO TO 50 TEMP = X(J,J) X(J,J) = QRAUX(J) T = -SDOT(N-J+1,X(J,J),1,QY(J),1)/X(J,J) CALL SAXPY(N-J+1, T, X(J,J), 1, QY(J), 1) X(J,J) = TEMP 50 CONTINUE 60 CONTINUE 70 CONTINUE IF (.NOT.CQTY) GO TO 100 C C COMPUTE TRANS(Q)*Y. C DO 90 J=1,JU IF (QRAUX(J).EQ.0.0E0) GO TO 80 TEMP = X(J,J) X(J,J) = QRAUX(J) T = -SDOT(N-J+1,X(J,J),1,QTY(J),1)/X(J,J) CALL SAXPY(N-J+1, T, X(J,J), 1, QTY(J), 1) X(J,J) = TEMP 80 CONTINUE 90 CONTINUE 100 CONTINUE C C SET UP TO COMPUTE B, RSD, OR XB. C IF (CB) CALL SCOPY(K, QTY, 1, B, 1) KP1 = K + 1 IF (CXB) CALL SCOPY(K, QTY, 1, XB, 1) IF (CR .AND. K.LT.N) CALL SCOPY(N-K, QTY(KP1), 1, RSD(KP1), 1) IF (.NOT.CXB .OR. KP1.GT.N) GO TO 120 DO 110 I=KP1,N XB(I) = 0.0E0 110 CONTINUE 120 CONTINUE IF (.NOT.CR) GO TO 140 DO 130 I=1,K RSD(I) = 0.0E0 130 CONTINUE 140 CONTINUE IF (.NOT.CB) GO TO 190 C C COMPUTE B. C DO 170 JJ=1,K J = K - JJ + 1 IF (X(J,J).NE.0.0E0) GO TO 150 INFO = J C ......EXIT GO TO 180 150 CONTINUE B(J) = B(J)/X(J,J) IF (J.EQ.1) GO TO 160 T = -B(J) CALL SAXPY(J-1, T, X(1,J), 1, B, 1) 160 CONTINUE 170 CONTINUE 180 CONTINUE 190 CONTINUE IF (.NOT.CR .AND. .NOT.CXB) GO TO 240 C C COMPUTE RSD OR XB AS REQUIRED. C DO 230 JJ=1,JU J = JU - JJ + 1 IF (QRAUX(J).EQ.0.0E0) GO TO 220 TEMP = X(J,J) X(J,J) = QRAUX(J) IF (.NOT.CR) GO TO 200 T = -SDOT(N-J+1,X(J,J),1,RSD(J),1)/X(J,J) CALL SAXPY(N-J+1, T, X(J,J), 1, RSD(J), 1) 200 CONTINUE IF (.NOT.CXB) GO TO 210 T = -SDOT(N-J+1,X(J,J),1,XB(J),1)/X(J,J) CALL SAXPY(N-J+1, T, X(J,J), 1, XB(J), 1) 210 CONTINUE X(J,J) = TEMP 220 CONTINUE 230 CONTINUE 240 CONTINUE 250 CONTINUE RETURN END SUBROUTINE SCOPY(N, SX, INCX, SY, INCY) SCO 10 C C COPY SINGLE PRECISION SX TO SINGLE PRECISION SY. C REAL SX(1), SY(1) C IF (N.LE.0) RETURN IF (INCX.EQ.INCY) IF (INCX-1) 10, 30, 70 10 CONTINUE C C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. C IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO 20 I=1,N SY(IY) = SX(IX) IX = IX + INCX IY = IY + INCY 20 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 7. C 30 M = N - (N/7)*7 IF (M.EQ.0) GO TO 50 DO 40 I=1,M SY(I) = SX(I) 40 CONTINUE IF (N.LT.7) RETURN 50 MP1 = M + 1 DO 60 I=MP1,N,7 SY(I) = SX(I) SY(I+1) = SX(I+1) SY(I+2) = SX(I+2) SY(I+3) = SX(I+3) SY(I+4) = SX(I+4) SY(I+5) = SX(I+5) SY(I+6) = SX(I+6) 60 CONTINUE RETURN C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C 70 CONTINUE NS = N*INCX DO 80 I=1,NS,INCX SY(I) = SX(I) 80 CONTINUE RETURN END SUBROUTINE SSWAP(N, SX, INCX, SY, INCY) SSW 10 C C INTERCHANGE SINGLE PRECISION SX AND SINGLE PRECISION SY. C REAL SX(1), SY(1), STEMP1, STEMP2, STEMP3 C IF (N.LE.0) RETURN IF (INCX.EQ.INCY) IF (INCX-1) 10, 30, 70 10 CONTINUE C C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. C IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO 20 I=1,N STEMP1 = SX(IX) SX(IX) = SY(IY) SY(IY) = STEMP1 IX = IX + INCX IY = IY + INCY 20 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 3. C 30 M = N - (N/3)*3 IF (M.EQ.0) GO TO 50 DO 40 I=1,M STEMP1 = SX(I) SX(I) = SY(I) SY(I) = STEMP1 40 CONTINUE IF (N.LT.3) RETURN 50 MP1 = M + 1 DO 60 I=MP1,N,3 STEMP1 = SX(I) STEMP2 = SX(I+1) STEMP3 = SX(I+2) SX(I) = SY(I) SX(I+1) = SY(I+1) SX(I+2) = SY(I+2) SY(I) = STEMP1 SY(I+1) = STEMP2 SY(I+2) = STEMP3 60 CONTINUE RETURN 70 CONTINUE C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C NS = N*INCX DO 80 I=1,NS,INCX STEMP1 = SX(I) SX(I) = SY(I) SY(I) = STEMP1 80 CONTINUE RETURN END REAL FUNCTION SNRM2(N, SX, INCX) SNR 10 INTEGER NEXT REAL SX(1), CUTLO, CUTHI, HITEST, SUM, XMAX, ZERO, ONE DATA ZERO, ONE /0.0E0,1.0E0/ C C EUCLIDEAN NORM OF THE N-VECTOR STORED IN SX() WITH STORAGE C INCREMENT INCX . C IF N .LE. 0 RETURN WITH RESULT = 0. C IF N .GE. 1 THEN INCX MUST BE .GE. 1 C C C.L.LAWSON, 1978 JAN 08 C C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE C HOPEFULLY APPLICABLE TO ALL MACHINES. C CUTLO = MAXIMUM OF SQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF SQRT(V) OVER ALL KNOWN MACHINES. C WHERE C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT) C V = LARGEST NO. (OVERFLOW LIMIT) C C BRIEF OUTLINE OF ALGORITHM.. C C PHASE 1 SCANS ZERO COMPONENTS. C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. C C VALUES FOR CUTLO AND CUTHI.. C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE C UNIVAC AND DEC AT 2**(-103) C THUS CUTLO = 2**(-51) = 4.44089E-16 C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. C THUS CUTHI = 2**(63.5) = 1.30438E19 C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. C THUS CUTLO = 2**(-33.5) = 8.23181D-11 C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19 C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / DATA CUTLO, CUTHI /4.441E-16,1.304E19/ C IF (N.GT.0) GO TO 10 SNRM2 = ZERO GO TO 140 C 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N*INCX C BEGIN MAIN LOOP I = 1 20 GO TO NEXT, (30, 40, 70, 80) 30 IF (ABS(SX(I)).GT.CUTLO) GO TO 110 ASSIGN 40 TO NEXT XMAX = ZERO C C PHASE 1. SUM IS ZERO C 40 IF (SX(I).EQ.ZERO) GO TO 130 IF (ABS(SX(I)).GT.CUTLO) GO TO 110 C C PREPARE FOR PHASE 2. ASSIGN 70 TO NEXT GO TO 60 C C PREPARE FOR PHASE 4. C 50 I = J ASSIGN 80 TO NEXT SUM = (SUM/SX(I))/SX(I) 60 XMAX = ABS(SX(I)) GO TO 90 C C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 IF (ABS(SX(I)).GT.CUTLO) GO TO 100 C C COMMON CODE FOR PHASES 2 AND 4. C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. C 80 IF (ABS(SX(I)).LE.XMAX) GO TO 90 SUM = ONE + SUM*(XMAX/SX(I))**2 XMAX = ABS(SX(I)) GO TO 130 C 90 SUM = SUM + (SX(I)/XMAX)**2 GO TO 130 C C C PREPARE FOR PHASE 3. C 100 SUM = (SUM*XMAX)*XMAX C C C FOR REAL OR D.P. SET HITEST = CUTHI/N C FOR COMPLEX SET HITEST = CUTHI/(2*N) C 110 HITEST = CUTHI/FLOAT(N) C C PHASE 3. SUM IS MID-RANGE. NO SCALING. C DO 120 J=I,NN,INCX IF (ABS(SX(J)).GE.HITEST) GO TO 50 SUM = SUM + SX(J)**2 120 CONTINUE SNRM2 = SQRT(SUM) GO TO 140 C 130 CONTINUE I = I + INCX IF (I.LE.NN) GO TO 20 C C END OF MAIN LOOP. C C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. C SNRM2 = XMAX*SQRT(SUM) 140 CONTINUE RETURN END REAL FUNCTION SASUM(N, SX, INCX) SAS 10 C C RETURNS SUM OF MAGNITUDES OF SINGLE PRECISION SX. C REAL SX(1) C SASUM = 0.0E0 IF (N.LE.0) RETURN IF (INCX.EQ.1) GO TO 20 C C CODE FOR INCREMENTS NOT EQUAL TO 1. C NS = N*INCX DO 10 I=1,NS,INCX SASUM = SASUM + ABS(SX(I)) 10 CONTINUE RETURN C C CODE FOR INCREMENTS EQUAL TO 1. C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 6. C 20 M = N - (N/6)*6 IF (M.EQ.0) GO TO 40 DO 30 I=1,M SASUM = SASUM + ABS(SX(I)) 30 CONTINUE IF (N.LT.6) RETURN 40 MP1 = M + 1 DO 50 I=MP1,N,6 SASUM = SASUM + ABS(SX(I)) + ABS(SX(I+1)) + ABS(SX(I+2)) + * ABS(SX(I+3)) + ABS(SX(I+4)) + ABS(SX(I+5)) 50 CONTINUE RETURN END SUBROUTINE SSCAL(N, SA, SX, INCX) SSC 10 C C REPLACE SINGLE PRECISION SX BY SINGLE PRECISION SA*SX. C REAL SA, SX(1) C IF (N.LE.0) RETURN IF (INCX.EQ.1) GO TO 20 C C CODE FOR INCREMENTS NOT EQUAL TO 1. C NS = N*INCX DO 10 I=1,NS,INCX SX(I) = SA*SX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENTS EQUAL TO 1. C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5. C 20 M = N - (N/5)*5 IF (M.EQ.0) GO TO 40 DO 30 I=1,M SX(I) = SA*SX(I) 30 CONTINUE IF (N.LT.5) RETURN 40 MP1 = M + 1 DO 50 I=MP1,N,5 SX(I) = SA*SX(I) SX(I+1) = SA*SX(I+1) SX(I+2) = SA*SX(I+2) SX(I+3) = SA*SX(I+3) SX(I+4) = SA*SX(I+4) 50 CONTINUE RETURN END INTEGER FUNCTION ISAMAX(N, SX, INCX) ISA 10 C C FIND SMALLEST INDEX OF MAXIMUM MAGNITUDE OF SINGLE PRECISION SX. C REAL SX(1), SMAX, XMAG C ISAMAX = 0 IF (N.LE.0) RETURN ISAMAX = 1 IF (N.LE.1) RETURN IF (INCX.EQ.1) GO TO 30 C C CODE FOR INCREMENTS NOT EQUAL TO 1. C SMAX = ABS(SX(1)) NS = N*INCX II = 1 DO 20 I=1,NS,INCX XMAG = ABS(SX(I)) IF (XMAG.LE.SMAX) GO TO 10 ISAMAX = II SMAX = XMAG 10 II = II + 1 20 CONTINUE RETURN C C CODE FOR INCREMENTS EQUAL TO 1. C 30 SMAX = ABS(SX(1)) DO 40 I=2,N XMAG = ABS(SX(I)) IF (XMAG.LE.SMAX) GO TO 40 ISAMAX = I SMAX = XMAG 40 CONTINUE RETURN END REAL FUNCTION SDOT(N, SX, INCX, SY, INCY) SDO 10 C C RETURNS THE DOT PRODUCT OF SINGLE PRECISION SX AND SY. C REAL SX(1), SY(1) SDOT = 0.0E0 IF (N.LE.0) RETURN IF (INCX.EQ.INCY) IF (INCX-1) 10, 30, 70 10 CONTINUE C C CODE FOR UNEQUAL INCREMENTS OR NONPOSITIVE INCREMENTS. C IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO 20 I=1,N SDOT = SDOT + SX(IX)*SY(IY) IX = IX + INCX IY = IY + INCY 20 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5. C 30 M = N - (N/5)*5 IF (M.EQ.0) GO TO 50 DO 40 I=1,M SDOT = SDOT + SX(I)*SY(I) 40 CONTINUE IF (N.LT.5) RETURN 50 MP1 = M + 1 DO 60 I=MP1,N,5 SDOT = SDOT + SX(I)*SY(I) + SX(I+1)*SY(I+1) + SX(I+2)*SY(I+2) + * SX(I+3)*SY(I+3) + SX(I+4)*SY(I+4) 60 CONTINUE RETURN C C CODE FOR POSITIVE EQUAL INCREMENTS .NE.1. C 70 CONTINUE NS = N*INCX DO 80 I=1,NS,INCX SDOT = SDOT + SX(I)*SY(I) 80 CONTINUE RETURN END SUBROUTINE SAXPY(N, SA, SX, INCX, SY, INCY) SAX 10 C C OVERWRITE SINGLE PRECISION SY WITH SINGLE PRECISION SA*SX +SY. C REAL SX(1), SY(1), SA IF (N.LE.0 .OR. SA.EQ.0.E0) RETURN IF (INCX.EQ.INCY) IF (INCX-1) 10, 30, 70 10 CONTINUE C C CODE FOR NONEQUAL OR NONPOSITIVE INCREMENTS. C IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO 20 I=1,N SY(IY) = SY(IY) + SA*SX(IX) IX = IX + INCX IY = IY + INCY 20 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 4. C 30 M = N - (N/4)*4 IF (M.EQ.0) GO TO 50 DO 40 I=1,M SY(I) = SY(I) + SA*SX(I) 40 CONTINUE IF (N.LT.4) RETURN 50 MP1 = M + 1 DO 60 I=MP1,N,4 SY(I) = SY(I) + SA*SX(I) SY(I+1) = SY(I+1) + SA*SX(I+1) SY(I+2) = SY(I+2) + SA*SX(I+2) SY(I+3) = SY(I+3) + SA*SX(I+3) 60 CONTINUE RETURN C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C 70 CONTINUE NS = N*INCX DO 80 I=1,NS,INCX SY(I) = SA*SX(I) + SY(I) 80 CONTINUE RETURN END SUBROUTINE QZHES(NM, N, A, B, MATZ, Z) QZH 10 C INTEGER I, J, K, L, N, LB, L1, NM, NK1, NM1, NM2 REAL A(NM,N), B(NM,N), Z(NM,N) REAL R, S, T, U1, U2, V1, V2, RHO C REAL SQRT,ABS,SIGN LOGICAL MATZ C C THIS SUBROUTINE IS THE FIRST STEP OF THE QZ ALGORITHM C FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, C SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. C C THIS SUBROUTINE ACCEPTS A PAIR OF REAL GENERAL MATRICES AND C REDUCES ONE OF THEM TO UPPER HESSENBERG FORM AND THE OTHER C TO UPPER TRIANGULAR FORM USING ORTHOGONAL TRANSFORMATIONS. C IT IS USUALLY FOLLOWED BY QZIT, QZVAL AND, POSSIBLY, QZVEC. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRICES, C C A CONTAINS A REAL GENERAL MATRIX, C C B CONTAINS A REAL GENERAL MATRIX, C C MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS C ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING C EIGENVECTORS, AND TO .FALSE. OTHERWISE. C C ON OUTPUT- C C A HAS BEEN REDUCED TO UPPER HESSENBERG FORM. THE ELEMENTS C BELOW THE FIRST SUBDIAGONAL HAVE BEEN SET TO ZERO, C C B HAS BEEN REDUCED TO UPPER TRIANGULAR FORM. THE ELEMENTS C BELOW THE MAIN DIAGONAL HAVE BEEN SET TO ZERO, C C Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS IF C MATZ HAS BEEN SET TO .TRUE. OTHERWISE, Z IS NOT REFERENCED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C C ********** INITIALIZE Z ********** IF (.NOT.MATZ) GO TO 30 C DO 20 I=1,N C DO 10 J=1,N Z(I,J) = 0.0 10 CONTINUE C Z(I,I) = 1.0 20 CONTINUE C ********** REDUCE B TO UPPER TRIANGULAR FORM ********** 30 IF (N.LE.1) GO TO 210 NM1 = N - 1 C DO 130 L=1,NM1 L1 = L + 1 S = 0.0 C DO 40 I=L1,N S = S + ABS(B(I,L)) 40 CONTINUE C IF (S.EQ.0.0) GO TO 130 S = S + ABS(B(L,L)) R = 0.0 C DO 50 I=L,N B(I,L) = B(I,L)/S R = R + B(I,L)**2 50 CONTINUE C R = SIGN(SQRT(R),B(L,L)) B(L,L) = B(L,L) + R RHO = R*B(L,L) C DO 80 J=L1,N T = 0.0 C DO 60 I=L,N T = T + B(I,L)*B(I,J) 60 CONTINUE C T = -T/RHO C DO 70 I=L,N B(I,J) = B(I,J) + T*B(I,L) 70 CONTINUE C 80 CONTINUE C DO 110 J=1,N T = 0.0 C DO 90 I=L,N T = T + B(I,L)*A(I,J) 90 CONTINUE C T = -T/RHO C DO 100 I=L,N A(I,J) = A(I,J) + T*B(I,L) 100 CONTINUE C 110 CONTINUE C B(L,L) = -S*R C DO 120 I=L1,N B(I,L) = 0.0 120 CONTINUE C 130 CONTINUE C ********** REDUCE A TO UPPER HESSENBERG FORM, WHILE C KEEPING B TRIANGULAR ********** IF (N.EQ.2) GO TO 210 NM2 = N - 2 C DO 200 K=1,NM2 NK1 = NM1 - K C ********** FOR L=N-1 STEP -1 UNTIL K+1 DO -- ********** DO 190 LB=1,NK1 L = N - LB L1 = L + 1 C ********** ZERO A(L+1,K) ********** S = ABS(A(L,K)) + ABS(A(L1,K)) IF (S.EQ.0.0) GO TO 190 U1 = A(L,K)/S U2 = A(L1,K)/S R = SIGN(SQRT(U1*U1+U2*U2),U1) V1 = -(U1+R)/R V2 = -U2/R U2 = V2/V1 C DO 140 J=K,N T = A(L,J) + U2*A(L1,J) A(L,J) = A(L,J) + T*V1 A(L1,J) = A(L1,J) + T*V2 140 CONTINUE C A(L1,K) = 0.0 C DO 150 J=L,N T = B(L,J) + U2*B(L1,J) B(L,J) = B(L,J) + T*V1 B(L1,J) = B(L1,J) + T*V2 150 CONTINUE C ********** ZERO B(L+1,L) ********** S = ABS(B(L1,L1)) + ABS(B(L1,L)) IF (S.EQ.0.0) GO TO 190 U1 = B(L1,L1)/S U2 = B(L1,L)/S R = SIGN(SQRT(U1*U1+U2*U2),U1) V1 = -(U1+R)/R V2 = -U2/R U2 = V2/V1 C DO 160 I=1,L1 T = B(I,L1) + U2*B(I,L) B(I,L1) = B(I,L1) + T*V1 B(I,L) = B(I,L) + T*V2 160 CONTINUE C B(L1,L) = 0.0 C DO 170 I=1,N T = A(I,L1) + U2*A(I,L) A(I,L1) = A(I,L1) + T*V1 A(I,L) = A(I,L) + T*V2 170 CONTINUE C IF (.NOT.MATZ) GO TO 190 C DO 180 I=1,N T = Z(I,L1) + U2*Z(I,L) Z(I,L1) = Z(I,L1) + T*V1 Z(I,L) = Z(I,L) + T*V2 180 CONTINUE C 190 CONTINUE C 200 CONTINUE C 210 RETURN C ********** LAST CARD OF QZHES ********** END SUBROUTINE QZIT(NM, N, A, B, EPS1, MATZ, Z, IERR) QZI 10 C INTEGER I, J, K, L, N, EN, K1, K2, LD, LL, L1, NA, NM, ISH, ITS, * KM1, LM1, ENM2, IERR, LOR1, ENORN REAL A(NM,N), B(NM,N), Z(NM,N) REAL R, S, T, A1, A2, A3, EP, SH, U1, U2, U3, V1, V2, V3, ANI, * A11, A12, A21, A22, A33, A34, A43, A44, BNI, B11, B12, B22, B33, * B34, B44, EPSA, EPSB, EPS1, ANORM, BNORM C REAL SQRT,ABS,SIGN C INTEGER MAX0,MIN0 LOGICAL MATZ, NOTLAS C C THIS SUBROUTINE IS THE SECOND STEP OF THE QZ ALGORITHM C FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, C SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART, C AS MODIFIED IN TECHNICAL NOTE NASA TN E-7305(1973) BY WARD. C C THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM C IN UPPER HESSENBERG FORM AND THE OTHER IN UPPER TRIANGULAR FORM. C IT REDUCES THE HESSENBERG MATRIX TO QUASI-TRIANGULAR FORM USING C ORTHOGONAL TRANSFORMATIONS WHILE MAINTAINING THE TRIANGULAR FORM C OF THE OTHER MATRIX. IT IS USUALLY PRECEDED BY QZHES AND C FOLLOWED BY QZVAL AND, POSSIBLY, QZVEC. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRICES, C C A CONTAINS A REAL UPPER HESSENBERG MATRIX, C C B CONTAINS A REAL UPPER TRIANGULAR MATRIX, C C EPS1 IS A TOLERANCE USED TO DETERMINE NEGLIGIBLE ELEMENTS. C EPS1 = 0.0 (OR NEGATIVE) MAY BE INPUT, IN WHICH CASE AN C ELEMENT WILL BE NEGLECTED ONLY IF IT IS LESS THAN ROUNDOFF C ERROR TIMES THE NORM OF ITS MATRIX. IF THE INPUT EPS1 IS C POSITIVE, THEN AN ELEMENT WILL BE CONSIDERED NEGLIGIBLE C IF IT IS LESS THAN EPS1 TIMES THE NORM OF ITS MATRIX. A C POSITIVE VALUE OF EPS1 MAY RESULT IN FASTER EXECUTION, C BUT LESS ACCURATE RESULTS, C C MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS C ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING C EIGENVECTORS, AND TO .FALSE. OTHERWISE, C C Z CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE C TRANSFORMATION MATRIX PRODUCED IN THE REDUCTION C BY QZHES, IF PERFORMED, OR ELSE THE IDENTITY MATRIX. C IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED. C C ON OUTPUT- C C A HAS BEEN REDUCED TO QUASI-TRIANGULAR FORM. THE ELEMENTS C BELOW THE FIRST SUBDIAGONAL ARE STILL ZERO AND NO TWO C CONSECUTIVE SUBDIAGONAL ELEMENTS ARE NONZERO, C C B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS C HAVE BEEN ALTERED. THE LOCATION B(N,1) IS USED TO STORE C EPS1 TIMES THE NORM OF B FOR LATER USE BY QZVAL AND QZVEC, C C Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS C (FOR BOTH STEPS) IF MATZ HAS BEEN SET TO .TRUE., C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF NEITHER A(J,J-1) NOR A(J-1,J-2) HAS BECOME C ZERO AFTER 50 ITERATIONS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C IERR = 0 C ********** COMPUTE EPSA,EPSB ********** ANORM = 0.0 BNORM = 0.0 C DO 20 I=1,N ANI = 0.0 IF (I.NE.1) ANI = ABS(A(I,I-1)) BNI = 0.0 C DO 10 J=I,N ANI = ANI + ABS(A(I,J)) BNI = BNI + ABS(B(I,J)) 10 CONTINUE C IF (ANI.GT.ANORM) ANORM = ANI IF (BNI.GT.BNORM) BNORM = BNI 20 CONTINUE C IF (ANORM.EQ.0.0) ANORM = 1.0 IF (BNORM.EQ.0.0) BNORM = 1.0 EP = EPS1 IF (EP.GT.0.0) GO TO 40 C ********** COMPUTE ROUNDOFF LEVEL IF EPS1 IS ZERO ********** EP = 1.0 30 EP = EP/2.0 IF (1.0+EP.GT.1.0) GO TO 30 40 EPSA = EP*ANORM EPSB = EP*BNORM C ********** REDUCE A TO QUASI-TRIANGULAR FORM, WHILE C KEEPING B TRIANGULAR ********** LOR1 = 1 ENORN = N EN = N C ********** BEGIN QZ STEP ********** 50 IF (EN.LE.2) GO TO 310 IF (.NOT.MATZ) ENORN = EN ITS = 0 NA = EN - 1 ENM2 = NA - 1 60 ISH = 2 C ********** CHECK FOR CONVERGENCE OR REDUCIBILITY. C FOR L=EN STEP -1 UNTIL 1 DO -- ********** DO 70 LL=1,EN LM1 = EN - LL L = LM1 + 1 IF (L.EQ.1) GO TO 90 IF (ABS(A(L,LM1)).LE.EPSA) GO TO 80 70 CONTINUE C 80 A(L,LM1) = 0.0 IF (L.LT.NA) GO TO 90 C ********** 1-BY-1 OR 2-BY-2 BLOCK ISOLATED ********** EN = LM1 GO TO 50 C ********** CHECK FOR SMALL TOP OF B ********** 90 LD = L 100 L1 = L + 1 B11 = B(L,L) IF (ABS(B11).GT.EPSB) GO TO 120 B(L,L) = 0.0 S = ABS(A(L,L)) + ABS(A(L1,L)) U1 = A(L,L)/S U2 = A(L1,L)/S R = SIGN(SQRT(U1*U1+U2*U2),U1) V1 = -(U1+R)/R V2 = -U2/R U2 = V2/V1 C DO 110 J=L,ENORN T = A(L,J) + U2*A(L1,J) A(L,J) = A(L,J) + T*V1 A(L1,J) = A(L1,J) + T*V2 T = B(L,J) + U2*B(L1,J) B(L,J) = B(L,J) + T*V1 B(L1,J) = B(L1,J) + T*V2 110 CONTINUE C IF (L.NE.1) A(L,LM1) = -A(L,LM1) LM1 = L L = L1 GO TO 80 120 A11 = A(L,L)/B11 A21 = A(L1,L)/B11 IF (ISH.EQ.1) GO TO 140 C ********** ITERATION STRATEGY ********** IF (ITS.EQ.50) GO TO 300 IF (ITS.EQ.10) GO TO 160 C ********** DETERMINE TYPE OF SHIFT ********** B22 = B(L1,L1) IF (ABS(B22).LT.EPSB) B22 = EPSB B33 = B(NA,NA) IF (ABS(B33).LT.EPSB) B33 = EPSB B44 = B(EN,EN) IF (ABS(B44).LT.EPSB) B44 = EPSB A33 = A(NA,NA)/B33 A34 = A(NA,EN)/B44 A43 = A(EN,NA)/B33 A44 = A(EN,EN)/B44 B34 = B(NA,EN)/B44 T = 0.5*(A43*B34-A33-A44) R = T*T + A34*A43 - A33*A44 IF (R.LT.0.0) GO TO 150 C ********** DETERMINE SINGLE SHIFT ZEROTH COLUMN OF A ********** ISH = 1 R = SQRT(R) SH = -T + R S = -T - R IF (ABS(S-A44).LT.ABS(SH-A44)) SH = S C ********** LOOK FOR TWO CONSECUTIVE SMALL C SUB-DIAGONAL ELEMENTS OF A. C FOR L=EN-2 STEP -1 UNTIL LD DO -- ********** DO 130 LL=LD,ENM2 L = ENM2 + LD - LL IF (L.EQ.LD) GO TO 140 LM1 = L - 1 L1 = L + 1 T = A(L,L) IF (ABS(B(L,L)).GT.EPSB) T = T - SH*B(L,L) IF (ABS(A(L,LM1)).LE.ABS(T/A(L1,L))*EPSA) GO TO 100 130 CONTINUE C 140 A1 = A11 - SH A2 = A21 IF (L.NE.LD) A(L,LM1) = -A(L,LM1) GO TO 170 C ********** DETERMINE DOUBLE SHIFT ZEROTH COLUMN OF A ********** 150 A12 = A(L,L1)/B22 A22 = A(L1,L1)/B22 B12 = B(L,L1)/B22 A1 = ((A33-A11)*(A44-A11)-A34*A43+A43*B34*A11)/A21 + A12 - A11*B12 A2 = (A22-A11) - A21*B12 - (A33-A11) - (A44-A11) + A43*B34 A3 = A(L1+1,L1)/B22 GO TO 170 C ********** AD HOC SHIFT ********** 160 A1 = 0.0 A2 = 1.0 A3 = 1.1605 170 ITS = ITS + 1 IF (.NOT.MATZ) LOR1 = LD C ********** MAIN LOOP ********** DO 290 K=L,NA NOTLAS = K.NE.NA .AND. ISH.EQ.2 K1 = K + 1 K2 = K + 2 KM1 = MAX0(K-1,L) LL = MIN0(EN,K1+ISH) IF (NOTLAS) GO TO 200 C ********** ZERO A(K+1,K-1) ********** IF (K.EQ.L) GO TO 180 A1 = A(K,KM1) A2 = A(K1,KM1) 180 S = ABS(A1) + ABS(A2) IF (S.EQ.0.0) GO TO 60 U1 = A1/S U2 = A2/S R = SIGN(SQRT(U1*U1+U2*U2),U1) V1 = -(U1+R)/R V2 = -U2/R U2 = V2/V1 C DO 190 J=KM1,ENORN T = A(K,J) + U2*A(K1,J) A(K,J) = A(K,J) + T*V1 A(K1,J) = A(K1,J) + T*V2 T = B(K,J) + U2*B(K1,J) B(K,J) = B(K,J) + T*V1 B(K1,J) = B(K1,J) + T*V2 190 CONTINUE C IF (K.NE.L) A(K1,KM1) = 0.0 GO TO 260 C ********** ZERO A(K+1,K-1) AND A(K+2,K-1) ********** 200 IF (K.EQ.L) GO TO 210 A1 = A(K,KM1) A2 = A(K1,KM1) A3 = A(K2,KM1) 210 S = ABS(A1) + ABS(A2) + ABS(A3) IF (S.EQ.0.0) GO TO 290 U1 = A1/S U2 = A2/S U3 = A3/S R = SIGN(SQRT(U1*U1+U2*U2+U3*U3),U1) V1 = -(U1+R)/R V2 = -U2/R V3 = -U3/R U2 = V2/V1 U3 = V3/V1 C DO 220 J=KM1,ENORN T = A(K,J) + U2*A(K1,J) + U3*A(K2,J) A(K,J) = A(K,J) + T*V1 A(K1,J) = A(K1,J) + T*V2 A(K2,J) = A(K2,J) + T*V3 T = B(K,J) + U2*B(K1,J) + U3*B(K2,J) B(K,J) = B(K,J) + T*V1 B(K1,J) = B(K1,J) + T*V2 B(K2,J) = B(K2,J) + T*V3 220 CONTINUE C IF (K.EQ.L) GO TO 230 A(K1,KM1) = 0.0 A(K2,KM1) = 0.0 C ********** ZERO B(K+2,K+1) AND B(K+2,K) ********** 230 S = ABS(B(K2,K2)) + ABS(B(K2,K1)) + ABS(B(K2,K)) IF (S.EQ.0.0) GO TO 260 U1 = B(K2,K2)/S U2 = B(K2,K1)/S U3 = B(K2,K)/S R = SIGN(SQRT(U1*U1+U2*U2+U3*U3),U1) V1 = -(U1+R)/R V2 = -U2/R V3 = -U3/R U2 = V2/V1 U3 = V3/V1 C DO 240 I=LOR1,LL T = A(I,K2) + U2*A(I,K1) + U3*A(I,K) A(I,K2) = A(I,K2) + T*V1 A(I,K1) = A(I,K1) + T*V2 A(I,K) = A(I,K) + T*V3 T = B(I,K2) + U2*B(I,K1) + U3*B(I,K) B(I,K2) = B(I,K2) + T*V1 B(I,K1) = B(I,K1) + T*V2 B(I,K) = B(I,K) + T*V3 240 CONTINUE C B(K2,K) = 0.0 B(K2,K1) = 0.0 IF (.NOT.MATZ) GO TO 260 C DO 250 I=1,N T = Z(I,K2) + U2*Z(I,K1) + U3*Z(I,K) Z(I,K2) = Z(I,K2) + T*V1 Z(I,K1) = Z(I,K1) + T*V2 Z(I,K) = Z(I,K) + T*V3 250 CONTINUE C ********** ZERO B(K+1,K) ********** 260 S = ABS(B(K1,K1)) + ABS(B(K1,K)) IF (S.EQ.0.0) GO TO 290 U1 = B(K1,K1)/S U2 = B(K1,K)/S R = SIGN(SQRT(U1*U1+U2*U2),U1) V1 = -(U1+R)/R V2 = -U2/R U2 = V2/V1 C DO 270 I=LOR1,LL T = A(I,K1) + U2*A(I,K) A(I,K1) = A(I,K1) + T*V1 A(I,K) = A(I,K) + T*V2 T = B(I,K1) + U2*B(I,K) B(I,K1) = B(I,K1) + T*V1 B(I,K) = B(I,K) + T*V2 270 CONTINUE C B(K1,K) = 0.0 IF (.NOT.MATZ) GO TO 290 C DO 280 I=1,N T = Z(I,K1) + U2*Z(I,K) Z(I,K1) = Z(I,K1) + T*V1 Z(I,K) = Z(I,K) + T*V2 280 CONTINUE C 290 CONTINUE C ********** END QZ STEP ********** GO TO 60 C ********** SET ERROR -- NEITHER BOTTOM SUBDIAGONAL ELEMENT C HAS BECOME NEGLIGIBLE AFTER 50 ITERATIONS ********** 300 IERR = EN C ********** SAVE EPSB FOR USE BY QZVAL AND QZVEC ********** 310 IF (N.GT.1) B(N,1) = EPSB RETURN C ********** LAST CARD OF QZIT ********** END SUBROUTINE QZVAL(NM, N, A, B, ALFR, ALFI, BETA, MATZ, Z) QZV 10 C INTEGER I, J, N, EN, NA, NM, NN, ISW REAL A(NM,N), B(NM,N), ALFR(N), ALFI(N), BETA(N), Z(NM,N) REAL C, D, E, R, S, T, AN, A1, A2, BN, CQ, CZ, DI, DR, EI, TI, * TR, U1, U2, V1, V2, A1I, A11, A12, A2I, A21, A22, B11, B12, B22, * SQI, SQR, SSI, SSR, SZI, SZR, A11I, A11R, A12I, A12R, A22I, * A22R, EPSB C REAL SQRT,ABS,SIGN LOGICAL MATZ C C THIS SUBROUTINE IS THE THIRD STEP OF THE QZ ALGORITHM C FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, C SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. C C THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM C IN QUASI-TRIANGULAR FORM AND THE OTHER IN UPPER TRIANGULAR FORM. C IT REDUCES THE QUASI-TRIANGULAR MATRIX FURTHER, SO THAT ANY C REMAINING 2-BY-2 BLOCKS CORRESPOND TO PAIRS OF COMPLEX C EIGENVALUES, AND RETURNS QUANTITIES WHOSE RATIOS GIVE THE C GENERALIZED EIGENVALUES. IT IS USUALLY PRECEDED BY QZHES C AND QZIT AND MAY BE FOLLOWED BY QZVEC. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRICES, C C A CONTAINS A REAL UPPER QUASI-TRIANGULAR MATRIX, C C B CONTAINS A REAL UPPER TRIANGULAR MATRIX. IN ADDITION, C LOCATION B(N,1) CONTAINS THE TOLERANCE QUANTITY (EPSB) C COMPUTED AND SAVED IN QZIT, C C MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS C ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING C EIGENVECTORS, AND TO .FALSE. OTHERWISE, C C Z CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE C TRANSFORMATION MATRIX PRODUCED IN THE REDUCTIONS BY QZHES C AND QZIT, IF PERFORMED, OR ELSE THE IDENTITY MATRIX. C IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED. C C ON OUTPUT- C C A HAS BEEN REDUCED FURTHER TO A QUASI-TRIANGULAR MATRIX C IN WHICH ALL NONZERO SUBDIAGONAL ELEMENTS CORRESPOND TO C PAIRS OF COMPLEX EIGENVALUES, C C B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS C HAVE BEEN ALTERED. B(N,1) IS UNALTERED, C C ALFR AND ALFI CONTAIN THE REAL AND IMAGINARY PARTS OF THE C DIAGONAL ELEMENTS OF THE TRIANGULAR MATRIX THAT WOULD BE C OBTAINED IF A WERE REDUCED COMPLETELY TO TRIANGULAR FORM C BY UNITARY TRANSFORMATIONS. NON-ZERO VALUES OF ALFI OCCUR C IN PAIRS, THE FIRST MEMBER POSITIVE AND THE SECOND NEGATIVE, C C BETA CONTAINS THE DIAGONAL ELEMENTS OF THE CORRESPONDING B, C NORMALIZED TO BE REAL AND NON-NEGATIVE. THE GENERALIZED C EIGENVALUES ARE THEN THE RATIOS ((ALFR+I*ALFI)/BETA), C C Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS C (FOR ALL THREE STEPS) IF MATZ HAS BEEN SET TO .TRUE. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C EPSB = B(N,1) ISW = 1 C ********** FIND EIGENVALUES OF QUASI-TRIANGULAR MATRICES. C FOR EN=N STEP -1 UNTIL 1 DO -- ********** DO 260 NN=1,N EN = N + 1 - NN NA = EN - 1 IF (ISW.EQ.2) GO TO 250 IF (EN.EQ.1) GO TO 10 IF (A(EN,NA).NE.0.0) GO TO 20 C ********** 1-BY-1 BLOCK, ONE REAL ROOT ********** 10 ALFR(EN) = A(EN,EN) IF (B(EN,EN).LT.0.0) ALFR(EN) = -ALFR(EN) BETA(EN) = ABS(B(EN,EN)) ALFI(EN) = 0.0 GO TO 260 C ********** 2-BY-2 BLOCK ********** 20 IF (ABS(B(NA,NA)).LE.EPSB) GO TO 100 IF (ABS(B(EN,EN)).GT.EPSB) GO TO 30 A1 = A(EN,EN) A2 = A(EN,NA) BN = 0.0 GO TO 60 30 AN = ABS(A(NA,NA)) + ABS(A(NA,EN)) + ABS(A(EN,NA)) + * ABS(A(EN,EN)) BN = ABS(B(NA,NA)) + ABS(B(NA,EN)) + ABS(B(EN,EN)) A11 = A(NA,NA)/AN A12 = A(NA,EN)/AN A21 = A(EN,NA)/AN A22 = A(EN,EN)/AN B11 = B(NA,NA)/BN B12 = B(NA,EN)/BN B22 = B(EN,EN)/BN E = A11/B11 EI = A22/B22 S = A21/(B11*B22) T = (A22-E*B22)/B22 IF (ABS(E).LE.ABS(EI)) GO TO 40 E = EI T = (A11-E*B11)/B11 40 C = 0.5*(T-S*B12) D = C*C + S*(A12-E*B12) IF (D.LT.0.0) GO TO 140 C ********** TWO REAL ROOTS. C ZERO BOTH A(EN,NA) AND B(EN,NA) ********** E = E + (C+SIGN(SQRT(D),C)) A11 = A11 - E*B11 A12 = A12 - E*B12 A22 = A22 - E*B22 IF (ABS(A11)+ABS(A12).LT.ABS(A21)+ABS(A22)) GO TO 50 A1 = A12 A2 = A11 GO TO 60 50 A1 = A22 A2 = A21 C ********** CHOOSE AND APPLY REAL Z ********** 60 S = ABS(A1) + ABS(A2) U1 = A1/S U2 = A2/S R = SIGN(SQRT(U1*U1+U2*U2),U1) V1 = -(U1+R)/R V2 = -U2/R U2 = V2/V1 C DO 70 I=1,EN T = A(I,EN) + U2*A(I,NA) A(I,EN) = A(I,EN) + T*V1 A(I,NA) = A(I,NA) + T*V2 T = B(I,EN) + U2*B(I,NA) B(I,EN) = B(I,EN) + T*V1 B(I,NA) = B(I,NA) + T*V2 70 CONTINUE C IF (.NOT.MATZ) GO TO 90 C DO 80 I=1,N T = Z(I,EN) + U2*Z(I,NA) Z(I,EN) = Z(I,EN) + T*V1 Z(I,NA) = Z(I,NA) + T*V2 80 CONTINUE C 90 IF (BN.EQ.0.0) GO TO 130 IF (AN.LT.ABS(E)*BN) GO TO 100 A1 = B(NA,NA) A2 = B(EN,NA) GO TO 110 100 A1 = A(NA,NA) A2 = A(EN,NA) C ********** CHOOSE AND APPLY REAL Q ********** 110 S = ABS(A1) + ABS(A2) IF (S.EQ.0.0) GO TO 130 U1 = A1/S U2 = A2/S R = SIGN(SQRT(U1*U1+U2*U2),U1) V1 = -(U1+R)/R V2 = -U2/R U2 = V2/V1 C DO 120 J=NA,N T = A(NA,J) + U2*A(EN,J) A(NA,J) = A(NA,J) + T*V1 A(EN,J) = A(EN,J) + T*V2 T = B(NA,J) + U2*B(EN,J) B(NA,J) = B(NA,J) + T*V1 B(EN,J) = B(EN,J) + T*V2 120 CONTINUE C 130 A(EN,NA) = 0.0 B(EN,NA) = 0.0 ALFR(NA) = A(NA,NA) ALFR(EN) = A(EN,EN) IF (B(NA,NA).LT.0.0) ALFR(NA) = -ALFR(NA) IF (B(EN,EN).LT.0.0) ALFR(EN) = -ALFR(EN) BETA(NA) = ABS(B(NA,NA)) BETA(EN) = ABS(B(EN,EN)) ALFI(EN) = 0.0 ALFI(NA) = 0.0 GO TO 250 C ********** TWO COMPLEX ROOTS ********** 140 E = E + C EI = SQRT(-D) A11R = A11 - E*B11 A11I = EI*B11 A12R = A12 - E*B12 A12I = EI*B12 A22R = A22 - E*B22 A22I = EI*B22 IF (ABS(A11R)+ABS(A11I)+ABS(A12R)+ABS(A12I).LT.ABS(A21) * +ABS(A22R)+ABS(A22I)) GO TO 150 A1 = A12R A1I = A12I A2 = -A11R A2I = -A11I GO TO 160 150 A1 = A22R A1I = A22I A2 = -A21 A2I = 0.0 C ********** CHOOSE COMPLEX Z ********** 160 CZ = SQRT(A1*A1+A1I*A1I) IF (CZ.EQ.0.0) GO TO 170 SZR = (A1*A2+A1I*A2I)/CZ SZI = (A1*A2I-A1I*A2)/CZ R = SQRT(CZ*CZ+SZR*SZR+SZI*SZI) CZ = CZ/R SZR = SZR/R SZI = SZI/R GO TO 180 170 SZR = 1.0 SZI = 0.0 180 IF (AN.LT.(ABS(E)+EI)*BN) GO TO 190 A1 = CZ*B11 + SZR*B12 A1I = SZI*B12 A2 = SZR*B22 A2I = SZI*B22 GO TO 200 190 A1 = CZ*A11 + SZR*A12 A1I = SZI*A12 A2 = CZ*A21 + SZR*A22 A2I = SZI*A22 C ********** CHOOSE COMPLEX Q ********** 200 CQ = SQRT(A1*A1+A1I*A1I) IF (CQ.EQ.0.0) GO TO 210 SQR = (A1*A2+A1I*A2I)/CQ SQI = (A1*A2I-A1I*A2)/CQ R = SQRT(CQ*CQ+SQR*SQR+SQI*SQI) CQ = CQ/R SQR = SQR/R SQI = SQI/R GO TO 220 210 SQR = 1.0 SQI = 0.0 C ********** COMPUTE DIAGONAL ELEMENTS THAT WOULD RESULT C IF TRANSFORMATIONS WERE APPLIED ********** 220 SSR = SQR*SZR + SQI*SZI SSI = SQR*SZI - SQI*SZR I = 1 TR = CQ*CZ*A11 + CQ*SZR*A12 + SQR*CZ*A21 + SSR*A22 TI = CQ*SZI*A12 - SQI*CZ*A21 + SSI*A22 DR = CQ*CZ*B11 + CQ*SZR*B12 + SSR*B22 DI = CQ*SZI*B12 + SSI*B22 GO TO 240 230 I = 2 TR = SSR*A11 - SQR*CZ*A12 - CQ*SZR*A21 + CQ*CZ*A22 TI = -SSI*A11 - SQI*CZ*A12 + CQ*SZI*A21 DR = SSR*B11 - SQR*CZ*B12 + CQ*CZ*B22 DI = -SSI*B11 - SQI*CZ*B12 240 T = TI*DR - TR*DI J = NA IF (T.LT.0.0) J = EN R = SQRT(DR*DR+DI*DI) BETA(J) = BN*R ALFR(J) = AN*(TR*DR+TI*DI)/R ALFI(J) = AN*T/R IF (I.EQ.1) GO TO 230 250 ISW = 3 - ISW 260 CONTINUE C RETURN C ********** LAST CARD OF QZVAL ********** END SUBROUTINE QZVEC(NM, N, A, B, ALFR, ALFI, BETA, Z) QZV 10 C INTEGER I, J, K, M, N, EN, II, JJ, NA, NM, NN, ISW, ENM2 REAL A(NM,N), B(NM,N), ALFR(N), ALFI(N), BETA(N), Z(NM,N) REAL D, Q, R, S, T, W, X, Y, DI, DR, RA, RR, SA, TI, TR, T1, T2, * W1, X1, ZZ, Z1, ALFM, ALMI, ALMR, BETM, EPSB C REAL SQRT,ABS C C THIS SUBROUTINE IS THE OPTIONAL FOURTH STEP OF THE QZ ALGORITHM C FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, C SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. C C THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM IN C QUASI-TRIANGULAR FORM (IN WHICH EACH 2-BY-2 BLOCK CORRESPONDS TO C A PAIR OF COMPLEX EIGENVALUES) AND THE OTHER IN UPPER TRIANGULAR C FORM. IT COMPUTES THE EIGENVECTORS OF THE TRIANGULAR PROBLEM AND C TRANSFORMS THE RESULTS BACK TO THE ORIGINAL COORDINATE SYSTEM. C IT IS USUALLY PRECEDED BY QZHES, QZIT, AND QZVAL. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRICES, C C A CONTAINS A REAL UPPER QUASI-TRIANGULAR MATRIX, C C B CONTAINS A REAL UPPER TRIANGULAR MATRIX. IN ADDITION, C LOCATION B(N,1) CONTAINS THE TOLERANCE QUANTITY (EPSB) C COMPUTED AND SAVED IN QZIT, C C ALFR, ALFI, AND BETA ARE VECTORS WITH COMPONENTS WHOSE C RATIOS ((ALFR+I*ALFI)/BETA) ARE THE GENERALIZED C EIGENVALUES. THEY ARE USUALLY OBTAINED FROM QZVAL, C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTIONS BY QZHES, QZIT, AND QZVAL, IF PERFORMED. C IF THE EIGENVECTORS OF THE TRIANGULAR PROBLEM ARE C DESIRED, Z MUST CONTAIN THE IDENTITY MATRIX. C C ON OUTPUT- C C A IS UNALTERED. ITS SUBDIAGONAL ELEMENTS PROVIDE INFORMATION C ABOUT THE STORAGE OF THE COMPLEX EIGENVECTORS, C C B HAS BEEN DESTROYED, C C ALFR, ALFI, AND BETA ARE UNALTERED, C C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. C IF ALFI(I) .EQ. 0.0, THE I-TH EIGENVALUE IS REAL AND C THE I-TH COLUMN OF Z CONTAINS ITS EIGENVECTOR. C IF ALFI(I) .NE. 0.0, THE I-TH EIGENVALUE IS COMPLEX. C IF ALFI(I) .GT. 0.0, THE EIGENVALUE IS THE FIRST OF C A COMPLEX PAIR AND THE I-TH AND (I+1)-TH COLUMNS C OF Z CONTAIN ITS EIGENVECTOR. C IF ALFI(I) .LT. 0.0, THE EIGENVALUE IS THE SECOND OF C A COMPLEX PAIR AND THE (I-1)-TH AND I-TH COLUMNS C OF Z CONTAIN THE CONJUGATE OF ITS EIGENVECTOR. C EACH EIGENVECTOR IS NORMALIZED SO THAT THE MODULUS C OF ITS LARGEST COMPONENT IS 1.0 . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C EPSB = B(N,1) ISW = 1 C ********** FOR EN=N STEP -1 UNTIL 1 DO -- ********** DO 190 NN=1,N EN = N + 1 - NN NA = EN - 1 IF (ISW.EQ.2) GO TO 180 IF (ALFI(EN).NE.0.0) GO TO 70 C ********** REAL VECTOR ********** M = EN B(EN,EN) = 1.0 IF (NA.EQ.0) GO TO 190 ALFM = ALFR(M) BETM = BETA(M) C ********** FOR I=EN-1 STEP -1 UNTIL 1 DO -- ********** DO 60 II=1,NA I = EN - II W = BETM*A(I,I) - ALFM*B(I,I) R = 0.0 C DO 10 J=M,EN R = R + (BETM*A(I,J)-ALFM*B(I,J))*B(J,EN) 10 CONTINUE C IF (I.EQ.1 .OR. ISW.EQ.2) GO TO 20 IF (BETM*A(I,I-1).EQ.0.0) GO TO 20 ZZ = W S = R GO TO 50 20 M = I IF (ISW.EQ.2) GO TO 30 C ********** REAL 1-BY-1 BLOCK ********** T = W IF (W.EQ.0.0) T = EPSB B(I,EN) = -R/T GO TO 60 C ********** REAL 2-BY-2 BLOCK ********** 30 X = BETM*A(I,I+1) - ALFM*B(I,I+1) Y = BETM*A(I+1,I) Q = W*ZZ - X*Y T = (X*S-ZZ*R)/Q B(I,EN) = T IF (ABS(X).LE.ABS(ZZ)) GO TO 40 B(I+1,EN) = (-R-W*T)/X GO TO 50 40 B(I+1,EN) = (-S-Y*T)/ZZ 50 ISW = 3 - ISW 60 CONTINUE C ********** END REAL VECTOR ********** GO TO 190 C ********** COMPLEX VECTOR ********** 70 M = NA ALMR = ALFR(M) ALMI = ALFI(M) BETM = BETA(M) C ********** LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT C EIGENVECTOR MATRIX IS TRIANGULAR ********** Y = BETM*A(EN,NA) B(NA,NA) = -ALMI*B(EN,EN)/Y B(NA,EN) = (ALMR*B(EN,EN)-BETM*A(EN,EN))/Y B(EN,NA) = 0.0 B(EN,EN) = 1.0 ENM2 = NA - 1 IF (ENM2.EQ.0) GO TO 180 C ********** FOR I=EN-2 STEP -1 UNTIL 1 DO -- ********** DO 170 II=1,ENM2 I = NA - II W = BETM*A(I,I) - ALMR*B(I,I) W1 = -ALMI*B(I,I) RA = 0.0 SA = 0.0 C DO 80 J=M,EN X = BETM*A(I,J) - ALMR*B(I,J) X1 = -ALMI*B(I,J) RA = RA + X*B(J,NA) - X1*B(J,EN) SA = SA + X*B(J,EN) + X1*B(J,NA) 80 CONTINUE C IF (I.EQ.1 .OR. ISW.EQ.2) GO TO 90 IF (BETM*A(I,I-1).EQ.0.0) GO TO 90 ZZ = W Z1 = W1 R = RA S = SA ISW = 2 GO TO 170 90 M = I IF (ISW.EQ.2) GO TO 130 C ********** COMPLEX 1-BY-1 BLOCK ********** TR = -RA TI = -SA 100 DR = W DI = W1 C ********** COMPLEX DIVIDE (T1,T2) = (TR,TI) / (DR,DI) ********** 110 IF (ABS(DI).GT.ABS(DR)) GO TO 120 RR = DI/DR D = DR + DI*RR T1 = (TR+TI*RR)/D T2 = (TI-TR*RR)/D GO TO (160, 140), ISW 120 RR = DR/DI D = DR*RR + DI T1 = (TR*RR+TI)/D T2 = (TI*RR-TR)/D GO TO (160, 140), ISW C ********** COMPLEX 2-BY-2 BLOCK ********** 130 X = BETM*A(I,I+1) - ALMR*B(I,I+1) X1 = -ALMI*B(I,I+1) Y = BETM*A(I+1,I) TR = Y*RA - W*R + W1*S TI = Y*SA - W*S - W1*R DR = W*ZZ - W1*Z1 - X*Y DI = W*Z1 + W1*ZZ - X1*Y IF (DR.EQ.0.0 .AND. DI.EQ.0.0) DR = EPSB GO TO 110 140 B(I+1,NA) = T1 B(I+1,EN) = T2 ISW = 1 IF (ABS(Y).GT.ABS(W)+ABS(W1)) GO TO 150 TR = -RA - X*B(I+1,NA) + X1*B(I+1,EN) TI = -SA - X*B(I+1,EN) - X1*B(I+1,NA) GO TO 100 150 T1 = (-R-ZZ*B(I+1,NA)+Z1*B(I+1,EN))/Y T2 = (-S-ZZ*B(I+1,EN)-Z1*B(I+1,NA))/Y 160 B(I,NA) = T1 B(I,EN) = T2 170 CONTINUE C ********** END COMPLEX VECTOR ********** 180 ISW = 3 - ISW 190 CONTINUE C ********** END BACK SUBSTITUTION. C TRANSFORM TO ORIGINAL COORDINATE SYSTEM. C FOR J=N STEP -1 UNTIL 1 DO -- ********** DO 220 JJ=1,N J = N + 1 - JJ C DO 210 I=1,N ZZ = 0.0 C DO 200 K=1,J ZZ = ZZ + Z(I,K)*B(K,J) 200 CONTINUE C Z(I,J) = ZZ 210 CONTINUE 220 CONTINUE C ********** NORMALIZE SO THAT MODULUS OF LARGEST C COMPONENT OF EACH VECTOR IS 1. C (ISW IS 1 INITIALLY FROM BEFORE) ********** DO 290 J=1,N D = 0.0 IF (ISW.EQ.2) GO TO 250 IF (ALFI(J).NE.0.0) GO TO 280 C DO 230 I=1,N IF (ABS(Z(I,J)).GT.D) D = ABS(Z(I,J)) 230 CONTINUE C DO 240 I=1,N Z(I,J) = Z(I,J)/D 240 CONTINUE C GO TO 290 C 250 DO 260 I=1,N R = ABS(Z(I,J-1)) + ABS(Z(I,J)) IF (R.NE.0.0) R = R*SQRT((Z(I,J-1)/R)**2+(Z(I,J)/R)**2) IF (R.GT.D) D = R 260 CONTINUE C DO 270 I=1,N Z(I,J-1) = Z(I,J-1)/D Z(I,J) = Z(I,J)/D 270 CONTINUE C 280 ISW = 3 - ISW 290 CONTINUE C RETURN C ********** LAST CARD OF QZVEC ********** END 62 1 GAUSSIAN ELIMINATION---DATA 1 2 4 1 101 3 5 6 2 102 101 3 9 10 2 104 101 3 13 14 2 106 101 3 17 18 2 108 3 4 1 110 3 5 7 2 111 110 3 9 11 2 113 110 3 13 15 2 115 110 3 17 19 2 117 4 4 1 119 3 5 8 2 120 119 3 9 12 2 122 119 3 13 16 2 124 119 3 17 20 2 126 112 4 103 128 3 105 114 2 129 128 3 107 116 2 131 128 3 109 118 2 133 121 4 103 135 3 105 123 2 136 135 3 107 125 2 138 135 3 109 127 2 140 137 4 130 142 3 132 139 2 143 142 3 134 141 2 145 146 4 144 132 3 147 134 2 148 149 4 130 105 3 150 107 3 147 151 1 152 109 2 153 154 4 103 5 3 155 9 3 150 156 1 157 13 3 147 158 1 159 17 2 160 161 4 1 4 62 55 50 47 1 0.0 20 DATA 1 2.908 1.212 4.552 -5.809 -2.253 1.995 5.681 -5.030 6.775 2.266 8.850 0.099 3.970 8.008 1.302 7.832 6.291 7.219 5.730 9.574 20 DATA 2 1.0 0.5 0.25 0.125 1.0 0.25 0.0625 0.015625 1.0 0.125 0.015625 0.0019531 1.0 0.015625 0.00390625 0.00024414 10.0 1.625 0.4375 0.16308594 20 DATA 3 6.0 4.0 4.0 1.0 4.0 6.0 1.0 4.0 4.0 1.0 6.0 4.0 1.0 4.0 4.0 6.0 30.0 35.0 40.0 45.0 20 DATA 4 4.0 2.0 1.0 3.0 1.0 -5.0 2.0 1.0 1.0 1.0 6.0 2.0 1.0 1.0 1.0 -7.0 13.0 -1.0 27.0 -17.0 20 DATA 5 10.0 3.0 2.0 1.0 3.0 11.0 3.0 2.0 2.0 3.0 12.0 3.0 1.0 2.0 3.0 13.0 26.0 42.0 56.0 66.0 -1 76 1 GAUSS-JORDAN---DATA 1 2 4 1 101 3 5 6 2 102 101 3 9 10 2 104 101 3 13 14 2 106 101 3 17 18 2 108 3 4 1 110 3 5 7 2 111 110 3 9 11 2 113 110 3 13 15 2 115 110 3 17 19 2 117 4 4 1 119 3 5 8 2 120 119 3 9 12 2 122 119 3 13 16 2 124 119 3 17 20 2 126 5 4 103 128 3 105 9 2 129 128 3 107 13 2 131 128 3 109 17 2 133 112 4 103 135 3 105 114 2 136 135 3 107 116 2 138 135 3 109 118 2 140 121 4 103 142 3 105 123 2 143 142 3 107 125 2 145 142 3 109 127 2 147 130 4 137 149 3 139 132 2 150 149 3 141 134 2 152 105 4 137 154 3 139 107 2 155 154 3 141 109 2 157 144 4 137 159 3 139 146 2 160 159 3 141 148 2 162 151 4 161 164 3 163 153 2 165 156 4 161 167 3 163 158 2 168 139 4 161 170 3 163 141 2 171 166 4 1 169 4 103 172 4 137 163 4 161 4 73 74 75 76 0 20 DATA 1 2.908 1.212 4.552 -5.809 -2.253 1.995 5.681 -5.030 6.775 2.266 8.850 0.099 3.970 8.008 1.302 7.832 6.291 7.219 5.730 9.574 20 DATA 2 1.0 0.5 0.25 0.125 1.0 0.25 0.0625 0.015625 1.0 0.125 0.015625 0.0019531 1.0 0.015625 0.00390625 0.00024414 10.0 1.625 0.4375 0.16308594 20 DATA 3 6.0 4.0 4.0 1.0 4.0 6.0 1.0 4.0 4.0 1.0 6.0 4.0 1.0 4.0 4.0 6.0 30.0 35.0 40.0 45.0 20 DATA 4 4.0 2.0 1.0 3.0 1.0 -5.0 2.0 1.0 1.0 1.0 6.0 2.0 1.0 1.0 1.0 -7.0 13.0 -1.0 27.0 -17.0 20 DATA 5 10.0 3.0 2.0 1.0 3.0 11.0 3.0 2.0 2.0 3.0 12.0 3.0 1.0 2.0 3.0 13.0 26.0 42.0 56.0 66.0 -1 168 1 HOUSEHOLDER'S METHOD---DATA 1 1 3 1 2 3 2 101 1 102 3 3 3 103 1 104 4 3 4 105 1 106 107 5 0 108 6 0 1 1 108 108 3 110 110 3 5 2 3 6 112 1 113 3 3 7 114 1 115 4 3 8 116 1 117 118 4 111 119 3 110 5 2 120 119 3 2 6 2 122 119 3 3 7 2 124 119 3 4 8 2 126 110 3 9 2 3 10 128 1 129 3 3 11 130 1 131 4 3 12 132 1 133 134 4 111 135 3 110 9 2 136 135 3 2 10 2 138 135 3 3 11 2 140 135 3 4 12 2 142 110 3 13 2 3 14 144 1 145 3 3 15 146 1 147 4 3 16 148 1 149 150 4 111 151 3 110 13 2 152 151 3 2 14 2 154 151 3 3 15 2 156 151 3 4 16 2 158 110 3 17 2 3 18 160 1 161 3 3 19 162 1 163 4 3 20 164 1 165 166 4 111 167 3 110 17 2 168 167 3 2 18 2 170 167 3 3 19 2 172 167 3 4 20 2 174 123 3 123 125 3 125 176 1 177 127 3 127 178 1 179 180 5 0 181 6 0 123 1 181 181 3 183 183 3 139 125 3 141 185 1 186 127 3 143 187 1 188 189 4 184 190 3 183 139 2 191 190 3 125 141 2 193 190 3 127 143 2 195 183 3 155 125 3 157 197 1 198 127 3 159 199 1 200 201 4 184 202 3 183 155 2 203 202 3 125 157 2 205 202 3 127 159 2 207 183 3 171 125 3 173 209 1 210 127 3 175 211 1 212 213 4 184 214 3 183 171 2 215 214 3 125 173 2 217 214 3 127 175 2 219 194 3 194 196 3 196 221 1 222 223 5 0 224 6 0 194 1 224 224 3 226 226 3 206 196 3 208 228 1 229 230 4 227 231 3 226 206 2 232 231 3 196 208 2 234 226 3 218 196 3 220 236 1 237 238 4 227 239 3 226 218 2 240 239 3 196 220 2 242 235 3 235 244 5 0 245 6 0 235 1 245 245 3 247 247 3 243 249 4 248 250 3 247 243 2 251 252 4 246 233 3 253 241 2 254 255 4 225 192 3 256 204 3 253 257 1 258 216 2 259 260 4 182 121 3 261 137 3 256 262 1 263 153 3 253 264 1 265 169 2 266 267 4 109 4 168 161 156 153 1 0.0 20 DATA 1 2.908 1.212 4.552 -5.809 -2.253 1.995 5.681 -5.030 6.775 2.266 8.850 0.099 3.970 8.008 1.302 7.832 6.291 7.219 5.730 9.574 20 DATA 2 1.0 0.5 0.25 0.125 1.0 0.25 0.0625 0.015625 1.0 0.125 0.015625 0.0019531 1.0 0.015625 0.00390625 0.00024414 10.0 1.625 0.4375 0.16308594 20 DATA 3 6.0 4.0 4.0 1.0 4.0 6.0 1.0 4.0 4.0 1.0 6.0 4.0 1.0 4.0 4.0 6.0 30.0 35.0 40.0 45.0 20 DATA 4 4.0 2.0 1.0 3.0 1.0 -5.0 2.0 1.0 1.0 1.0 6.0 2.0 1.0 1.0 1.0 -7.0 13.0 -1.0 27.0 -17.0 20 DATA 5 10.0 3.0 2.0 1.0 3.0 11.0 3.0 2.0 2.0 3.0 12.0 3.0 1.0 2.0 3.0 13.0 26.0 42.0 56.0 66.0 -1 138 1 GRAM-SCHMIDT---DATA 1 1 3 1 2 3 2 101 1 102 3 3 3 103 1 104 4 3 4 105 1 106 107 5 0 1 4 108 2 4 108 3 4 108 4 4 108 5 3 109 6 3 110 113 1 114 7 3 111 115 1 116 8 3 112 117 1 118 119 3 109 119 3 110 119 3 111 119 3 112 5 2 120 6 2 121 7 2 122 8 2 123 124 3 124 125 3 125 128 1 129 126 3 126 130 1 131 127 3 127 132 1 133 134 5 0 124 4 135 125 4 135 126 4 135 127 4 135 9 3 109 10 3 110 140 1 141 11 3 111 142 1 143 12 3 112 144 1 145 146 3 109 146 3 110 146 3 111 146 3 112 9 3 136 10 3 137 151 1 152 11 3 138 153 1 154 12 3 139 155 1 156 157 3 136 147 1 158 157 3 137 148 1 160 157 3 138 149 1 162 157 3 139 150 1 164 9 2 159 10 2 161 11 2 163 12 2 165 166 3 166 167 3 167 170 1 171 168 3 168 172 1 173 169 3 169 174 1 175 176 5 0 166 4 177 167 4 177 168 4 177 169 4 177 13 3 109 14 3 110 182 1 183 15 3 111 184 1 185 16 3 112 186 1 187 188 3 109 188 3 110 188 3 111 188 3 112 13 3 136 14 3 137 193 1 194 15 3 138 195 1 196 16 3 139 197 1 198 199 3 136 189 1 200 199 3 137 190 1 202 199 3 138 191 1 204 199 3 139 192 1 206 13 3 178 14 3 179 208 1 209 15 3 180 210 1 211 16 3 181 212 1 213 214 3 178 201 1 215 214 3 179 203 1 217 214 3 180 205 1 219 214 3 181 207 1 221 13 2 216 14 2 218 15 2 220 16 2 222 223 3 223 224 3 224 227 1 228 225 3 225 229 1 230 226 3 226 231 1 232 233 5 0 223 4 234 224 4 234 225 4 234 226 4 234 16 9 10 11 12 36 37 38 39 78 79 80 81 135 136 137 138 1 0.0 16 DATA 1 2.908 1.212 4.552 -5.809 -2.253 1.995 5.681 -5.030 6.775 2.266 8.850 0.099 3.970 8.008 1.302 7.832 16 DATA 2 1.0 0.5 0.25 0.125 1.0 0.25 0.0625 0.015625 1.0 0.125 0.015625 0.0019531 1.0 0.015625 0.00390625 0.00024414 16 DATA 3 6.0 4.0 4.0 1.0 4.0 6.0 1.0 4.0 4.0 1.0 6.0 4.0 1.0 4.0 4.0 6.0 16 DATA 4 4.0 2.0 1.0 3.0 1.0 -5.0 2.0 1.0 1.0 1.0 6.0 2.0 1.0 1.0 1.0 -7.0 16 DATA 5 10.0 3.0 2.0 1.0 3.0 11.0 3.0 2.0 2.0 3.0 12.0 3.0 1.0 2.0 3.0 13.0 -1 138 1 MODIFIED GRAM-SCHMIDT---DATA 1 1 3 1 2 3 2 101 1 102 3 3 3 103 1 104 4 3 4 105 1 106 107 5 0 1 4 108 2 4 108 3 4 108 4 4 108 109 3 5 110 3 6 113 1 114 111 3 7 115 1 116 112 3 8 117 1 118 119 3 109 5 2 120 119 3 110 6 2 122 119 3 111 7 2 124 119 3 112 8 2 126 109 3 9 110 3 10 128 1 129 111 3 11 130 1 131 112 3 12 132 1 133 134 3 109 9 2 135 134 3 110 10 2 137 134 3 111 11 2 139 134 3 112 12 2 141 109 3 13 110 3 14 143 1 144 111 3 15 145 1 146 112 3 16 147 1 148 149 3 109 13 2 150 149 3 110 14 2 152 149 3 111 15 2 154 149 3 112 16 2 156 121 3 121 123 3 123 158 1 159 125 3 125 160 1 161 127 3 127 162 1 163 164 5 0 121 4 165 123 4 165 125 4 165 127 4 165 166 3 136 167 3 138 170 1 171 168 3 140 172 1 173 169 3 142 174 1 175 176 3 166 136 2 177 176 3 167 138 2 179 176 3 168 140 2 181 176 3 169 142 2 183 166 3 151 167 3 153 185 1 186 168 3 155 187 1 188 169 3 157 189 1 190 191 3 166 151 2 192 191 3 167 153 2 194 191 3 168 155 2 196 191 3 169 157 2 198 178 3 178 180 3 180 200 1 201 182 3 182 202 1 203 184 3 184 204 1 205 206 5 0 178 4 207 180 4 207 182 4 207 184 4 207 208 3 193 209 3 195 212 1 213 210 3 197 214 1 215 211 3 199 216 1 217 218 3 208 193 2 219 218 3 209 195 2 221 218 3 210 197 2 223 218 3 211 199 2 225 220 3 220 222 3 222 227 1 228 224 3 224 229 1 230 226 3 226 231 1 232 233 5 0 220 4 234 222 4 234 224 4 234 226 4 234 16 9 10 11 12 66 67 68 69 108 109 110 111 135 136 137 138 0 16 DATA 1 2.908 1.212 4.552 -5.809 -2.253 1.995 5.681 -5.030 6.775 2.266 8.850 0.099 3.970 8.008 1.302 7.832 16 DATA 2 1.0 0.5 0.25 0.125 1.0 0.25 0.0625 0.015625 1.0 0.125 0.015625 0.0019531 1.0 0.015625 0.00390625 0.00024414 16 DATA 3 6.0 4.0 4.0 1.0 4.0 6.0 1.0 4.0 4.0 1.0 6.0 4.0 1.0 4.0 4.0 6.0 16 DATA 4 4.0 2.0 1.0 3.0 1.0 -5.0 2.0 1.0 1.0 1.0 6.0 2.0 1.0 1.0 1.0 -7.0 16 DATA 5 10.0 3.0 2.0 1.0 3.0 11.0 3.0 2.0 2.0 3.0 12.0 3.0 1.0 2.0 3.0 13.0 -1 172 1 GIVEN'S METHOD---DATA 1 1 3 1 2 3 2 101 1 102 103 5 0 1 4 104 2 4 104 105 3 17 106 3 18 107 1 108 105 3 18 106 3 17 110 2 111 105 3 5 106 3 6 113 1 114 105 3 6 106 3 5 116 2 117 105 3 9 106 3 10 119 1 120 105 3 10 106 3 9 122 2 123 105 3 13 106 3 14 125 1 126 105 3 14 106 3 13 128 2 129 104 3 104 3 3 3 131 1 132 133 5 0 104 4 134 3 4 134 135 3 109 136 3 19 137 1 138 135 3 19 136 3 109 140 2 141 135 3 115 136 3 7 143 1 144 135 3 7 136 3 115 146 2 147 135 3 121 136 3 11 149 1 150 135 3 11 136 3 121 152 2 153 135 3 127 136 3 15 155 1 156 135 3 15 136 3 127 158 2 159 134 3 134 4 3 4 161 1 162 163 5 0 134 4 164 4 4 164 165 3 139 166 3 20 167 1 168 165 3 20 166 3 139 170 2 171 165 3 145 166 3 8 173 1 174 165 3 8 166 3 145 176 2 177 165 3 151 166 3 12 179 1 180 165 3 12 166 3 151 182 2 183 165 3 157 166 3 16 185 1 186 165 3 16 166 3 157 188 2 189 118 3 118 148 3 148 191 1 192 193 5 0 118 4 194 148 4 194 195 3 112 196 3 142 197 1 198 195 3 142 196 3 112 200 2 201 195 3 124 196 3 154 203 1 204 195 3 154 196 3 124 206 2 207 195 3 130 196 3 160 209 1 210 195 3 160 196 3 130 212 2 213 194 3 194 178 3 178 215 1 216 217 5 0 194 4 218 178 4 218 219 3 199 220 3 172 221 1 222 219 3 172 220 3 199 224 2 225 219 3 205 220 3 184 227 1 228 219 3 184 220 3 205 230 2 231 219 3 211 220 3 190 233 1 234 219 3 190 220 3 211 236 2 237 208 3 208 232 3 232 239 1 240 241 5 0 208 4 242 232 4 242 243 3 202 244 3 226 245 1 246 243 3 226 244 3 202 248 2 249 243 3 214 244 3 238 251 1 252 243 3 238 244 3 214 254 2 255 250 4 256 257 3 253 247 2 258 259 4 242 260 3 229 257 3 235 261 1 262 223 2 263 264 4 218 265 3 175 260 3 181 266 1 267 257 3 187 268 1 269 169 2 270 271 4 164 4 172 165 160 157 0 20 DATA 1 2.908 1.212 4.552 -5.809 -2.253 1.995 5.681 -5.030 6.775 2.266 8.850 0.099 3.970 8.008 1.302 7.832 6.291 7.219 5.730 9.574 20 DATA 2 1.0 0.5 0.25 0.125 1.0 0.25 0.0625 0.015625 1.0 0.125 0.015625 0.0019531 1.0 0.015625 0.00390625 0.00024414 10.0 1.625 0.4375 .16308594 20 DATA 3 6.0 4.0 4.0 1.0 4.0 6.0 1.0 4.0 4.0 1.0 6.0 4.0 1.0 4.0 4.0 6.0 30.0 35.0 40.0 45.0 20 DATA 4 4.0 2.0 1.0 3.0 1.0 -5.0 2.0 1.0 1.0 1.0 6.0 2.0 1.0 1.0 1.0 -7.0 13.0 -1.0 27.0 -17.0 20 DATA 5 10.0 3.0 2.0 1.0 3.0 11.0 3.0 2.0 2.0 3.0 12.0 3.0 1.0 2.0 3.0 13.0 26.0 42.0 56.0 66.0 -1 -1