C ALGORITHM 702 , COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 18, NO. 1, MARCH, 1992, PP. 46-70. C Plus remark to appear in 1999. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Doc/ # Doc/Makefile # Drivers/ # Drivers/RES1 # Drivers/RES2 # Drivers/RES3 # Drivers/alg566.f # Drivers/tn1.f # Drivers/tn2.f # Drivers/tn3.f # Src/ # Src/blas1.f # Src/tnpack.f # This archive created: Tue Oct 26 10:22:34 1999 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'Makefile' then echo shar: will not over-write existing file "'Makefile'" else cat << SHAR_EOF > 'Makefile' #---------------------------------------------------------- # Makefile for TNPACK tests #--------------------------------------------------------- F77 = epcf90 FFLAGS = epcf90 -sloppy -C -d5 -g -temp=/tmp -u .SUFFIXES: .f .o .f.o:;$(F77) $(FFLAGS) -c .f test1: tn1.o tnpack.o blas1.o $(F77) $(LFLAGS) -o tn1.o tnpack.o blas1.o test2: tn2.o tnpack.o blas1.o $(F77) $(LFLAGS) -o tn2.o tnpack.o blas1.o test3: tn3.o tnpack.o alg566.o blas1.o $(F77) $(LFLAGS) -o tn3.o tnpack.o alg566.o blas1.o SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'RES1' then echo shar: will not over-write existing file "'RES1'" else cat << SHAR_EOF > 'RES1' PARAMETER INFO, ENTERING TNPACK N = 1000 OPLIST: 1) IPCG = 1 (preconditioning option) 2) IPRINT = 0 (printing option) 3) IPFREQ = 0 (printing frequency) 4) MXITN = 100000 (max Newton itns.) 5) MXITCG = 100 (max PCG itns.) 6) MAXNF = 100000 (max F&G evals.) 7) IORDER = 0 (M-reordering option - calc. per.) 8) IPERPR = 0 (M-reordering printing option) 9) IPKNOW = 0 (M-reordering option - per. known) 10) IHD = 0 (Numeric HD option) 11) MP = 6 (printing unit) 12) IEXIT = 1 (descent option for neg. curvature) 13) ITT = 0 (truncation test option) 14) MAXFEV = 20 (max F evals. in line search) 15) TA = 1 (descent direction test option by PCG) 16) SRLS = 2 (stopping rule for line search) 17) MC = 2 (modified Cholesky) PLIST: 1) EPSF = 1.000E-10 (controls conv. test, F accuracy) 2) EPSG = 1.000E-08 (controls conv. test, G accuracy) 3) ETAR = 1.000E-01 (controls res.-based truncation) 4) ETAQ = 5.000E-01 (controls quad.-model based truncation) 5) PRODK = 1.000E-10 (controls PCG neg curv. test) 6) FTOL = 1.000E-04 (controls F test in line search) 7) GTOL = 9.000E-01 (controls G test in line search) 8) TAU = 1.000E+01 (parameter in UMC) ***** upper-M has 1000 nonzeros ( 0.20 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 1.0242E+05 8.45E+02 1 2.9712E+03 6.14E+01 1 0 1.67E+03 1.00E+00 2 2 2.3331E+03 2.07E+00 1 0 4.63E+01 1.00E+00 3 3 2.3280E+03 6.70E+00 2 2 8.64E+02 1.77E-01 5 4 2.2896E+03 6.86E+00 10 2 2.23E+03 2.49E-03 10 5 2.2343E+03 1.10E+01 26 2 9.71E+03 7.86E-04 15 6 2.0452E+03 2.53E+01 33 2 4.50E+03 4.78E-03 19 7 1.4969E+03 9.50E+00 20 0 1.03E+01 1.00E+00 20 8 1.3740E+03 9.61E+00 41 0 3.66E+00 1.91E-01 22 9 1.1540E+03 1.21E+01 39 0 2.76E+00 5.32E-01 24 10 8.9806E+02 7.10E+00 22 0 3.23E+00 1.00E+00 25 11 7.4635E+02 1.09E+01 26 0 1.96E+00 1.00E+00 26 12 5.4498E+02 5.37E+00 22 0 2.19E+00 1.00E+00 27 13 4.4745E+02 4.94E+00 40 0 1.19E+00 3.20E-01 29 14 3.3057E+02 4.75E+00 35 0 1.01E+00 1.00E+00 30 15 3.2641E+02 1.23E+01 40 0 9.53E-01 1.00E+00 31 16 2.6634E+02 1.17E+01 31 0 1.78E+00 1.00E+00 32 17 1.4484E+02 1.05E+01 25 0 1.95E+00 1.00E+00 33 18 9.2739E+01 7.40E+00 15 0 1.25E+00 3.26E-01 35 19 4.1542E+01 2.61E+00 10 0 8.63E-01 1.00E+00 36 20 2.4740E+01 2.85E+00 16 0 2.20E-01 1.00E+00 37 21 2.1255E+01 4.61E+00 9 0 3.69E-01 1.00E+00 38 22 6.7698E+00 2.50E+00 9 0 4.46E-01 1.00E+00 39 23 4.3075E+00 2.73E+00 9 0 2.12E-01 1.00E+00 40 24 1.7983E-01 4.96E-01 5 0 2.80E-01 1.00E+00 41 25 5.7636E-03 6.80E-02 5 0 2.95E-02 1.00E+00 42 26 5.6278E-05 1.02E-02 3 0 3.11E-03 1.00E+00 43 27 1.3402E-09 5.07E-06 3 0 4.48E-05 1.00E+00 44 28 4.3512E-18 2.82E-09 2 0 5.80E-09 1.00E+00 45 ____________________________ #CG,tot 500 INFORM = 1 CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED CONV. TEST 4 SATISFIED ____________________________ The CPU time in seconds 1.360000 SHAR_EOF fi # end of overwriting check if test -f 'RES2' then echo shar: will not over-write existing file "'RES2'" else cat << SHAR_EOF > 'RES2' PARAMETER INFO, ENTERING TNPACK N = 1000 OPLIST: 1) IPCG = 1 (preconditioning option) 2) IPRINT = 0 (printing option) 3) IPFREQ = 0 (printing frequency) 4) MXITN = 100000 (max Newton itns.) 5) MXITCG = 40 (max PCG itns.) 6) MAXNF = 10000 (max F&G evals.) 7) IORDER = 0 (M-reordering option - calc. per.) 8) IPERPR = 0 (M-reordering printing option) 9) IPKNOW = 0 (M-reordering option - per. known) 10) IHD = 0 (Numeric HD option) 11) MP = 6 (printing unit) 12) IEXIT = 1 (descent option for neg. curvature) 13) ITT = 0 (truncation test option) 14) MAXFEV = 20 (max F evals. in line search) 15) TA = 2 (descent direction test option by PCG) 16) SRLS = 2 (stopping rule for line search) 17) MC = 2 (modified Cholesky) PLIST: 1) EPSF = 1.000E-10 (controls conv. test, F accuracy) 2) EPSG = 1.490E-08 (controls conv. test, G accuracy) 3) ETAR = 5.000E-01 (controls res.-based truncation) 4) ETAQ = 5.000E-01 (controls quad.-model based truncation) 5) PRODK = 1.000E-10 (controls PCG neg curv. test) 6) FTOL = 1.000E-04 (controls F test in line search) 7) GTOL = 9.000E-01 (controls G test in line search) 8) TAU = 5.000E-01 (parameter in UMC) ***** upper-M has 1002 nonzeros ( 0.20 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 2.4882E+05 7.34E+03 1 5.1715E+04 2.35E+03 1 0 2.22E+04 1.00E+00 2 2 1.1001E+04 7.75E+02 1 0 8.46E+03 1.00E+00 3 3 2.3935E+03 2.62E+02 1 0 3.18E+03 1.00E+00 4 4 4.7219E+02 7.75E+01 2 0 6.14E+00 1.00E+00 5 5 9.3201E+01 2.30E+01 2 0 1.85E+00 1.00E+00 6 6 1.8395E+01 6.80E+00 2 0 5.71E-01 1.00E+00 7 7 3.6283E+00 2.01E+00 2 0 1.85E-01 1.00E+00 8 8 7.1451E-01 5.96E-01 2 0 6.57E-02 1.00E+00 9 9 1.4019E-01 1.76E-01 2 0 2.64E-02 1.00E+00 10 10 2.7283E-02 5.20E-02 2 0 1.22E-02 1.00E+00 11 11 5.2219E-03 1.53E-02 2 0 6.30E-03 1.00E+00 12 12 9.6963E-04 4.44E-03 2 0 3.42E-03 1.00E+00 13 13 1.7283E-04 1.26E-03 2 0 1.82E-03 1.00E+00 14 14 1.2755E-04 1.04E-03 5 4 3.22E-03 2.12E-01 16 15 2.5497E-05 3.76E-04 4 4 1.38E-03 1.00E+00 17 16 4.7662E-06 1.35E-04 5 4 2.57E-04 1.00E+00 18 17 6.8607E-07 4.05E-05 6 0 1.54E-05 1.00E+00 19 18 6.2964E-08 9.84E-06 6 0 4.89E-07 1.00E+00 20 19 3.3559E-09 1.98E-06 7 0 3.28E-08 1.00E+00 21 20 6.3769E-11 2.42E-07 8 0 1.38E-09 1.00E+00 22 21 1.1215E-13 9.43E-09 9 5 9.79E-11 1.00E+00 23 ____________________________ #CG,tot 73 INFORM = 1 CONV. TEST 1 SATISFIED CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED CONV. TEST 4 SATISFIED ____________________________ The CPU time in seconds 41.62000 SHAR_EOF fi # end of overwriting check if test -f 'RES3' then echo shar: will not over-write existing file "'RES3'" else cat << SHAR_EOF > 'RES3' ======================= Problem 1 ======================= PARAMETER INFO, ENTERING TNPACK N = 3 OPLIST: 1) IPCG = 1 (preconditioning option) 2) IPRINT = 0 (printing option) 3) IPFREQ = 0 (printing frequency) 4) MXITN = 10000 (max Newton itns.) 5) MXITCG = 100 (max PCG itns.) 6) MAXNF = 4000 (max F&G evals.) 7) IORDER = 0 (M-reordering option - calc. per.) 8) IPERPR = 0 (M-reordering printing option) 9) IPKNOW = 0 (M-reordering option - per. known) 10) IHD = 0 (Numeric HD option) 11) MP = 6 (printing unit) 12) IEXIT = 1 (descent option for neg. curvature) 13) ITT = 0 (truncation test option) 14) MAXFEV = 20 (max F evals. in line search) 15) TA = 2 (descent direction test option by PCG) 16) SRLS = 2 (stopping rule for line search) 17) MC = 2 (modified Cholesky) PLIST: 1) EPSF = 1.000E-10 (controls conv. test, F accuracy) 2) EPSG = 1.000E-08 (controls conv. test, G accuracy) 3) ETAR = 1.000E+00 (controls res.-based truncation) 4) ETAQ = 5.000E-01 (controls quad.-model based truncation) 5) PRODK = 1.000E-10 (controls PCG neg curv. test) 6) FTOL = 1.000E-04 (controls F test in line search) 7) GTOL = 9.000E-01 (controls G test in line search) 8) TAU = 1.000E+01 (parameter in UMC) ***** upper-M has 3 nonzeros ( 50.00 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 2.5000E+03 1.09E+03 1 1.7343E+02 1.74E+02 2 4 2.54E+03 1.00E+00 2 2 5.0479E+01 9.61E+01 1 0 1.37E+02 1.00E+00 3 3 1.3952E+01 4.42E+01 1 0 5.23E+01 1.00E+00 4 4 7.7547E+00 7.68E+00 2 0 8.91E+00 1.00E+00 5 5 7.3822E+00 4.29E+00 3 4 7.51E+00 1.00E+00 6 6 5.2902E+00 7.69E+00 3 0 1.10E-13 3.81E-01 8 7 3.1074E+00 9.95E+00 3 0 4.09E-14 1.00E+00 9 8 1.4921E+00 4.55E+00 3 0 1.36E-13 1.00E+00 10 9 9.1426E-01 9.08E+00 3 0 6.26E-14 1.00E+00 11 10 3.7209E-01 6.16E-01 2 0 1.36E+00 1.00E+00 12 11 1.3083E-01 1.89E+00 3 0 7.32E-15 4.62E-01 14 12 2.0904E-02 9.70E-01 3 0 1.88E-14 1.00E+00 15 13 1.1509E-03 2.15E-01 3 0 4.25E-14 1.00E+00 16 14 5.8748E-06 1.90E-02 3 0 8.85E-16 1.00E+00 17 15 1.8241E-10 9.21E-05 3 0 9.24E-16 1.00E+00 18 16 1.7884E-19 3.35E-09 3 0 4.70E-18 1.00E+00 19 ____________________________ #CG,tot 41 INFORM = 1 CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED CONV. TEST 4 SATISFIED ____________________________ ======================= Problem 2 ======================= ***** upper-M has 6 nonzeros ( 28.57 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 7.7907E-01 1.04E+00 1 5.9650E-01 5.55E-01 1 0 1.31E+00 1.00E+00 2 2 3.1288E-01 1.87E-01 2 0 9.76E-02 1.00E+00 3 3 2.9374E-01 2.30E-02 2 0 3.48E-02 1.00E+00 4 4 2.9004E-01 4.30E-02 3 4 1.08E-01 1.00E+00 5 5 2.8541E-01 2.24E-02 3 4 4.75E-02 1.00E+00 6 6 2.8130E-01 4.13E-02 3 4 1.04E-01 1.00E+00 7 7 2.7664E-01 2.62E-02 3 4 5.68E-02 1.00E+00 8 8 2.7210E-01 4.05E-02 3 4 1.02E-01 1.00E+00 9 9 2.6716E-01 3.04E-02 3 4 6.65E-02 1.00E+00 10 10 2.6216E-01 3.99E-02 3 4 1.01E-01 1.00E+00 11 11 2.5685E-01 3.47E-02 3 4 7.63E-02 1.00E+00 12 12 2.5135E-01 3.91E-02 3 4 9.89E-02 1.00E+00 13 13 2.4558E-01 3.92E-02 3 4 8.65E-02 1.00E+00 14 14 2.3954E-01 3.82E-02 3 4 9.66E-02 1.00E+00 15 15 2.3320E-01 4.39E-02 3 4 9.73E-02 1.00E+00 16 16 2.2655E-01 3.71E-02 3 4 9.38E-02 1.00E+00 17 17 2.1950E-01 4.91E-02 3 4 1.09E-01 1.00E+00 18 18 2.1213E-01 3.57E-02 3 4 9.01E-02 1.00E+00 19 19 2.0413E-01 5.49E-02 3 4 1.22E-01 1.00E+00 20 20 1.9586E-01 3.40E-02 3 4 8.55E-02 1.00E+00 21 21 1.3853E-01 7.22E-02 4 0 6.60E-12 1.86E-02 23 22 9.2647E-02 3.07E-01 3 0 3.60E-04 1.00E+00 24 23 4.8813E-02 1.04E-02 2 0 1.87E-02 1.00E+00 25 24 4.3142E-02 4.16E-02 4 4 1.74E-02 4.35E-01 27 25 3.6266E-02 8.87E-02 4 0 9.64E-13 5.16E-01 29 26 2.5315E-02 6.02E-02 4 0 1.17E-10 1.00E+00 30 27 1.6396E-02 9.74E-02 4 0 1.36E-10 1.00E+00 31 28 1.3678E-02 2.79E-03 1 0 5.87E-03 1.00E+00 32 29 9.9249E-03 9.59E-02 4 0 1.11E-12 1.00E+00 33 30 6.7620E-03 2.10E-03 3 0 1.19E-03 1.00E+00 34 31 5.8356E-03 2.12E-03 4 0 4.25E-11 1.00E+00 35 32 5.6760E-03 1.50E-03 4 0 1.83E-09 1.00E+00 36 33 5.6563E-03 3.05E-04 4 0 9.30E-08 1.00E+00 37 34 5.6556E-03 1.32E-05 5 4 6.76E-06 1.00E+00 38 35 5.6556E-03 3.89E-06 3 4 9.54E-06 1.00E+00 39 36 5.6556E-03 1.36E-05 3 4 3.34E-05 1.00E+00 40 37 5.6556E-03 4.06E-06 3 4 9.95E-06 1.00E+00 41 38 5.6556E-03 1.44E-05 3 4 3.53E-05 1.00E+00 42 39 5.6556E-03 4.30E-06 3 4 1.05E-05 1.00E+00 43 40 5.6556E-03 1.53E-05 3 4 3.76E-05 1.00E+00 44 41 5.6556E-03 4.59E-06 3 4 1.12E-05 1.00E+00 45 42 5.6556E-03 1.64E-05 3 4 4.03E-05 1.00E+00 46 43 5.6556E-03 4.93E-06 3 4 1.21E-05 1.00E+00 47 44 5.6556E-03 1.77E-05 3 4 4.34E-05 1.00E+00 48 45 5.6556E-03 5.31E-06 3 4 1.30E-05 1.00E+00 49 46 5.6556E-03 1.91E-05 3 4 4.67E-05 1.00E+00 50 47 5.6556E-03 5.73E-06 3 4 1.40E-05 1.00E+00 51 48 5.6556E-03 2.06E-05 3 4 5.05E-05 1.00E+00 52 49 5.6556E-03 6.18E-06 3 4 1.51E-05 1.00E+00 53 50 5.6556E-03 2.22E-05 3 4 5.45E-05 1.00E+00 54 51 5.6556E-03 6.68E-06 3 4 1.64E-05 1.00E+00 55 52 5.6556E-03 2.40E-05 3 4 5.89E-05 1.00E+00 56 53 5.6556E-03 7.22E-06 3 4 1.77E-05 1.00E+00 57 54 5.6556E-03 2.60E-05 3 4 6.36E-05 1.00E+00 58 55 5.6556E-03 7.81E-06 3 4 1.91E-05 1.00E+00 59 56 5.6556E-03 2.80E-05 3 4 6.87E-05 1.00E+00 60 57 5.6556E-03 8.44E-06 3 4 2.07E-05 1.00E+00 61 58 5.6556E-03 3.03E-05 3 4 7.42E-05 1.00E+00 62 59 5.6556E-03 9.13E-06 3 4 2.24E-05 1.00E+00 63 60 5.6556E-03 3.27E-05 3 4 8.01E-05 1.00E+00 64 61 5.6556E-03 9.87E-06 3 4 2.42E-05 1.00E+00 65 62 5.6556E-03 3.53E-05 3 4 8.65E-05 1.00E+00 66 63 5.6556E-03 1.07E-05 3 4 2.61E-05 1.00E+00 67 64 5.6556E-03 3.82E-05 3 4 9.34E-05 1.00E+00 68 65 5.6556E-03 1.15E-05 3 4 2.83E-05 1.00E+00 69 66 5.6556E-03 4.12E-05 3 4 1.01E-04 1.00E+00 70 67 5.6556E-03 1.25E-05 3 4 3.06E-05 1.00E+00 71 68 5.6556E-03 4.45E-05 3 4 1.09E-04 1.00E+00 72 69 5.6556E-03 1.35E-05 3 4 3.30E-05 1.00E+00 73 70 5.6556E-03 4.80E-05 3 4 1.18E-04 1.00E+00 74 71 5.6556E-03 1.46E-05 3 4 3.57E-05 1.00E+00 75 72 5.6556E-03 5.18E-05 3 4 1.27E-04 1.00E+00 76 73 5.6556E-03 1.58E-05 3 4 3.86E-05 1.00E+00 77 74 5.6556E-03 5.59E-05 3 4 1.37E-04 1.00E+00 78 75 5.6556E-03 1.70E-05 3 4 4.17E-05 1.00E+00 79 76 5.6556E-03 6.03E-05 3 4 1.48E-04 1.00E+00 80 77 5.6556E-03 1.84E-05 3 4 4.51E-05 1.00E+00 81 78 5.6556E-03 6.50E-05 3 4 1.59E-04 1.00E+00 82 79 5.6556E-03 1.99E-05 3 4 4.88E-05 1.00E+00 83 80 5.6556E-03 7.00E-05 3 4 1.72E-04 1.00E+00 84 81 5.6556E-03 2.15E-05 3 4 5.27E-05 1.00E+00 85 82 5.6556E-03 7.55E-05 3 4 1.85E-04 1.00E+00 86 83 5.6555E-03 2.33E-05 3 4 5.70E-05 1.00E+00 87 84 5.6555E-03 8.13E-05 3 4 1.99E-04 1.00E+00 88 85 5.6555E-03 2.51E-05 3 4 6.16E-05 1.00E+00 89 86 5.6555E-03 8.76E-05 3 4 2.14E-04 1.00E+00 90 87 5.6555E-03 2.72E-05 3 4 6.65E-05 1.00E+00 91 88 5.6555E-03 9.42E-05 3 4 2.31E-04 1.00E+00 92 89 5.6555E-03 2.93E-05 3 4 7.19E-05 1.00E+00 93 90 5.6555E-03 1.01E-04 3 4 2.48E-04 1.00E+00 94 91 5.6555E-03 3.17E-05 3 4 7.76E-05 1.00E+00 95 92 5.6554E-03 1.09E-04 3 4 2.67E-04 1.00E+00 96 93 5.6554E-03 3.42E-05 3 4 8.38E-05 1.00E+00 97 94 5.6554E-03 1.17E-04 3 4 2.87E-04 1.00E+00 98 95 5.6554E-03 3.70E-05 3 4 9.05E-05 1.00E+00 99 96 5.6554E-03 1.26E-04 3 4 3.08E-04 1.00E+00 100 97 5.6554E-03 3.99E-05 3 4 9.78E-05 1.00E+00 101 98 5.6553E-03 1.35E-04 3 4 3.31E-04 1.00E+00 102 99 5.6553E-03 4.31E-05 3 4 1.06E-04 1.00E+00 103 100 5.6553E-03 1.45E-04 3 4 3.55E-04 1.00E+00 104 101 5.6552E-03 4.65E-05 3 4 1.14E-04 1.00E+00 105 102 5.6552E-03 1.56E-04 3 4 3.81E-04 1.00E+00 106 103 5.6552E-03 5.02E-05 3 4 1.23E-04 1.00E+00 107 104 5.6551E-03 1.67E-04 3 4 4.08E-04 1.00E+00 108 105 5.6551E-03 5.42E-05 3 4 1.33E-04 1.00E+00 109 106 5.6551E-03 1.78E-04 3 4 4.37E-04 1.00E+00 110 107 5.6550E-03 5.84E-05 3 4 1.43E-04 1.00E+00 111 108 5.6550E-03 1.91E-04 3 4 4.67E-04 1.00E+00 112 109 5.6549E-03 6.30E-05 3 4 1.54E-04 1.00E+00 113 110 5.6549E-03 2.04E-04 3 4 4.99E-04 1.00E+00 114 111 5.6548E-03 6.79E-05 3 4 1.66E-04 1.00E+00 115 112 5.6547E-03 2.18E-04 3 4 5.33E-04 1.00E+00 116 113 5.6547E-03 7.32E-05 3 4 1.79E-04 1.00E+00 117 114 5.6546E-03 2.32E-04 3 4 5.68E-04 1.00E+00 118 115 5.6545E-03 7.89E-05 3 4 1.93E-04 1.00E+00 119 116 5.6544E-03 2.47E-04 3 4 6.05E-04 1.00E+00 120 117 5.6543E-03 8.50E-05 3 4 2.08E-04 1.00E+00 121 118 5.6542E-03 2.63E-04 3 4 6.44E-04 1.00E+00 122 119 5.6541E-03 9.16E-05 3 4 2.24E-04 1.00E+00 123 120 5.6540E-03 2.79E-04 3 4 6.84E-04 1.00E+00 124 121 5.6539E-03 9.86E-05 3 4 2.41E-04 1.00E+00 125 122 5.6537E-03 2.96E-04 3 4 7.25E-04 1.00E+00 126 123 5.6536E-03 1.06E-04 3 4 2.60E-04 1.00E+00 127 124 5.6534E-03 3.14E-04 3 4 7.68E-04 1.00E+00 128 125 5.6532E-03 1.14E-04 3 4 2.79E-04 1.00E+00 129 126 5.6531E-03 3.32E-04 3 4 8.12E-04 1.00E+00 130 127 5.6529E-03 1.23E-04 3 4 3.01E-04 1.00E+00 131 128 5.6527E-03 3.50E-04 3 4 8.57E-04 1.00E+00 132 129 5.6525E-03 1.32E-04 3 4 3.23E-04 1.00E+00 133 130 5.6522E-03 3.69E-04 3 4 9.03E-04 1.00E+00 134 131 5.6520E-03 1.42E-04 3 4 3.47E-04 1.00E+00 135 132 5.6517E-03 3.88E-04 3 4 9.50E-04 1.00E+00 136 133 5.6514E-03 1.52E-04 3 4 3.73E-04 1.00E+00 137 134 5.6511E-03 4.08E-04 3 4 9.98E-04 1.00E+00 138 135 5.6508E-03 1.63E-04 3 4 4.00E-04 1.00E+00 139 136 5.6504E-03 4.27E-04 3 4 1.05E-03 1.00E+00 140 137 5.6500E-03 1.75E-04 3 4 4.29E-04 1.00E+00 141 138 5.6496E-03 4.47E-04 3 4 1.09E-03 1.00E+00 142 139 5.6492E-03 1.88E-04 3 4 4.60E-04 1.00E+00 143 140 5.6487E-03 4.67E-04 3 4 1.14E-03 1.00E+00 144 141 5.6483E-03 2.01E-04 3 4 4.93E-04 1.00E+00 145 142 5.6477E-03 4.87E-04 3 4 1.19E-03 1.00E+00 146 143 5.6472E-03 2.16E-04 3 4 5.28E-04 1.00E+00 147 144 5.6466E-03 5.07E-04 3 4 1.24E-03 1.00E+00 148 145 5.6460E-03 2.31E-04 3 4 5.66E-04 1.00E+00 149 146 5.6453E-03 5.27E-04 3 4 1.29E-03 1.00E+00 150 147 5.6446E-03 2.47E-04 3 4 6.05E-04 1.00E+00 151 148 5.6438E-03 5.46E-04 3 4 1.34E-03 1.00E+00 152 149 5.6430E-03 2.64E-04 3 4 6.48E-04 1.00E+00 153 150 5.6422E-03 5.66E-04 3 4 1.38E-03 1.00E+00 154 151 5.6412E-03 2.83E-04 3 4 6.92E-04 1.00E+00 155 152 5.6403E-03 5.86E-04 3 4 1.43E-03 1.00E+00 156 153 5.6392E-03 3.02E-04 3 4 7.40E-04 1.00E+00 157 154 5.6381E-03 6.06E-04 3 4 1.48E-03 1.00E+00 158 155 5.6370E-03 3.22E-04 3 4 7.90E-04 1.00E+00 159 156 5.6357E-03 6.27E-04 3 4 1.53E-03 1.00E+00 160 157 5.6344E-03 3.44E-04 3 4 8.43E-04 1.00E+00 161 158 5.6330E-03 6.47E-04 3 4 1.58E-03 1.00E+00 162 159 5.6316E-03 3.67E-04 3 4 8.99E-04 1.00E+00 163 160 5.6300E-03 6.69E-04 3 4 1.63E-03 1.00E+00 164 161 5.6283E-03 3.91E-04 3 4 9.58E-04 1.00E+00 165 162 5.6266E-03 6.90E-04 3 4 1.69E-03 1.00E+00 166 163 5.6247E-03 4.17E-04 3 4 1.02E-03 1.00E+00 167 164 5.6227E-03 7.13E-04 3 4 1.74E-03 1.00E+00 168 165 5.6207E-03 4.44E-04 3 4 1.09E-03 1.00E+00 169 166 5.6184E-03 7.36E-04 3 4 1.80E-03 1.00E+00 170 167 5.6161E-03 4.72E-04 3 4 1.16E-03 1.00E+00 171 168 5.6136E-03 7.60E-04 3 4 1.86E-03 1.00E+00 172 169 5.6110E-03 5.03E-04 3 4 1.23E-03 1.00E+00 173 170 5.6083E-03 7.85E-04 3 4 1.92E-03 1.00E+00 174 171 5.6053E-03 5.34E-04 3 4 1.31E-03 1.00E+00 175 172 5.6023E-03 8.11E-04 3 4 1.98E-03 1.00E+00 176 173 5.5990E-03 5.68E-04 3 4 1.39E-03 1.00E+00 177 174 5.5956E-03 8.39E-04 3 4 2.05E-03 1.00E+00 178 175 5.5920E-03 6.03E-04 3 4 1.48E-03 1.00E+00 179 176 5.5881E-03 8.67E-04 3 4 2.12E-03 1.00E+00 180 177 5.5841E-03 6.40E-04 3 4 1.57E-03 1.00E+00 181 178 5.5799E-03 8.97E-04 3 4 2.19E-03 1.00E+00 182 179 5.5754E-03 6.78E-04 3 4 1.66E-03 1.00E+00 183 180 5.5707E-03 9.27E-04 3 4 2.27E-03 1.00E+00 184 181 5.5658E-03 7.19E-04 3 4 1.76E-03 1.00E+00 185 182 5.5606E-03 9.59E-04 3 4 2.34E-03 1.00E+00 186 183 5.5552E-03 7.61E-04 3 4 1.86E-03 1.00E+00 187 184 5.5495E-03 9.91E-04 3 4 2.42E-03 1.00E+00 188 185 5.5435E-03 8.05E-04 3 4 1.97E-03 1.00E+00 189 186 5.5372E-03 1.02E-03 3 4 2.50E-03 1.00E+00 190 187 5.5306E-03 8.51E-04 3 4 2.08E-03 1.00E+00 191 188 5.5237E-03 1.06E-03 3 4 2.58E-03 1.00E+00 192 189 5.5165E-03 8.99E-04 3 4 2.20E-03 1.00E+00 193 190 5.5089E-03 1.09E-03 3 4 2.67E-03 1.00E+00 194 191 5.5010E-03 9.49E-04 3 4 2.32E-03 1.00E+00 195 192 5.4928E-03 1.13E-03 3 4 2.75E-03 1.00E+00 196 193 5.4842E-03 1.00E-03 3 4 2.45E-03 1.00E+00 197 194 5.4752E-03 1.16E-03 3 4 2.83E-03 1.00E+00 198 195 5.4658E-03 1.05E-03 3 4 2.58E-03 1.00E+00 199 196 5.4561E-03 1.19E-03 3 4 2.91E-03 1.00E+00 200 197 5.4459E-03 1.11E-03 3 4 2.71E-03 1.00E+00 201 198 5.4354E-03 1.22E-03 3 4 2.99E-03 1.00E+00 202 199 5.4244E-03 1.17E-03 3 4 2.85E-03 1.00E+00 203 200 5.4130E-03 1.26E-03 3 4 3.07E-03 1.00E+00 204 201 5.4013E-03 1.22E-03 3 4 2.99E-03 1.00E+00 205 202 5.3890E-03 1.29E-03 3 4 3.15E-03 1.00E+00 206 203 5.3764E-03 1.28E-03 3 4 3.13E-03 1.00E+00 207 204 5.3633E-03 1.32E-03 3 4 3.22E-03 1.00E+00 208 205 5.3497E-03 1.35E-03 3 4 3.28E-03 1.00E+00 209 206 5.3357E-03 1.34E-03 3 4 3.29E-03 1.00E+00 210 207 5.3213E-03 1.41E-03 3 4 3.43E-03 1.00E+00 211 208 5.3064E-03 1.37E-03 3 4 3.35E-03 1.00E+00 212 209 5.2911E-03 1.47E-03 3 4 3.58E-03 1.00E+00 213 210 5.2753E-03 1.39E-03 3 4 3.41E-03 1.00E+00 214 211 5.2591E-03 1.54E-03 3 4 3.73E-03 1.00E+00 215 212 5.2424E-03 1.41E-03 3 4 3.46E-03 1.00E+00 216 213 5.2254E-03 1.60E-03 3 4 3.89E-03 1.00E+00 217 214 5.2078E-03 1.43E-03 3 4 3.51E-03 1.00E+00 218 215 5.1899E-03 1.67E-03 3 4 4.05E-03 1.00E+00 219 216 5.1714E-03 1.45E-03 3 4 3.55E-03 1.00E+00 220 217 5.1527E-03 1.74E-03 3 4 4.21E-03 1.00E+00 221 218 5.1334E-03 1.46E-03 3 4 3.58E-03 1.00E+00 222 219 5.1138E-03 1.80E-03 3 4 4.36E-03 1.00E+00 223 220 5.0938E-03 1.47E-03 3 4 3.60E-03 1.00E+00 224 221 5.0735E-03 1.87E-03 3 4 4.52E-03 1.00E+00 225 222 5.0526E-03 1.48E-03 3 4 3.62E-03 1.00E+00 226 223 5.0316E-03 1.94E-03 3 4 4.68E-03 1.00E+00 227 224 5.0101E-03 1.48E-03 3 4 3.63E-03 1.00E+00 228 225 4.9884E-03 2.00E-03 3 4 4.83E-03 1.00E+00 229 226 4.9662E-03 1.48E-03 3 4 3.64E-03 1.00E+00 230 227 4.9439E-03 2.07E-03 3 4 4.98E-03 1.00E+00 231 228 4.9094E-03 3.58E-02 4 4 3.44E-02 5.04E-02 234 229 3.4231E-03 4.68E-03 5 4 4.90E-03 1.00E+00 235 230 2.8494E-03 1.62E-02 6 4 6.89E-04 3.40E-02 239 231 2.5941E-03 4.71E-03 5 4 1.21E-03 1.00E+00 240 232 2.3752E-03 6.28E-03 5 4 1.88E-03 5.22E-01 242 233 2.2062E-03 4.26E-03 5 4 1.88E-03 5.07E-01 244 234 1.8153E-03 1.79E-02 6 4 4.83E-04 5.99E-02 247 235 1.4097E-03 6.50E-03 6 4 9.00E-04 1.00E+00 248 236 8.8479E-04 9.70E-03 6 4 4.65E-03 1.11E-01 250 237 6.4422E-04 5.81E-03 5 4 9.12E-04 1.00E+00 251 238 5.1053E-04 2.20E-02 5 0 4.09E-05 2.16E-01 253 239 2.6532E-04 1.86E-02 7 0 1.52E-10 4.25E-01 255 240 1.0474E-04 3.82E-03 7 4 3.00E-04 1.00E+00 256 241 6.0857E-05 6.14E-03 7 4 2.92E-04 1.00E+00 257 242 4.5866E-05 5.54E-05 4 0 5.40E-05 1.00E+00 258 243 4.0753E-05 3.94E-03 7 4 3.85E-05 1.38E-01 260 244 3.0314E-05 2.39E-03 5 0 7.43E-06 1.00E+00 261 245 2.6151E-05 2.84E-03 7 0 2.14E-09 4.85E-01 263 246 2.2027E-05 1.18E-03 5 0 2.35E-05 1.00E+00 264 247 1.8321E-05 2.29E-03 7 0 1.35E-09 5.16E-01 266 248 1.4092E-05 1.79E-03 7 0 2.49E-11 1.00E+00 267 249 1.0650E-05 2.04E-03 7 4 1.33E-04 1.00E+00 268 250 8.3516E-06 5.21E-04 5 0 1.65E-05 1.00E+00 269 251 6.6844E-06 1.09E-03 7 4 1.61E-04 2.87E-01 271 252 5.0034E-06 1.87E-03 7 4 1.86E-04 1.00E+00 272 253 3.9939E-06 5.72E-06 4 0 1.27E-05 1.00E+00 273 254 2.9832E-06 8.05E-04 8 0 1.53E-13 1.86E-01 275 255 2.0259E-06 8.25E-04 7 4 1.77E-04 1.00E+00 276 256 1.3230E-06 4.58E-04 7 0 6.79E-09 1.00E+00 277 257 8.9452E-07 7.14E-04 7 0 2.14E-09 1.00E+00 278 258 5.2317E-07 1.97E-04 7 0 1.41E-09 1.00E+00 279 259 3.9583E-07 7.29E-04 7 4 1.78E-04 1.00E+00 280 260 2.0033E-07 9.47E-07 4 0 2.23E-06 1.00E+00 281 261 1.2384E-07 1.35E-04 8 5 1.80E-10 2.51E-01 283 262 6.1715E-08 2.52E-04 7 4 1.03E-04 1.00E+00 284 263 1.9829E-08 4.07E-05 7 0 2.72E-09 1.00E+00 285 264 7.8033E-09 1.31E-04 7 0 1.30E-10 1.00E+00 286 265 1.5731E-09 3.02E-07 6 4 3.34E-07 1.00E+00 287 266 1.4754E-10 2.05E-05 7 4 5.61E-07 1.00E+00 288 267 1.9696E-13 2.88E-09 5 5 7.05E-09 1.00E+00 289 ____________________________ #CG,tot 930 INFORM = 1 CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED CONV. TEST 4 SATISFIED ____________________________ ======================= Problem 3 ======================= ***** upper-M has 3 nonzeros ( 50.00 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 3.8881E-06 4.30E-03 1 1.2943E-08 1.25E-05 1 0 2.37E-05 1.00E+00 2 2 1.1279E-08 6.74E-09 2 0 8.23E-21 1.00E+00 3 ____________________________ #CG,tot 3 INFORM = 1 CONV. TEST 3 SATISFIED CONV. TEST 4 SATISFIED ____________________________ ======================= Problem 4 ======================= ***** upper-M has 2 nonzeros ( 66.67 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 1.1353E+00 1.41E+04 1 1.3484E-01 6.03E+01 1 0 3.09E+01 1.00E+00 2 2 1.2233E-01 5.85E+03 2 0 3.64E-12 1.00E+00 3 3 4.7900E-02 1.19E+01 1 0 8.35E+00 1.00E+00 4 4 2.9613E-02 3.14E+03 2 0 3.30E-12 1.00E+00 5 5 1.7518E-02 1.49E+00 1 0 1.40E+00 1.00E+00 6 6 1.0189E-02 2.18E+03 2 0 3.64E-12 1.00E+00 7 7 6.4346E-03 2.92E-01 1 0 3.01E-01 1.00E+00 8 8 3.9072E-03 1.68E+03 2 0 1.36E-12 1.00E+00 9 9 2.3638E-03 6.92E-02 1 0 7.24E-02 1.00E+00 10 10 1.6146E-03 1.36E+03 2 0 1.14E-13 1.00E+00 11 11 8.6782E-04 1.87E-02 1 0 1.91E-02 1.00E+00 12 12 7.2182E-04 1.14E+03 2 0 1.68E-12 1.00E+00 13 13 3.1823E-04 5.70E-03 1 0 5.50E-03 1.00E+00 14 14 2.0693E-04 2.38E+02 2 0 1.14E-13 5.06E-01 16 15 1.9140E-04 1.92E-03 1 0 2.63E-03 1.00E+00 17 16 1.2922E-04 1.76E+02 2 0 2.59E-12 4.52E-01 19 17 1.2152E-04 1.08E-03 1 0 1.49E-03 1.00E+00 20 18 8.4814E-05 1.34E+02 2 0 9.58E-12 4.06E-01 22 19 8.0734E-05 6.53E-04 1 0 9.05E-04 1.00E+00 23 20 5.7986E-05 1.05E+02 2 0 3.35E-12 3.69E-01 25 21 5.5675E-05 4.14E-04 1 0 5.76E-04 1.00E+00 26 22 4.0979E-05 8.42E+01 2 0 5.24E-12 3.38E-01 28 23 3.9588E-05 2.73E-04 1 0 3.82E-04 1.00E+00 29 24 2.9751E-05 6.92E+01 2 0 4.19E-12 3.13E-01 31 25 2.8867E-05 1.87E-04 1 0 2.61E-04 1.00E+00 32 26 2.2080E-05 5.81E+01 2 0 2.16E-12 2.92E-01 34 27 2.1490E-05 1.31E-04 1 0 1.84E-04 1.00E+00 35 28 1.6686E-05 4.97E+01 2 0 2.47E-12 2.74E-01 37 29 1.6276E-05 9.42E-05 1 0 1.32E-04 1.00E+00 38 30 1.2801E-05 4.32E+01 2 0 1.06E-11 2.60E-01 40 31 1.2505E-05 6.89E-05 1 0 9.68E-05 1.00E+00 41 32 9.9438E-06 3.81E+01 2 0 6.10E-12 2.47E-01 43 33 9.7234E-06 5.13E-05 1 0 7.20E-05 1.00E+00 44 34 7.8061E-06 3.40E+01 2 0 1.50E-11 2.37E-01 46 35 7.6372E-06 3.86E-05 1 0 5.43E-05 1.00E+00 47 36 7.6372E-06 1.28E-05 1 1 5.47E-05 1.70E-10 52 ____________________________ #CG,tot 53 INFORM = 1 CONV. TEST 1 SATISFIED CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED ____________________________ ======================= Problem 5 ======================= ***** upper-M has 3 nonzeros ( 50.00 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 1.0312E+03 8.62E+01 1 4.6036E+02 7.51E+01 1 1 1.49E+02 2.61E-02 7 2 5.2247E+01 7.09E+01 2 0 1.57E-02 1.00E+00 8 3 2.2510E+01 2.29E+01 1 0 5.47E+00 1.00E+00 9 4 1.2861E+01 6.59E+00 1 0 8.72E+00 1.00E+00 10 5 2.4015E+00 4.38E+00 2 0 4.62E-02 1.00E+00 11 6 7.3565E-01 1.53E+00 2 0 3.35E-02 1.00E+00 12 7 1.9620E-01 5.81E-01 2 0 1.68E-02 1.00E+00 13 8 4.0810E-02 2.09E-01 2 0 7.19E-03 1.00E+00 14 9 5.3958E-03 6.55E-02 2 0 1.91E-03 1.00E+00 15 10 6.1868E-04 1.29E-02 2 0 3.31E-04 1.00E+00 16 11 2.6148E-05 2.07E-03 3 0 2.28E-14 1.00E+00 17 12 8.7855E-08 1.38E-04 3 0 2.05E-14 1.00E+00 18 13 1.4722E-12 4.13E-07 3 0 1.39E-15 1.00E+00 19 14 5.6077E-13 1.85E-08 3 5 3.20E-08 1.00E+00 20 ____________________________ #CG,tot 29 INFORM = 1 CONV. TEST 1 SATISFIED CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED ____________________________ ======================= Problem 6 ======================= ***** upper-M has 3 nonzeros ( 50.00 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 4.9760E+02 9.00E+02 1 1.0118E+02 2.68E+02 1 0 1.24E+00 1.00E+00 2 2 2.1651E+01 8.09E+01 1 0 2.01E+00 1.00E+00 3 3 5.4612E+00 2.51E+01 1 0 2.38E+00 1.00E+00 4 4 2.0545E+00 8.21E+00 1 0 2.37E+00 1.00E+00 5 5 1.2456E+00 2.94E+00 1 0 2.08E+00 1.00E+00 6 6 7.7702E-02 7.85E-01 2 0 4.66E-01 1.00E+00 7 7 1.0223E-03 7.16E-02 2 0 5.97E-02 1.00E+00 8 8 1.2181E-05 4.10E-03 2 0 7.09E-03 1.00E+00 9 9 3.2357E-22 8.04E-11 3 0 2.42E-19 1.00E+00 10 ____________________________ #CG,tot 14 INFORM = 1 CONV. TEST 3 SATISFIED CONV. TEST 4 SATISFIED ____________________________ ======================= Problem 7 ======================= ***** upper-M has 3 nonzeros ( 50.00 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 3.0000E+01 4.90E+01 1 7.2941E+00 1.76E+01 1 0 3.12E+01 1.00E+00 2 2 1.1688E+00 3.19E+00 1 0 6.15E+00 1.00E+00 3 3 9.2465E-01 1.01E+00 1 0 1.79E+00 1.00E+00 4 4 5.2021E-01 1.23E+00 2 0 1.79E-01 1.00E+00 5 5 4.8428E-01 2.47E-01 1 0 4.00E-01 1.00E+00 6 6 4.7157E-01 8.06E-02 3 0 6.21E-17 1.00E+00 7 7 4.7142E-01 8.64E-03 1 0 1.50E-02 1.00E+00 8 8 4.7140E-01 1.46E-04 3 0 2.00E-17 1.00E+00 9 9 4.7140E-01 1.19E-09 3 0 2.74E-19 1.00E+00 10 ____________________________ #CG,tot 16 INFORM = 1 CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED CONV. TEST 4 SATISFIED ____________________________ ======================= Problem 8 ======================= ***** upper-M has 3 nonzeros ( 50.00 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 1.8906E+02 1.19E+02 1 3.8740E+01 3.66E+01 1 0 1.54E+01 1.00E+00 2 2 7.8067E+00 1.13E+01 1 0 5.01E+00 1.00E+00 3 3 1.4805E+00 3.40E+00 1 0 1.42E+00 1.00E+00 4 4 2.4173E-01 9.78E-01 1 0 2.98E-01 1.00E+00 5 5 2.8774E-02 2.54E-01 1 0 4.36E-02 1.00E+00 6 6 1.6588E-03 5.04E-02 1 0 4.56E-03 1.00E+00 7 7 3.5152E-05 4.51E-03 1 0 2.65E-04 1.00E+00 8 8 2.0131E-05 5.31E-05 1 0 2.59E-05 1.00E+00 9 9 1.9834E-05 4.49E-04 3 0 3.16E-13 1.75E-01 11 10 1.9395E-05 2.44E-04 2 0 2.15E-07 1.00E+00 12 11 1.9199E-05 7.30E-04 3 0 3.52E-14 1.00E+00 13 12 1.8633E-05 8.01E-05 2 0 2.41E-07 1.00E+00 14 13 1.8416E-05 2.51E-04 3 0 3.08E-13 2.03E-01 16 14 1.8168E-05 5.41E-04 3 0 2.64E-13 1.00E+00 17 15 1.7778E-05 1.10E-04 2 0 2.40E-07 1.00E+00 18 16 1.7605E-05 2.28E-04 3 0 1.40E-13 2.73E-01 20 17 1.7400E-05 4.95E-04 3 0 1.94E-13 1.00E+00 21 18 1.7076E-05 9.83E-05 2 0 2.14E-07 1.00E+00 22 19 1.6935E-05 2.05E-04 3 0 9.06E-14 2.74E-01 24 20 1.6767E-05 4.41E-04 3 0 1.78E-13 1.00E+00 25 21 1.6509E-05 8.81E-05 2 0 1.74E-07 1.00E+00 26 22 1.6396E-05 1.83E-04 3 0 2.90E-13 2.82E-01 28 23 1.6261E-05 3.79E-04 3 0 3.43E-13 1.00E+00 29 24 1.6065E-05 8.03E-05 2 0 1.30E-07 1.00E+00 30 25 1.5978E-05 1.63E-04 3 0 5.06E-14 3.03E-01 32 26 1.5870E-05 3.07E-04 3 0 5.21E-14 1.00E+00 33 27 1.5732E-05 7.69E-05 2 0 8.84E-08 1.00E+00 34 28 1.5666E-05 1.49E-04 3 0 1.20E-13 3.58E-01 36 29 1.5581E-05 2.18E-04 3 0 1.00E-13 1.00E+00 37 30 1.5493E-05 8.61E-05 2 0 5.46E-08 1.00E+00 38 31 1.5482E-05 3.65E-04 3 5 4.62E-08 1.00E+00 39 32 1.5362E-05 1.88E-05 2 0 3.08E-08 1.00E+00 40 33 1.5331E-05 8.41E-05 3 5 2.77E-08 2.05E-01 42 34 1.5299E-05 1.80E-04 3 5 2.22E-08 1.00E+00 43 35 1.5257E-05 3.19E-05 2 0 1.45E-08 1.00E+00 44 36 1.5240E-05 7.06E-05 3 5 1.17E-08 3.62E-01 46 37 1.5221E-05 9.28E-05 2 0 8.68E-09 1.00E+00 47 38 1.5204E-05 3.64E-05 2 0 5.26E-09 1.00E+00 48 39 1.5196E-05 9.86E-05 2 0 3.55E-09 1.00E+00 49 40 1.5186E-05 9.32E-06 2 0 1.47E-09 1.00E+00 50 41 1.5183E-05 2.85E-05 3 5 1.02E-09 4.74E-01 52 42 1.5181E-05 1.76E-05 2 0 4.63E-10 1.00E+00 53 43 1.5180E-05 1.01E-05 2 0 1.60E-10 1.00E+00 54 44 1.5179E-05 3.44E-06 2 0 3.47E-11 1.00E+00 55 45 1.5179E-05 3.49E-08 2 5 6.05E-08 1.00E+00 56 ____________________________ #CG,tot 101 INFORM = 1 CONV. TEST 1 SATISFIED CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED ____________________________ ======================= Problem 9 ======================= ***** upper-M has 3 nonzeros ( 50.00 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 3.4000E-01 2.45E+00 1 5.2485E-02 4.94E-01 1 0 2.26E-01 1.00E+00 2 2 3.5429E-02 1.35E-01 1 0 2.15E-01 1.00E+00 3 3 1.0907E-02 1.56E-01 2 0 2.54E-02 2.22E-01 5 4 1.1279E-03 1.07E-01 2 0 7.51E-04 1.00E+00 6 5 2.1022E-05 4.81E-03 1 0 6.14E-03 1.00E+00 7 6 3.2039E-06 1.92E-04 2 0 1.91E-08 1.00E+00 8 7 3.2001E-06 4.93E-08 2 0 7.96E-08 1.00E+00 9 8 3.2000E-06 1.14E-05 3 0 1.09E-15 3.53E-02 12 9 3.2000E-06 4.46E-08 3 5 7.73E-08 1.00E+00 13 ____________________________ #CG,tot 17 INFORM = 1 CONV. TEST 1 SATISFIED CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED ____________________________ ======================= Problem 10 ======================= ***** upper-M has 2 nonzeros ( 66.67 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 1.0000E+12 1.41E+06 1 5.0000E+11 3.54E+11 1 0 0.00E+00 1.00E+00 2 2 2.5000E+11 7.07E+05 1 1 5.00E+11 2.00E-12 7 3 3.9999E+00 2.83E+06 2 0 8.02E-12 1.00E+00 8 4 1.9722E-31 6.28E-10 1 1 4.00E+06 5.00E-13 14 ____________________________ #CG,tot 5 INFORM = 1 CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED CONV. TEST 4 SATISFIED ____________________________ ======================= Problem 11 ======================= ***** upper-M has 4 nonzeros ( 40.00 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 7.9267E+06 1.07E+06 1 2.4743E+06 2.87E+05 1 0 1.88E+05 1.00E+00 2 2 1.5973E+06 1.07E+05 1 0 2.27E+05 1.00E+00 3 3 5.7661E+05 3.12E+04 2 0 5.71E+04 1.00E+00 4 4 2.3078E+05 1.21E+04 3 0 3.95E+03 1.00E+00 5 5 1.0245E+05 6.17E+03 4 0 3.82E-11 1.00E+00 6 6 8.6899E+04 1.22E+03 4 0 6.38E-11 1.00E+00 7 7 8.5833E+04 1.10E+02 3 0 1.55E+02 1.00E+00 8 8 8.5822E+04 1.66E+00 4 0 5.06E-14 1.00E+00 9 9 8.5822E+04 8.97E-02 2 0 1.79E-01 1.00E+00 10 10 8.5822E+04 8.96E-03 3 0 1.79E-02 1.00E+00 11 ____________________________ #CG,tot 27 INFORM = 1 CONV. TEST 1 SATISFIED CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED ____________________________ ======================= Problem 12 ======================= ***** upper-M has 3 nonzeros ( 50.00 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 1.2111E+01 2.29E+01 1 7.7272E+00 1.35E+01 1 0 3.90E-01 1.58E-01 3 2 6.7366E+00 4.34E+00 1 0 2.83E-01 1.00E+00 4 3 6.6386E+00 2.11E-01 1 0 5.63E-01 1.00E+00 5 4 6.2797E+00 3.42E+00 3 4 1.17E-01 3.54E-01 7 5 6.2198E+00 1.47E-01 1 0 3.74E-01 1.00E+00 8 6 5.7640E+00 3.72E+00 2 4 1.77E-01 2.87E-01 10 7 5.6955E+00 8.68E-02 1 0 2.49E-01 1.00E+00 11 8 1.1051E+00 5.97E+00 3 4 2.27E+00 1.21E-01 13 9 7.2737E-01 4.79E-01 1 0 3.57E-01 1.00E+00 14 10 3.7606E-01 6.32E-01 2 0 3.65E-02 9.87E-02 16 11 3.4757E-01 1.73E+00 2 0 9.49E-03 1.00E+00 17 12 1.7056E-02 2.71E-01 2 0 2.64E-03 1.00E+00 18 13 5.3882E-03 4.03E-02 2 0 6.23E-04 1.00E+00 19 14 5.3614E-03 1.16E-03 1 0 2.04E-03 1.00E+00 20 15 3.4708E-03 1.25E-01 3 0 2.93E-13 4.20E-01 22 16 3.2185E-03 4.98E-04 1 0 5.84E-04 1.00E+00 23 17 1.9786E-03 9.37E-02 3 0 4.84E-13 3.78E-01 25 18 1.8366E-03 9.26E-04 1 0 1.44E-03 1.00E+00 26 19 1.0227E-03 8.23E-02 3 0 7.33E-13 4.58E-01 28 20 9.1320E-04 1.03E-03 1 0 1.66E-03 1.00E+00 29 21 4.4429E-04 6.55E-02 3 0 1.34E-13 5.22E-01 31 22 3.7502E-04 8.87E-04 1 0 1.46E-03 1.00E+00 32 23 2.7177E-04 1.23E-01 3 0 4.43E-13 1.00E+00 33 24 3.0445E-05 8.28E-04 1 0 1.15E-03 1.00E+00 34 25 3.1637E-06 1.32E-02 3 0 3.38E-12 1.00E+00 35 26 3.9823E-07 7.99E-05 1 0 1.36E-04 1.00E+00 36 27 8.3256E-10 2.15E-04 3 0 9.70E-14 1.00E+00 37 28 7.9992E-11 2.70E-08 2 0 4.67E-08 1.00E+00 38 29 7.9990E-11 1.40E-07 1 1 4.67E-08 1.00E+00 39 ____________________________ #CG,tot 53 INFORM = 1 CONV. TEST 1 SATISFIED CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED ____________________________ ======================= Problem 13 ======================= ***** upper-M has 3 nonzeros ( 50.00 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 1.4165E-02 7.40E-02 1 6.5610E-03 3.19E-02 2 4 8.71E-02 3.66E-01 3 2 6.0759E-03 5.66E-02 2 0 2.18E-03 1.00E+00 4 3 3.4019E-03 9.75E-03 2 0 4.92E-03 1.00E+00 5 4 3.2145E-03 5.41E-03 3 4 9.44E-03 1.00E+00 6 5 3.1654E-03 2.44E-02 3 4 2.53E-02 1.00E+00 7 6 2.5854E-03 2.91E-03 3 0 1.37E-16 1.00E+00 8 7 2.5737E-03 1.38E-04 3 0 4.76E-18 1.00E+00 9 8 2.5737E-03 1.55E-07 3 0 7.83E-20 1.00E+00 10 9 2.5737E-03 1.10E-12 3 0 4.86E-22 1.00E+00 11 ____________________________ #CG,tot 24 INFORM = 1 CONV. TEST 1 SATISFIED CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED CONV. TEST 4 SATISFIED ____________________________ ======================= Problem 14 ======================= ***** upper-M has 2 nonzeros ( 66.67 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 2.4200E+01 1.65E+02 1 4.5032E+00 9.10E+00 1 0 6.06E+00 1.00E+00 2 2 4.4456E+00 1.45E+00 1 0 1.98E+00 1.00E+00 3 3 4.3695E+00 1.08E+01 2 4 1.51E+01 1.00E+00 4 4 4.2819E+00 1.48E+00 1 0 2.00E+00 1.00E+00 5 5 4.1548E+00 1.58E+01 2 4 2.11E+01 1.00E+00 6 6 3.9288E+00 1.57E+00 1 0 1.96E+00 1.00E+00 7 7 3.2481E+00 7.46E+00 2 4 1.46E+02 3.40E-02 10 8 2.7941E+00 1.03E+01 2 0 4.58E-14 4.43E-01 12 9 2.2043E+00 8.55E+00 2 0 8.66E-14 1.00E+00 13 10 1.6873E+00 5.60E+00 2 0 1.07E-14 1.00E+00 14 11 1.2928E+00 5.65E+00 2 0 9.92E-15 1.00E+00 15 12 9.2551E-01 2.62E+00 2 0 5.09E-15 1.00E+00 16 13 7.5306E-01 7.09E+00 2 0 2.29E-16 1.00E+00 17 14 4.3837E-01 8.38E-01 2 0 1.15E-14 1.00E+00 18 15 3.3226E-01 2.09E+00 2 0 6.56E-15 3.01E-01 20 16 2.2060E-01 4.66E+00 2 0 3.92E-14 1.00E+00 21 17 1.1558E-01 8.97E-01 2 0 7.20E-14 1.00E+00 22 18 1.0216E-01 7.17E+00 2 0 5.95E-14 1.00E+00 23 19 3.2292E-02 1.34E-01 1 0 2.61E-01 1.00E+00 24 20 1.5244E-02 1.21E+00 2 0 4.25E-15 4.16E-01 26 21 4.2622E-03 9.82E-01 2 0 1.37E-13 1.00E+00 27 22 6.2667E-04 3.02E-01 2 0 8.16E-14 1.00E+00 28 23 2.7993E-05 1.05E-01 2 0 4.20E-14 1.00E+00 29 24 1.6491E-05 3.23E-03 1 0 4.58E-03 1.00E+00 30 25 2.7010E-08 5.18E-03 2 0 2.63E-15 1.00E+00 31 26 1.5213E-10 9.65E-06 1 0 1.37E-05 1.00E+00 32 27 2.3172E-18 4.80E-08 2 0 9.47E-22 1.00E+00 33 28 1.3433E-20 9.08E-11 2 5 1.28E-10 1.00E+00 34 ____________________________ #CG,tot 49 INFORM = 1 CONV. TEST 1 SATISFIED CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED CONV. TEST 4 SATISFIED ____________________________ ======================= Problem 15 ======================= ***** upper-M has 4 nonzeros ( 40.00 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 2.1500E+02 2.29E+02 1 4.8587E+01 6.34E+01 1 0 9.04E+01 1.00E+00 2 2 1.6032E+01 1.88E+01 1 0 4.37E+01 1.00E+00 3 3 9.0118E+00 6.53E+00 2 0 1.03E+01 1.00E+00 4 4 8.0751E-01 3.38E+00 4 0 4.43E-15 1.00E+00 5 5 1.5951E-01 1.00E+00 4 0 1.73E-15 1.00E+00 6 6 3.1919E-02 3.27E-01 3 0 1.82E-01 1.00E+00 7 7 6.2831E-03 8.84E-02 4 0 2.85E-16 1.00E+00 8 8 1.2411E-03 2.62E-02 4 0 3.47E-17 1.00E+00 9 9 2.4516E-04 7.76E-03 4 0 2.44E-16 1.00E+00 10 10 4.8426E-05 2.30E-03 4 0 3.66E-16 1.00E+00 11 11 9.5656E-06 6.81E-04 4 0 2.74E-16 1.00E+00 12 12 1.8895E-06 2.02E-04 4 0 1.13E-15 1.00E+00 13 13 3.7323E-07 5.98E-05 4 0 2.23E-14 1.00E+00 14 14 7.3725E-08 1.77E-05 4 0 1.62E-14 1.00E+00 15 15 1.4563E-08 5.25E-06 4 0 1.42E-13 1.00E+00 16 16 2.8767E-09 1.56E-06 4 0 8.99E-13 1.00E+00 17 17 5.6823E-10 4.61E-07 4 0 8.79E-13 1.00E+00 18 18 1.1224E-10 1.37E-07 4 0 5.76E-13 1.00E+00 19 19 2.2171E-11 4.05E-08 5 5 3.66E-12 1.00E+00 20 20 8.7411E-12 1.49E-08 4 5 3.52E-08 1.00E+00 21 21 3.0425E-12 1.05E-08 4 5 1.57E-08 1.00E+00 22 22 1.4061E-12 3.34E-09 4 5 7.95E-09 1.00E+00 23 ____________________________ #CG,tot 80 INFORM = 1 CONV. TEST 1 SATISFIED CONV. TEST 3 SATISFIED CONV. TEST 4 SATISFIED ____________________________ ======================= Problem 16 ======================= ***** upper-M has 2 nonzeros ( 66.67 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 1.4203E+01 1.96E+01 1 7.1404E+00 8.14E+00 1 0 1.12E+01 1.00E+00 2 2 2.9940E+00 7.58E+00 2 0 1.08E-14 2.85E-01 4 3 5.5221E-01 5.32E+00 1 0 2.34E+00 1.00E+00 5 4 9.1321E-03 4.36E-01 1 0 6.90E-01 1.00E+00 6 5 3.8615E-03 4.88E-02 1 0 7.86E-02 1.00E+00 7 6 1.7094E-04 5.74E-02 2 0 2.78E-17 1.00E+00 8 7 5.4757E-07 5.87E-04 2 0 1.99E-16 1.00E+00 9 8 8.4770E-12 1.13E-05 2 0 4.77E-18 1.00E+00 10 9 2.0461E-21 5.97E-11 2 0 4.42E-20 1.00E+00 11 ____________________________ #CG,tot 14 INFORM = 1 CONV. TEST 1 SATISFIED CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED CONV. TEST 4 SATISFIED ____________________________ ======================= Problem 17 ======================= ***** upper-M has 4 nonzeros ( 40.00 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 1.9192E+04 8.20E+03 1 8.9658E+02 1.04E+03 1 0 1.30E+03 1.00E+00 2 2 5.5820E+02 2.31E+02 1 0 3.25E+02 1.00E+00 3 3 1.0083E+02 2.64E+02 1 0 1.01E+02 1.00E+00 4 4 4.0143E+01 4.33E+01 1 0 5.27E+01 1.00E+00 5 5 1.1109E+01 3.60E+01 2 0 4.17E+00 1.00E+00 6 6 9.0119E+00 5.30E+00 1 0 9.02E+00 1.00E+00 7 7 7.8878E+00 2.33E+00 2 0 1.26E-01 1.00E+00 8 8 7.8763E+00 2.83E-01 1 0 5.56E-01 1.00E+00 9 9 7.8719E+00 1.86E-02 2 0 4.28E-02 1.00E+00 10 10 7.8719E+00 6.36E-02 2 4 1.27E-01 1.00E+00 11 11 7.8698E+00 6.60E-01 3 4 1.70E+00 1.00E+00 12 12 7.8687E+00 2.83E-02 1 0 5.59E-02 1.00E+00 13 13 7.8683E+00 1.91E-01 4 4 3.98E-01 1.00E+00 14 14 7.8679E+00 3.04E-02 4 4 5.97E-02 1.00E+00 15 15 7.8676E+00 1.75E-01 4 4 3.65E-01 1.00E+00 16 16 7.8672E+00 3.28E-02 4 4 6.42E-02 1.00E+00 17 17 7.8668E+00 1.63E-01 4 4 3.41E-01 1.00E+00 18 18 7.8664E+00 3.55E-02 4 4 6.95E-02 1.00E+00 19 19 7.8659E+00 1.54E-01 4 4 3.22E-01 1.00E+00 20 20 7.8655E+00 3.88E-02 4 4 7.57E-02 1.00E+00 21 21 7.8651E+00 1.46E-01 4 4 3.07E-01 1.00E+00 22 22 7.8646E+00 4.28E-02 4 4 8.34E-02 1.00E+00 23 23 7.8641E+00 1.40E-01 4 4 2.94E-01 1.00E+00 24 24 7.8637E+00 4.79E-02 4 4 9.32E-02 1.00E+00 25 25 7.8632E+00 1.35E-01 4 4 2.84E-01 1.00E+00 26 26 7.8627E+00 5.46E-02 4 4 1.06E-01 1.00E+00 27 27 7.8621E+00 1.32E-01 4 4 2.76E-01 1.00E+00 28 28 7.8616E+00 6.36E-02 4 4 1.24E-01 1.00E+00 29 29 7.8610E+00 1.29E-01 4 4 2.70E-01 1.00E+00 30 30 7.8604E+00 7.65E-02 4 4 1.49E-01 1.00E+00 31 31 7.8598E+00 1.27E-01 4 4 2.65E-01 1.00E+00 32 32 7.8592E+00 9.56E-02 4 4 1.86E-01 1.00E+00 33 33 7.8585E+00 1.26E-01 4 4 2.63E-01 1.00E+00 34 34 7.8578E+00 1.27E-01 4 4 2.47E-01 1.00E+00 35 35 7.8571E+00 1.27E-01 4 4 2.62E-01 1.00E+00 36 36 7.8563E+00 1.84E-01 4 4 3.59E-01 1.00E+00 37 37 7.8555E+00 1.29E-01 4 4 2.65E-01 1.00E+00 38 38 7.8546E+00 3.19E-01 4 4 6.24E-01 1.00E+00 39 39 7.8536E+00 1.34E-01 4 4 2.72E-01 1.00E+00 40 40 7.8508E+00 1.16E+00 4 4 2.21E+00 1.00E+00 41 41 7.8473E+00 1.52E-01 4 4 3.06E-01 1.00E+00 42 42 7.8462E+00 8.88E-02 3 4 1.81E-01 1.00E+00 43 43 7.8450E+00 1.59E-01 3 4 3.17E-01 1.00E+00 44 44 7.8433E+00 5.18E-01 4 4 1.05E+00 1.00E+00 45 45 7.8415E+00 1.73E-01 4 4 3.41E-01 1.00E+00 46 46 7.8399E+00 3.04E-01 4 4 6.19E-01 1.00E+00 47 47 7.8382E+00 1.89E-01 4 4 3.67E-01 1.00E+00 48 48 7.8364E+00 2.28E-01 4 4 4.67E-01 1.00E+00 49 49 7.8346E+00 2.12E-01 4 4 4.04E-01 1.00E+00 50 50 7.8326E+00 1.90E-01 4 4 3.93E-01 1.00E+00 51 51 7.8305E+00 2.42E-01 4 4 4.55E-01 1.00E+00 52 52 7.8283E+00 1.69E-01 4 4 3.53E-01 1.00E+00 53 53 7.8259E+00 2.83E-01 4 4 5.28E-01 1.00E+00 54 54 7.8233E+00 1.57E-01 4 4 3.30E-01 1.00E+00 55 55 7.8204E+00 3.45E-01 4 4 6.37E-01 1.00E+00 56 56 7.8172E+00 1.51E-01 4 4 3.20E-01 1.00E+00 57 57 7.8137E+00 4.43E-01 4 4 8.13E-01 1.00E+00 58 58 7.8096E+00 1.52E-01 4 4 3.20E-01 1.00E+00 59 59 7.8049E+00 6.26E-01 4 4 1.14E+00 1.00E+00 60 60 7.7989E+00 1.61E-01 4 4 3.36E-01 1.00E+00 61 61 7.7912E+00 1.14E+00 4 4 1.99E+00 1.00E+00 62 62 7.7785E+00 2.02E-01 4 4 3.95E-01 1.00E+00 63 63 7.7707E+00 4.05E-01 3 4 7.59E-01 1.00E+00 64 64 7.7597E+00 5.82E-01 4 4 1.29E+00 1.00E+00 65 65 7.7466E+00 5.69E-01 4 4 9.85E-01 1.00E+00 66 66 7.7275E+00 3.90E-01 4 4 9.26E-01 1.00E+00 67 67 7.7026E+00 9.81E-01 4 4 1.42E+00 1.00E+00 68 68 7.6382E+00 1.46E+00 4 4 1.54E+00 1.00E+00 69 69 7.5058E+00 2.20E+00 4 4 2.41E+00 1.00E+00 70 70 7.3476E+00 3.39E+00 4 0 9.48E-13 8.65E-02 73 71 7.0869E+00 5.89E+00 4 0 4.54E-14 1.00E+00 74 72 6.4024E+00 1.54E+00 4 0 8.29E-14 1.00E+00 75 73 6.0092E+00 4.97E+00 4 0 1.37E-13 1.80E-01 77 74 5.2592E+00 6.56E+00 4 0 2.12E-13 1.00E+00 78 75 4.4698E+00 8.45E+00 4 0 7.52E-13 1.00E+00 79 76 3.7490E+00 8.22E+00 4 0 5.46E-13 1.00E+00 80 77 3.1975E+00 7.22E+00 4 0 4.99E-13 1.00E+00 81 78 2.7604E+00 5.08E+00 4 0 5.32E-14 1.00E+00 82 79 2.4070E+00 4.85E+00 4 0 2.54E-13 1.00E+00 83 80 2.0488E+00 3.01E+00 4 0 1.03E-13 1.00E+00 84 81 1.7916E+00 4.91E+00 4 0 5.62E-14 1.00E+00 85 82 1.3635E+00 1.36E+00 4 0 8.64E-15 1.00E+00 86 83 1.1632E+00 2.21E+00 4 0 1.04E-14 2.55E-01 88 84 9.2601E-01 6.16E+00 4 0 4.69E-14 1.00E+00 89 85 5.1642E-01 9.75E-01 4 0 4.21E-14 1.00E+00 90 86 3.5861E-01 2.45E+00 4 0 7.40E-15 3.01E-01 92 87 1.8559E-01 4.77E+00 4 0 2.52E-13 1.00E+00 93 88 5.5997E-02 9.98E-01 4 0 4.95E-13 1.00E+00 94 89 1.7691E-02 2.59E+00 4 0 9.13E-14 1.00E+00 95 90 7.6297E-04 1.09E-01 4 0 1.78E-13 1.00E+00 96 91 8.0259E-06 5.93E-02 4 0 6.73E-14 1.00E+00 97 92 2.6383E-10 7.04E-05 4 0 5.68E-15 1.00E+00 98 93 1.0847E-18 2.15E-08 4 0 6.92E-17 1.00E+00 99 94 1.5576E-19 1.11E-09 2 5 2.23E-09 1.00E+00 100 ____________________________ #CG,tot 341 INFORM = 1 CONV. TEST 1 SATISFIED CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED CONV. TEST 4 SATISFIED ____________________________ ======================= Problem 18 ======================= ***** upper-M has 4 nonzeros ( 40.00 % ) ***** ____________________________ ITN. F |G| #CG mode |R| LAMBDA NFG 0 7.1184E-02 4.56E-01 1 3.7489E-02 2.84E-01 1 0 4.26E-01 1.00E+00 2 2 2.8140E-02 7.04E-01 1 1 5.69E-01 2.33E-01 4 3 8.6116E-04 7.31E-02 1 0 4.11E-01 1.00E+00 5 4 5.0772E-05 1.51E-02 2 0 3.76E-16 1.00E+00 6 5 6.0931E-08 5.86E-04 2 0 1.75E-16 1.00E+00 7 6 1.0716E-13 7.12E-07 2 0 2.08E-16 1.00E+00 8 7 3.3521E-25 1.33E-12 2 0 5.40E-17 1.00E+00 9 ____________________________ #CG,tot 11 INFORM = 1 CONV. TEST 1 SATISFIED CONV. TEST 2 SATISFIED CONV. TEST 3 SATISFIED CONV. TEST 4 SATISFIED ____________________________ SHAR_EOF fi # end of overwriting check if test -f 'alg566.f' then echo shar: will not over-write existing file "'alg566.f'" else cat << SHAR_EOF > 'alg566.f' C ===== 4. DOUBLE PRECISION TESTING AIDS FOR UNCONSTRAINED NONLINEAR C ===== OPTIMIZATION. SUBROUTINE INITPT(N,X,NPROB,FACTOR) C ********** C C SUBROUTINE INITPT C C THIS SUBROUTINE SPECIFIES THE STANDARD STARTING POINTS FOR THE C FUNCTIONS DEFINED BY SUBROUTINE OBJFCN. THE SUBROUTINE RETURNS C IN X A MULTIPLE (FACTOR) OF THE STANDARD STARTING POINT. FOR C THE SEVENTH FUNCTION THE STANDARD STARTING POINT IS ZERO, SO IN C THIS CASE, IF FACTOR IS NOT UNITY, THEN THE SUBROUTINE RETURNS C THE VECTOR X(J) = FACTOR, J=1,...,N. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE INITPT(N,X,NPROB,FACTOR) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE. C C X IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE STANDARD C STARTING POINT FOR PROBLEM NPROB MULTIPLIED BY FACTOR. C C NPROB IS A POSITIVE INTEGER INPUT VARIABLE WHICH DEFINES THE C NUMBER OF THE PROBLEM. NPROB MUST NOT EXCEED 18. C C FACTOR IS AN INPUT VARIABLE WHICH SPECIFIES THE MULTIPLE OF C THE STANDARD STARTING POINT. IF FACTOR IS UNITY, NO C MULTIPLICATION IS PERFORMED. C C MINPACK. VERSION OF JULY 1978. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** C .. Scalar Arguments .. DOUBLE PRECISION FACTOR INTEGER N,NPROB C .. C .. Array Arguments .. DOUBLE PRECISION X(N) C .. C .. Local Scalars .. DOUBLE PRECISION C1,C2,C3,C4,FIVE,H,HALF,ONE,TEN,THREE,TWENTY, + TWNTF,TWO,ZERO INTEGER IVAR,J C .. C .. Statement Functions .. DOUBLE PRECISION DFLOAT C .. C .. Data statements .. DATA ZERO,HALF,ONE,TWO,THREE,FIVE,TEN,TWENTY,TWNTF/0.0D0,0.5D0, + 1.0D0,2.0D0,3.0D0,5.0D0,1.0D1,2.0D1,2.5D1/ DATA C1,C2,C3,C4/4.0D-1,2.5D0,1.5D-1,1.2D0/ C .. C .. Statement Function definitions .. DFLOAT(IVAR) = IVAR C .. C C SELECTION OF INITIAL POINT. C GO TO (10,20,30,40,50,60,80,100,120, + 140,150,160,170,190,210,230,240, + 250) NPROB C C HELICAL VALLEY FUNCTION. C 10 CONTINUE X(1) = -ONE X(2) = ZERO X(3) = ZERO GO TO 270 C C BIGGS EXP6 FUNCTION. C 20 CONTINUE X(1) = ONE X(2) = TWO X(3) = ONE X(4) = ONE X(5) = ONE X(6) = ONE c do i = 1, 6 c X(i) = 5.0d0 c enddo GO TO 270 C C GAUSSIAN FUNCTION. C 30 CONTINUE X(1) = C1 X(2) = ONE X(3) = ZERO GO TO 270 C C POWELL BADLY SCALED FUNCTION. C 40 CONTINUE X(1) = ZERO X(2) = ONE GO TO 270 C C BOX 3-DIMENSIONAL FUNCTION. C 50 CONTINUE X(1) = ZERO X(2) = TEN X(3) = TWENTY GO TO 270 C C VARIABLY DIMENSIONED FUNCTION. C 60 CONTINUE H = ONE/DFLOAT(N) DO 70 J = 1,N X(J) = ONE - DFLOAT(J)*H 70 CONTINUE GO TO 270 C C WATSON FUNCTION. C 80 CONTINUE DO 90 J = 1,N X(J) = ZERO 90 CONTINUE GO TO 270 C C PENALTY FUNCTION I. C 100 CONTINUE DO 110 J = 1,N X(J) = DFLOAT(J) 110 CONTINUE GO TO 270 C C PENALTY FUNCTION II. C 120 CONTINUE DO 130 J = 1,N X(J) = HALF 130 CONTINUE GO TO 270 C C BROWN BADLY SCALED FUNCTION. C 140 CONTINUE X(1) = ONE X(2) = ONE GO TO 270 C C BROWN AND DENNIS FUNCTION. C 150 CONTINUE X(1) = TWNTF X(2) = FIVE X(3) = -FIVE X(4) = -ONE GO TO 270 C C GULF RESEARCH AND DEVELOPMENT FUNCTION. C 160 CONTINUE X(1) = FIVE X(2) = C2 X(3) = C3 GO TO 270 C C TRIGONOMETRIC FUNCTION. C 170 CONTINUE H = ONE/DFLOAT(N) DO 180 J = 1,N X(J) = H 180 CONTINUE GO TO 270 C C EXTENDED ROSENBROCK FUNCTION. C 190 CONTINUE DO 200 J = 1,N,2 X(J) = -C4 X(J+1) = ONE 200 CONTINUE GO TO 270 C C EXTENDED POWELL SINGULAR FUNCTION. C 210 CONTINUE DO 220 J = 1,N,4 X(J) = THREE X(J+1) = -ONE X(J+2) = ZERO X(J+3) = ONE 220 CONTINUE GO TO 270 C C BEALE FUNCTION. C 230 CONTINUE X(1) = ONE X(2) = ONE GO TO 270 C C WOOD FUNCTION. C 240 CONTINUE X(1) = -THREE X(2) = -ONE X(3) = -THREE X(4) = -ONE GO TO 270 C C CHEBYQUAD FUNCTION. C 250 CONTINUE H = ONE/DFLOAT(N+1) DO 260 J = 1,N X(J) = DFLOAT(J)*H 260 CONTINUE 270 CONTINUE C C COMPUTE MULTIPLE OF INITIAL POINT. C IF (FACTOR.EQ.ONE) GO TO 320 IF (NPROB.EQ.7) GO TO 290 DO 280 J = 1,N X(J) = FACTOR*X(J) 280 CONTINUE GO TO 310 290 CONTINUE DO 300 J = 1,N X(J) = FACTOR 300 CONTINUE 310 CONTINUE 320 CONTINUE RETURN C C LAST CARD OF SUBROUTINE INITPT. C END SUBROUTINE OBJFCN(N,X,F,NPROB) C ********** C C SUBROUTINE OBJFCN C C THIS SUBROUTINE DEFINES THE OBJECTIVE FUNCTIONS OF EIGHTEEN C NONLINEAR UNCONSTRAINED MINIMIZATION PROBLEMS. THE VALUES C OF N FOR FUNCTIONS 1,2,3,4,5,10,11,12,16 AND 17 ARE C 3,6,3,2,3,2,4,3,2 AND 4, RESPECTIVELY. C FOR FUNCTION 7, N MAY BE 2 OR GREATER BUT IS USUALLY 6 OR 9. C FOR FUNCTIONS 6,8,9,13,14,15 AND 18 N MAY BE VARIABLE, C HOWEVER IT MUST BE EVEN FOR FUNCTION 14, A MULTIPLE OF 4 FOR C FUNCTION 15, AND NOT GREATER THAN 50 FOR FUNCTION 18. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE OBJFCN(N,X,F,NPROB) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE. C C X IS AN INPUT ARRAY OF LENGTH N. C C F IS AN OUTPUT VARIABLE WHICH CONTAINS THE VALUE OF C THE NPROB OBJECTIVE FUNCTION EVALUATED AT X. C C NPROB IS A POSITIVE INTEGER INPUT VARIABLE WHICH DEFINES THE C NUMBER OF THE PROBLEM. NPROB MUST NOT EXCEED 18. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... DABS,DATAN,DCOS,DEXP,DLOG,DSIGN,DSIN, C DSQRT C C MINPACK. VERSION OF JULY 1978. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** C .. Scalar Arguments .. DOUBLE PRECISION F INTEGER N,NPROB C .. C .. Array Arguments .. DOUBLE PRECISION X(N) C .. C .. Local Scalars .. DOUBLE PRECISION AP,ARG,BP,C100,C10000,C1P5,C1PD6,C25,C29,C2P25, + C2P625,C2PDM6,C3P5,C90,CP0001,CP1,CP2,CP25,CP5, + D1,D2,EIGHT,FIFTY,FIVE,FOUR,ONE,R,S1,S2,S3,T,T1, + T2,T3,TEN,TH,THREE,TPI,TWO,ZERO INTEGER I,IEV,IVAR,J C .. C .. Local Arrays .. DOUBLE PRECISION FVEC(50),Y(15) C .. C .. Intrinsic Functions .. INTRINSIC DABS,DATAN,DCOS,DEXP,DLOG,DSIGN,DSIN,DSQRT C .. C .. Statement Functions .. DOUBLE PRECISION DFLOAT C .. C .. Data statements .. DATA ZERO,ONE,TWO,THREE,FOUR,FIVE,EIGHT,TEN,FIFTY/0.0D0,1.0D0, + 2.0D0,3.0D0,4.0D0,5.0D0,8.0D0,1.0D1,5.0D1/ DATA C2PDM6,CP0001,CP1,CP2,CP25,CP5,C1P5,C2P25,C2P625,C3P5,C25, + C29,C90,C100,C10000,C1PD6/2.0D-6,1.0D-4,1.0D-1,2.0D-1,2.5D-1, + 5.0D-1,1.5D0,2.25D0,2.625D0,3.5D0,2.5D1,2.9D1,9.0D1,1.0D2, + 1.0D4,1.0D6/ DATA AP,BP/1.0D-5,1.0D0/ DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8),Y(9),Y(10),Y(11), + Y(12),Y(13),Y(14),Y(15)/9.0D-4,4.4D-3,1.75D-2,5.4D-2, + 1.295D-1,2.42D-1,3.521D-1,3.989D-1,3.521D-1,2.42D-1,1.295D-1, + 5.4D-2,1.75D-2,4.4D-3,9.0D-4/ C .. C .. Statement Function definitions .. DFLOAT(IVAR) = IVAR C .. C C FUNCTION ROUTINE SELECTOR. C GO TO (10,20,40,60,70,90,110,150,170, + 200,210,230,250,280,300,320,330, + 340) NPROB C C HELICAL VALLEY FUNCTION. C 10 CONTINUE TPI = EIGHT*DATAN(ONE) TH = DSIGN(CP25,X(2)) IF (X(1).GT.ZERO) TH = DATAN(X(2)/X(1))/TPI IF (X(1).LT.ZERO) TH = DATAN(X(2)/X(1))/TPI + CP5 ARG = X(1)**2 + X(2)**2 R = DSQRT(ARG) T = X(3) - TEN*TH F = C100* (T**2+ (R-ONE)**2) + X(3)**2 GO TO 390 C C BIGGS EXP6 FUNCTION. C 20 CONTINUE F = ZERO DO 30 I = 1,13 c DO 30 I = 1, 6 D1 = DFLOAT(I)/TEN D2 = DEXP(-D1) - FIVE*DEXP(-TEN*D1) + THREE*DEXP(-FOUR*D1) S1 = DEXP(-D1*X(1)) S2 = DEXP(-D1*X(2)) S3 = DEXP(-D1*X(5)) T = X(3)*S1 - X(4)*S2 + X(6)*S3 - D2 F = F + T**2 30 CONTINUE GO TO 390 C C GAUSSIAN FUNCTION. C 40 CONTINUE F = ZERO DO 50 I = 1,15 D1 = CP5*DFLOAT(I-1) D2 = C3P5 - D1 - X(3) ARG = -CP5*X(2)*D2**2 R = DEXP(ARG) T = X(1)*R - Y(I) F = F + T**2 50 CONTINUE GO TO 390 C C POWELL BADLY SCALED FUNCTION. C 60 CONTINUE T1 = C10000*X(1)*X(2) - ONE S1 = DEXP(-X(1)) S2 = DEXP(-X(2)) T2 = S1 + S2 - ONE - CP0001 F = T1**2 + T2**2 GO TO 390 C C BOX 3-DIMENSIONAL FUNCTION. C 70 CONTINUE F = ZERO DO 80 I = 1,10 D1 = DFLOAT(I) D2 = D1/TEN S1 = DEXP(-D2*X(1)) S2 = DEXP(-D2*X(2)) S3 = DEXP(-D2) - DEXP(-D1) T = S1 - S2 - S3*X(3) F = F + T**2 80 CONTINUE GO TO 390 C C VARIABLY DIMENSIONED FUNCTION. C 90 CONTINUE T1 = ZERO T2 = ZERO DO 100 J = 1,N T1 = T1 + DFLOAT(J)* (X(J)-ONE) T2 = T2 + (X(J)-ONE)**2 100 CONTINUE F = T2 + T1**2* (ONE+T1**2) GO TO 390 C C WATSON FUNCTION. C 110 CONTINUE F = ZERO DO 140 I = 1,29 D1 = DFLOAT(I)/C29 S1 = ZERO D2 = ONE DO 120 J = 2,N S1 = S1 + DFLOAT(J-1)*D2*X(J) D2 = D1*D2 120 CONTINUE S2 = ZERO D2 = ONE DO 130 J = 1,N S2 = S2 + D2*X(J) D2 = D1*D2 130 CONTINUE T = S1 - S2**2 - ONE F = F + T**2 140 CONTINUE T1 = X(2) - X(1)**2 - ONE F = F + X(1)**2 + T1**2 GO TO 390 C C PENALTY FUNCTION I. C 150 CONTINUE T1 = -CP25 T2 = ZERO DO 160 J = 1,N T1 = T1 + X(J)**2 T2 = T2 + (X(J)-ONE)**2 160 CONTINUE F = AP*T2 + BP*T1**2 GO TO 390 C C PENALTY FUNCTION II. C 170 CONTINUE T1 = -ONE T2 = ZERO T3 = ZERO D1 = DEXP(CP1) D2 = ONE DO 190 J = 1,N T1 = T1 + DFLOAT(N-J+1)*X(J)**2 S1 = DEXP(X(J)/TEN) IF (J.EQ.1) GO TO 180 S3 = S1 + S2 - D2* (D1+ONE) T2 = T2 + S3**2 T3 = T3 + (S1-ONE/D1)**2 180 CONTINUE S2 = S1 D2 = D1*D2 190 CONTINUE F = AP* (T2+T3) + BP* (T1**2+ (X(1)-CP2)**2) GO TO 390 C C BROWN BADLY SCALED FUNCTION. C 200 CONTINUE T1 = X(1) - C1PD6 T2 = X(2) - C2PDM6 T3 = X(1)*X(2) - TWO F = T1**2 + T2**2 + T3**2 GO TO 390 C C BROWN AND DENNIS FUNCTION. C 210 CONTINUE F = ZERO DO 220 I = 1,20 D1 = DFLOAT(I)/FIVE D2 = DSIN(D1) T1 = X(1) + D1*X(2) - DEXP(D1) T2 = X(3) + D2*X(4) - DCOS(D1) T = T1**2 + T2**2 F = F + T**2 220 CONTINUE GO TO 390 C C GULF RESEARCH AND DEVELOPMENT FUNCTION. C 230 CONTINUE F = ZERO D1 = TWO/THREE DO 240 I = 1,99 ARG = DFLOAT(I)/C100 R = DABS((-FIFTY*DLOG(ARG))**D1+C25-X(2)) T1 = R**X(3)/X(1) T2 = DEXP(-T1) T = T2 - ARG F = F + T**2 240 CONTINUE GO TO 390 C C TRIGONOMETRIC FUNCTION. C 250 CONTINUE S1 = ZERO DO 260 J = 1,N S1 = S1 + DCOS(X(J)) 260 CONTINUE F = ZERO DO 270 J = 1,N T = DFLOAT(N+J) - DSIN(X(J)) - S1 - DFLOAT(J)*DCOS(X(J)) F = F + T**2 270 CONTINUE GO TO 390 C C EXTENDED ROSENBROCK FUNCTION. C 280 CONTINUE F = ZERO DO 290 J = 1,N,2 T1 = ONE - X(J) T2 = TEN* (X(J+1)-X(J)**2) F = F + T1**2 + T2**2 290 CONTINUE GO TO 390 C C EXTENDED POWELL FUNCTION. C 300 CONTINUE F = ZERO DO 310 J = 1,N,4 T = X(J) + TEN*X(J+1) T1 = X(J+2) - X(J+3) S1 = FIVE*T1 T2 = X(J+1) - TWO*X(J+2) S2 = T2**3 T3 = X(J) - X(J+3) S3 = TEN*T3**3 F = F + T**2 + S1*T1 + S2*T2 + S3*T3 310 CONTINUE GO TO 390 C C BEALE FUNCTION. C 320 CONTINUE S1 = ONE - X(2) T1 = C1P5 - X(1)*S1 S2 = ONE - X(2)**2 T2 = C2P25 - X(1)*S2 S3 = ONE - X(2)**3 T3 = C2P625 - X(1)*S3 F = T1**2 + T2**2 + T3**2 GO TO 390 C C WOOD FUNCTION. C 330 CONTINUE S1 = X(2) - X(1)**2 S2 = ONE - X(1) S3 = X(2) - ONE T1 = X(4) - X(3)**2 T2 = ONE - X(3) T3 = X(4) - ONE F = C100*S1**2 + S2**2 + C90*T1**2 + T2**2 + TEN* (S3+T3)**2 + + (S3-T3)**2/TEN GO TO 390 C C CHEBYQUAD FUNCTION. C 340 CONTINUE DO 350 I = 1,N FVEC(I) = ZERO 350 CONTINUE DO 370 J = 1,N T1 = ONE T2 = TWO*X(J) - ONE T = TWO*T2 DO 360 I = 1,N FVEC(I) = FVEC(I) + T2 TH = T*T2 - T1 T1 = T2 T2 = TH 360 CONTINUE 370 CONTINUE F = ZERO D1 = ONE/DFLOAT(N) IEV = -1 DO 380 I = 1,N T = D1*FVEC(I) IF (IEV.GT.0) T = T + ONE/ (DFLOAT(I)**2-ONE) F = F + T**2 IEV = -IEV 380 CONTINUE 390 CONTINUE RETURN C C LAST CARD OF SUBROUTINE OBJFCN. C END SUBROUTINE GRDFCN(N,X,G,NPROB) C ********** C C SUBROUTINE GRDFCN C C THIS SUBROUTINE DEFINES THE GRADIENT VECTORS OF EIGHTEEN C NONLINEAR UNCONSTRAINED MINIMIZATION PROBLEMS. THE PROBLEM C DIMENSIONS ARE AS DESCRIBED IN THE PROLOGUE COMMENTS OF OBJFCN. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE GRDFCN(N,X,G,NPROB) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE. C C X IS AN INPUT ARRAY OF LENGTH N. C C G IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE COMPONENTS C OF THE GRADIENT VECTOR OF THE NPROB OBJECTIVE FUNCTION C EVALUATED AT X. C C NPROB IS A POSITIVE INTEGER INPUT VARIABLE WHICH DEFINES THE C NUMBER OF THE PROBLEM. NPROB MUST NOT EXCEED 18. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... DABS,DATAN,DCOS,DEXP,DLOG,DSIGN,DSIN, C DSQRT C C MINPACK. VERSION OF JULY 1978. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** C .. Scalar Arguments .. INTEGER N,NPROB C .. C .. Array Arguments .. DOUBLE PRECISION G(N),X(N) C .. C .. Local Scalars .. DOUBLE PRECISION AP,ARG,BP,C100,C10000,C180,C19P8,C1P5,C1PD6,C200, + C20P2,C25,C29,C2P25,C2P625,C2PDM6,C3P5,CP0001, + CP1,CP2,CP25,CP5,D1,D2,EIGHT,FIFTY,FIVE,FOUR,ONE, + R,S1,S2,S3,T,T1,T2,T3,TEN,TH,THREE,TPI,TWENTY, + TWO,ZERO INTEGER I,IEV,IVAR,J C .. C .. Local Arrays .. DOUBLE PRECISION FVEC(50),Y(15) C .. C .. Intrinsic Functions .. INTRINSIC DABS,DATAN,DCOS,DEXP,DLOG,DSIGN,DSIN,DSQRT C .. C .. Statement Functions .. DOUBLE PRECISION DFLOAT C .. C .. Data statements .. DATA ZERO,ONE,TWO,THREE,FOUR,FIVE,EIGHT,TEN,TWENTY,FIFTY/0.0D0, + 1.0D0,2.0D0,3.0D0,4.0D0,5.0D0,8.0D0,1.0D1,2.0D1,5.0D1/ DATA C2PDM6,CP0001,CP1,CP2,CP25,CP5,C1P5,C2P25,C2P625,C3P5,C19P8, + C20P2,C25,C29,C100,C180,C200,C10000,C1PD6/2.0D-6,1.0D-4, + 1.0D-1,2.0D-1,2.5D-1,5.0D-1,1.5D0,2.25D0,2.625D0,3.5D0, + 1.98D1,2.02D1,2.5D1,2.9D1,1.0D2,1.8D2,2.0D2,1.0D4,1.0D6/ DATA AP,BP/1.0D-5,1.0D0/ DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8),Y(9),Y(10),Y(11), + Y(12),Y(13),Y(14),Y(15)/9.0D-4,4.4D-3,1.75D-2,5.4D-2, + 1.295D-1,2.42D-1,3.521D-1,3.989D-1,3.521D-1,2.42D-1,1.295D-1, + 5.4D-2,1.75D-2,4.4D-3,9.0D-4/ C .. C .. Statement Function definitions .. DFLOAT(IVAR) = IVAR C .. C C GRADIENT ROUTINE SELECTOR. C GO TO (10,20,50,70,80,100,130,190,220, + 260,270,290,310,350,370,390,400, + 410) NPROB C C HELICAL VALLEY FUNCTION. C 10 CONTINUE TPI = EIGHT*DATAN(ONE) TH = DSIGN(CP25,X(2)) IF (X(1).GT.ZERO) TH = DATAN(X(2)/X(1))/TPI IF (X(1).LT.ZERO) TH = DATAN(X(2)/X(1))/TPI + CP5 ARG = X(1)**2 + X(2)**2 R = DSQRT(ARG) T = X(3) - TEN*TH S1 = TEN*T/ (TPI*ARG) G(1) = C200* (X(1)-X(1)/R+X(2)*S1) G(2) = C200* (X(2)-X(2)/R-X(1)*S1) G(3) = TWO* (C100*T+X(3)) GO TO 490 C C BIGGS EXP6 FUNCTION. C 20 CONTINUE DO 30 J = 1,N G(J) = ZERO 30 CONTINUE DO 40 I = 1,13 c DO 40 I = 1, 6 D1 = DFLOAT(I)/TEN D2 = DEXP(-D1) - FIVE*DEXP(-TEN*D1) + THREE*DEXP(-FOUR*D1) S1 = DEXP(-D1*X(1)) S2 = DEXP(-D1*X(2)) S3 = DEXP(-D1*X(5)) T = X(3)*S1 - X(4)*S2 + X(6)*S3 - D2 TH = D1*T G(1) = G(1) - S1*TH G(2) = G(2) + S2*TH G(3) = G(3) + S1*T G(4) = G(4) - S2*T G(5) = G(5) - S3*TH G(6) = G(6) + S3*T 40 CONTINUE G(1) = TWO*X(3)*G(1) G(2) = TWO*X(4)*G(2) G(3) = TWO*G(3) G(4) = TWO*G(4) G(5) = TWO*X(6)*G(5) G(6) = TWO*G(6) GO TO 490 C C GAUSSIAN FUNCTION. C 50 CONTINUE G(1) = ZERO G(2) = ZERO G(3) = ZERO DO 60 I = 1,15 D1 = CP5*DFLOAT(I-1) D2 = C3P5 - D1 - X(3) ARG = -CP5*X(2)*D2**2 R = DEXP(ARG) T = X(1)*R - Y(I) S1 = R*T S2 = D2*S1 G(1) = G(1) + S1 G(2) = G(2) - D2*S2 G(3) = G(3) + S2 60 CONTINUE G(1) = TWO*G(1) G(2) = X(1)*G(2) G(3) = TWO*X(1)*X(2)*G(3) GO TO 490 C C POWELL BADLY SCALED FUNCTION. C 70 CONTINUE T1 = C10000*X(1)*X(2) - ONE S1 = DEXP(-X(1)) S2 = DEXP(-X(2)) T2 = S1 + S2 - ONE - CP0001 G(1) = TWO* (C10000*X(2)*T1-S1*T2) G(2) = TWO* (C10000*X(1)*T1-S2*T2) GO TO 490 C C BOX 3-DIMENSIONAL FUNCTION. C 80 CONTINUE G(1) = ZERO G(2) = ZERO G(3) = ZERO DO 90 I = 1,10 D1 = DFLOAT(I) D2 = D1/TEN S1 = DEXP(-D2*X(1)) S2 = DEXP(-D2*X(2)) S3 = DEXP(-D2) - DEXP(-D1) T = S1 - S2 - S3*X(3) TH = D2*T G(1) = G(1) - S1*TH G(2) = G(2) + S2*TH G(3) = G(3) - S3*T 90 CONTINUE G(1) = TWO*G(1) G(2) = TWO*G(2) G(3) = TWO*G(3) GO TO 490 C C VARIABLY DIMENSIONED FUNCTION. C 100 CONTINUE T1 = ZERO DO 110 J = 1,N T1 = T1 + DFLOAT(J)* (X(J)-ONE) 110 CONTINUE T = T1* (ONE+TWO*T1**2) DO 120 J = 1,N G(J) = TWO* (X(J)-ONE+DFLOAT(J)*T) 120 CONTINUE GO TO 490 C C WATSON FUNCTION. C 130 CONTINUE DO 140 J = 1,N G(J) = ZERO 140 CONTINUE DO 180 I = 1,29 D1 = DFLOAT(I)/C29 S1 = ZERO D2 = ONE DO 150 J = 2,N S1 = S1 + DFLOAT(J-1)*D2*X(J) D2 = D1*D2 150 CONTINUE S2 = ZERO D2 = ONE DO 160 J = 1,N S2 = S2 + D2*X(J) D2 = D1*D2 160 CONTINUE T = S1 - S2**2 - ONE S3 = TWO*D1*S2 D2 = TWO/D1 DO 170 J = 1,N G(J) = G(J) + D2* (DFLOAT(J-1)-S3)*T D2 = D1*D2 170 CONTINUE 180 CONTINUE T1 = X(2) - X(1)**2 - ONE G(1) = G(1) + X(1)* (TWO-FOUR*T1) G(2) = G(2) + TWO*T1 GO TO 490 C C PENALTY FUNCTION I. C 190 CONTINUE T1 = -CP25 DO 200 J = 1,N T1 = T1 + X(J)**2 200 CONTINUE D1 = TWO*AP TH = FOUR*BP*T1 DO 210 J = 1,N G(J) = D1* (X(J)-ONE) + X(J)*TH 210 CONTINUE GO TO 490 C C PENALTY FUNCTION II. C 220 CONTINUE T1 = -ONE DO 230 J = 1,N T1 = T1 + DFLOAT(N-J+1)*X(J)**2 230 CONTINUE D1 = DEXP(CP1) D2 = ONE TH = FOUR*BP*T1 DO 250 J = 1,N G(J) = DFLOAT(N-J+1)*X(J)*TH S1 = DEXP(X(J)/TEN) IF (J.EQ.1) GO TO 240 S3 = S1 + S2 - D2* (D1+ONE) G(J) = G(J) + AP*S1* (S3+S1-ONE/D1)/FIVE G(J-1) = G(J-1) + AP*S2*S3/FIVE 240 CONTINUE S2 = S1 D2 = D1*D2 250 CONTINUE G(1) = G(1) + TWO*BP* (X(1)-CP2) GO TO 490 C C BROWN BADLY SCALED FUNCTION. C 260 CONTINUE T1 = X(1) - C1PD6 T2 = X(2) - C2PDM6 T3 = X(1)*X(2) - TWO G(1) = TWO* (T1+X(2)*T3) G(2) = TWO* (T2+X(1)*T3) GO TO 490 C C BROWN AND DENNIS FUNCTION. C 270 CONTINUE G(1) = ZERO G(2) = ZERO G(3) = ZERO G(4) = ZERO DO 280 I = 1,20 D1 = DFLOAT(I)/FIVE D2 = DSIN(D1) T1 = X(1) + D1*X(2) - DEXP(D1) T2 = X(3) + D2*X(4) - DCOS(D1) T = T1**2 + T2**2 S1 = T1*T S2 = T2*T G(1) = G(1) + S1 G(2) = G(2) + D1*S1 G(3) = G(3) + S2 G(4) = G(4) + D2*S2 280 CONTINUE G(1) = FOUR*G(1) G(2) = FOUR*G(2) G(3) = FOUR*G(3) G(4) = FOUR*G(4) GO TO 490 C C GULF RESEARCH AND DEVELOPMENT FUNCTION. C 290 CONTINUE G(1) = ZERO G(2) = ZERO G(3) = ZERO D1 = TWO/THREE DO 300 I = 1,99 ARG = DFLOAT(I)/C100 R = DABS((-FIFTY*DLOG(ARG))**D1+C25-X(2)) T1 = R**X(3)/X(1) T2 = DEXP(-T1) T = T2 - ARG S1 = T1*T2*T G(1) = G(1) + S1 G(2) = G(2) + S1/R G(3) = G(3) - S1*DLOG(R) 300 CONTINUE G(1) = TWO*G(1)/X(1) G(2) = TWO*X(3)*G(2) G(3) = TWO*G(3) GO TO 490 C C TRIGONOMETRIC FUNCTION. C 310 CONTINUE S1 = ZERO DO 320 J = 1,N G(J) = DCOS(X(J)) S1 = S1 + G(J) 320 CONTINUE S2 = ZERO DO 330 J = 1,N TH = DSIN(X(J)) T = DFLOAT(N+J) - TH - S1 - DFLOAT(J)*G(J) S2 = S2 + T G(J) = (DFLOAT(J)*TH-G(J))*T 330 CONTINUE DO 340 J = 1,N G(J) = TWO* (G(J)+DSIN(X(J))*S2) 340 CONTINUE GO TO 490 C C EXTENDED ROSENBROCK FUNCTION. C 350 CONTINUE DO 360 J = 1,N,2 T1 = ONE - X(J) G(J+1) = C200* (X(J+1)-X(J)**2) G(J) = -TWO* (X(J)*G(J+1)+T1) 360 CONTINUE GO TO 490 C C EXTENDED POWELL FUNCTION. C 370 CONTINUE DO 380 J = 1,N,4 T = X(J) + TEN*X(J+1) T1 = X(J+2) - X(J+3) S1 = FIVE*T1 T2 = X(J+1) - TWO*X(J+2) S2 = FOUR*T2**3 T3 = X(J) - X(J+3) S3 = TWENTY*T3**3 G(J) = TWO* (T+S3) G(J+1) = TWENTY*T + S2 G(J+2) = TWO* (S1-S2) G(J+3) = -TWO* (S1+S3) 380 CONTINUE GO TO 490 C C BEALE FUNCTION. C 390 CONTINUE S1 = ONE - X(2) T1 = C1P5 - X(1)*S1 S2 = ONE - X(2)**2 T2 = C2P25 - X(1)*S2 S3 = ONE - X(2)**3 T3 = C2P625 - X(1)*S3 G(1) = -TWO* (S1*T1+S2*T2+S3*T3) G(2) = TWO*X(1)* (T1+X(2)* (TWO*T2+THREE*X(2)*T3)) GO TO 490 C C WOOD FUNCTION. C 400 CONTINUE S1 = X(2) - X(1)**2 S2 = ONE - X(1) S3 = X(2) - ONE T1 = X(4) - X(3)**2 T2 = ONE - X(3) T3 = X(4) - ONE G(1) = -TWO* (C200*X(1)*S1+S2) G(2) = C200*S1 + C20P2*S3 + C19P8*T3 G(3) = -TWO* (C180*X(3)*T1+T2) G(4) = C180*T1 + C20P2*T3 + C19P8*S3 GO TO 490 C C CHEBYQUAD FUNCTION. C 410 CONTINUE DO 420 I = 1,N FVEC(I) = ZERO 420 CONTINUE DO 440 J = 1,N T1 = ONE T2 = TWO*X(J) - ONE T = TWO*T2 DO 430 I = 1,N FVEC(I) = FVEC(I) + T2 TH = T*T2 - T1 T1 = T2 T2 = TH 430 CONTINUE 440 CONTINUE D1 = ONE/DFLOAT(N) IEV = -1 DO 450 I = 1,N FVEC(I) = D1*FVEC(I) IF (IEV.GT.0) FVEC(I) = FVEC(I) + ONE/ (DFLOAT(I)**2-ONE) IEV = -IEV 450 CONTINUE DO 470 J = 1,N G(J) = ZERO T1 = ONE T2 = TWO*X(J) - ONE T = TWO*T2 S1 = ZERO S2 = TWO DO 460 I = 1,N G(J) = G(J) + FVEC(I)*S2 TH = FOUR*T2 + T*S2 - S1 S1 = S2 S2 = TH TH = T*T2 - T1 T1 = T2 T2 = TH 460 CONTINUE 470 CONTINUE D2 = TWO*D1 DO 480 J = 1,N G(J) = D2*G(J) 480 CONTINUE 490 CONTINUE RETURN C C LAST CARD OF SUBROUTINE GRDFCN. C END SUBROUTINE HESFCN(N,X,HESD,HESL,NPROB) C ********** C C SUBROUTINE HESFCN C C THIS SUBROUTINE DEFINES THE HESSIAN MATRICES OF 18 C NONLINEAR UNCONSTRAINED MINIMIZATION PROBLEMS. THE PROBLEM C DIMENSIONS ARE AS DESCRIBED IN OBJFCN. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE HESFCN (N, X, HESD, HESL, NPROB) C C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE. C C X IS AN INPUT ARRAY OF LENGTH N. C C HESD IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE C DIAGONAL COMPONENTS OF THE HESSIAN MATRIX OF THE NPROB C OBJECTIVE FUNCTION EVALUATED AT X. C C HESL IS AN OUTPUT ARRAY OF LENGTH N*(N-1)/2 WHICH CONTAINS C THE LOWER TRIANGULAR PART OF THE HESSIAN MATRIX OF THE C NPROB OBJECTIVE FUNCTION EVALUATED AT X. C C NPROB IS A POSITIVE INTEGER INPUT VARIABLE WHICH DEFINES THE C NUMBER OF THE PROBLEM. NPROB MUST NOT EXCEED 18. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... ABS, ATAN, COS, EXP, LOG, SIGN, SIN, C SQRT C C INTEGER INLINE FUNCTION IX GIVES THE LOCATION OF A HESSIAN C ELEMENT (I,J), I>J, IN HESL C C VICTORIA Z. AVERBUKH, SAMUEL A. FIGUEROA, AND C TAMAR SCHLICK, 1993. C ********** C .. Parameters .. DOUBLE PRECISION ZERO,ONE,TWO,THREE,FOUR,FIVE,SIX,EIGHT,NINE,TEN, + FIFTY,CP0001,CP1,CP25,CP5,C1P5,C2P25,C2P625,C3P5, + C12,C19P8,C25,C29,C50,C100,C120,C200,C200P2,C202, + C220P2,C360,C400,C1000,C1080,C1200,C2000,C20000, + C2E8,C4E8,AP,BP,PI PARAMETER (ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,FOUR=4.0D0, + FIVE=5.0D0,SIX=6.0D0,EIGHT=8.0D0,NINE=9.0D0,TEN=1.0D1, + FIFTY=5.0D1,CP0001=1.0D-4,CP1=1.0D-1,CP25=2.5D-1, + CP5=5.0D-1,C1P5=1.5D0,C2P25=2.25D0,C2P625=2.625D0, + C3P5=3.5D0,C12=1.2D1,C19P8=1.98D1,C25=2.5D1,C29=2.9D1, + C50=5.0D1,C100=1.0D2,C120=1.2D2,C200=2.0D2, + C200P2=2.002D2,C202=2.02D2,C220P2=2.202D2,C360=3.6D2, + C400=4.0D2,C1000=1.0D3,C1080=1.08D3,C1200=1.2D3, + C2000=2.0D3,C20000=2.0D4,C2E8=2.0D8,C4E8=4.0D8, + AP=1.0D-5,BP=ONE,PI=3.141592653589793D0) C .. C .. Scalar Arguments .. INTEGER N,NPROB C .. C .. Array Arguments .. DOUBLE PRECISION HESD(N),HESL(*),X(N) C .. C .. Local Scalars .. DOUBLE PRECISION ARG,D1,D2,D3,LOGR,P1,P2,PIARG,PIARG2,R,R3INV,S1, + S1S2,S1S3,S2,S2S3,S3,SS1,SS2,T,T1,T2,T3,TH,TT, + TT1,TT2,TTH INTEGER I,II,IVAR,J,JJ,K,M LOGICAL IEV C .. C .. Local Arrays .. DOUBLE PRECISION FVEC(50),GVEC(50),Y(15) C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,COS,EXP,FLOAT,LOG,SIGN,SIN,SQRT C .. C .. Statement Functions .. DOUBLE PRECISION DFLOAT INTEGER IX C .. C .. Statement Function definitions .. IX(II,JJ) = (II-1)* (II-2)/2 + JJ DFLOAT(IVAR) = IVAR C .. C .. Data statements .. DATA Y/9.0D-4,4.4D-3,1.75D-2,5.4D-2,1.295D-1,2.42D-1,3.521D-1, + 3.989D-1,3.521D-1,2.42D-1,1.295D-1,5.4D-2,1.75D-2,4.4D-3, + 9.0D-4/ C .. C C HESSIAN ROUTINE SELECTOR. C GO TO (10,20,80,100,110,130,170,260,300, + 340,350,390,430,480,510,540,550, + 560) NPROB C C HELICAL VALLEY FUNCTION. C 10 CONTINUE C IF (X(1).EQ.ZERO) THEN TH = SIGN(CP25,X(2)) ELSE TH = ATAN(X(2)/X(1))/ (TWO*PI) IF (X(1).LT.ZERO) TH = TH + CP5 END IF ARG = X(1)**2 + X(2)**2 PIARG = PI*ARG PIARG2 = PIARG*ARG R3INV = ONE/SQRT(ARG)**3 T = X(3) - TEN*TH S1 = FIVE*T/PIARG P1 = C2000*X(1)*X(2)*T/PIARG2 P2 = (FIVE/PIARG)**2 HESD(1) = C200 - C200* (R3INV-P2)*X(2)**2 - P1 HESD(2) = C200 - C200* (R3INV-P2)*X(1)**2 + P1 HESD(3) = C202 HESL(1) = C200*X(1)*X(2)*R3INV + + C1000/PIARG2* (T* (X(1)**2-X(2)**2)-FIVE*X(1)*X(2)/PI) HESL(2) = C1000*X(2)/PIARG HESL(3) = -C1000*X(1)/PIARG RETURN C C BIGGS EXP6 FUNCTION. C 20 CONTINUE DO 30 I = 1,6 HESD(I) = ZERO 30 CONTINUE DO 40 I = 1,15 HESL(I) = ZERO 40 CONTINUE c DO 230 I = 1, 6 DO 50 I = 1,13 D1 = DFLOAT(I)/TEN D2 = EXP(-D1) - FIVE*EXP(-TEN*D1) + THREE*EXP(-FOUR*D1) S1 = EXP(-D1*X(1)) S2 = EXP(-D1*X(2)) S3 = EXP(-D1*X(5)) T = X(3)*S1 - X(4)*S2 + X(6)*S3 - D2 D2 = D1**2 S1S2 = S1*S2 S1S3 = S1*S3 S2S3 = S2*S3 HESD(1) = HESD(1) + D2*S1* (T+X(3)*S1) HESD(2) = HESD(2) - D2*S2* (T-X(4)*S2) HESD(3) = HESD(3) + S1**2 HESD(4) = HESD(4) + S2**2 HESD(5) = HESD(5) + D2*S3* (T+X(6)*S3) HESD(6) = HESD(6) + S3**2 HESL(1) = HESL(1) - D2*S1S2 HESL(2) = HESL(2) - D1*S1* (T+X(3)*S1) HESL(3) = HESL(3) + D1*S1S2 HESL(4) = HESL(4) + D1*S1S2 HESL(5) = HESL(5) + D1*S2* (T-X(4)*S2) HESL(6) = HESL(6) - S1S2 HESL(7) = HESL(7) + D2*S1S3 HESL(8) = HESL(8) - D2*S2S3 HESL(9) = HESL(9) - D1*S1S3 HESL(10) = HESL(10) + D1*S2S3 HESL(11) = HESL(11) - D1*S1S3 HESL(12) = HESL(12) + D1*S2S3 HESL(13) = HESL(13) + S1S3 HESL(14) = HESL(14) - S2S3 HESL(15) = HESL(15) - D1*S3* (T+X(6)*S3) 50 CONTINUE HESD(1) = X(3)*HESD(1) HESD(2) = X(4)*HESD(2) HESD(5) = X(6)*HESD(5) HESL(1) = X(3)*X(4)*HESL(1) HESL(3) = X(4)*HESL(3) HESL(4) = X(3)*HESL(4) HESL(7) = X(3)*X(6)*HESL(7) HESL(8) = X(4)*X(6)*HESL(8) HESL(9) = X(6)*HESL(9) HESL(10) = X(6)*HESL(10) HESL(11) = X(3)*HESL(11) HESL(12) = X(4)*HESL(12) DO 60 I = 1,6 HESD(I) = TWO*HESD(I) 60 CONTINUE DO 70 I = 1,15 HESL(I) = TWO*HESL(I) 70 CONTINUE RETURN C C GAUSSIAN FUNCTION. C 80 CONTINUE HESD(1) = ZERO HESD(2) = ZERO HESD(3) = ZERO HESL(1) = ZERO HESL(2) = ZERO HESL(3) = ZERO DO 90 I = 1,15 D1 = CP5*DFLOAT(I-1) D2 = C3P5 - D1 - X(3) ARG = CP5*X(2)*D2**2 R = EXP(-ARG) T = X(1)*R - Y(I) T1 = TWO*X(1)*R - Y(I) HESD(1) = HESD(1) + R**2 HESD(2) = HESD(2) + R*T1*D2**4 HESD(3) = HESD(3) + R* (X(2)*T1*D2**2-T) HESL(1) = HESL(1) - R*T1*D2**2 HESL(2) = HESL(2) + D2*R*T1 HESL(3) = HESL(3) + D2*R* (T-ARG*T1) 90 CONTINUE HESD(1) = TWO*HESD(1) HESD(2) = CP5*X(1)*HESD(2) HESD(3) = TWO*X(1)*X(2)*HESD(3) HESL(2) = TWO*X(2)*HESL(2) HESL(3) = TWO*X(1)*HESL(3) RETURN C C POWELL BADLY SCALED FUNCTION. C 100 CONTINUE S1 = EXP(-X(1)) S2 = EXP(-X(2)) T2 = S1 + S2 - ONE - CP0001 HESD(1) = C2E8*X(2)**2 + TWO*S1* (S1+T2) HESD(2) = C2E8*X(1)**2 + TWO*S2* (S2+T2) HESL(1) = C4E8*X(1)*X(2) + TWO*S1*S2 - C20000 RETURN C C BOX 3-DIMENSIONAL FUNCTION. C 110 CONTINUE HESD(1) = ZERO HESD(2) = ZERO HESD(3) = ZERO HESL(1) = ZERO HESL(2) = ZERO HESL(3) = ZERO DO 120 I = 1,10 D1 = DFLOAT(I) D2 = D1/TEN S1 = EXP(-D2*X(1)) S2 = EXP(-D2*X(2)) S3 = EXP(-D2) - EXP(-D1) T = S1 - S2 - S3*X(3) TH = T*D2**2 HESD(1) = HESD(1) + TH*S1 + (D2*S1)**2 HESD(2) = HESD(2) - TH*S2 + (D2*S2)**2 HESD(3) = HESD(3) + S3**2 HESL(1) = HESL(1) - S1*S2*D2**2 HESL(2) = HESL(2) + D2*S1*S3 HESL(3) = HESL(3) - D2*S2*S3 120 CONTINUE HESD(1) = TWO*HESD(1) HESD(2) = TWO*HESD(2) HESD(3) = TWO*HESD(3) HESL(1) = TWO*HESL(1) HESL(2) = TWO*HESL(2) HESL(3) = TWO*HESL(3) RETURN C C VARIABLY DIMENSIONED FUNCTION. C 130 CONTINUE T1 = ZERO DO 140 J = 1,N T1 = T1 + DFLOAT(J)* (X(J)-ONE) 140 CONTINUE T = ONE + SIX*T1**2 M = 0 DO 160 J = 1,N HESD(J) = TWO + TWO*T*DFLOAT(J)**2 DO 150 K = 1,J - 1 M = M + 1 HESL(M) = TWO*T*DFLOAT(J*K) 150 CONTINUE 160 CONTINUE RETURN C C WATSON FUNCTION. C 170 CONTINUE DO 180 J = 1,N HESD(J) = ZERO 180 CONTINUE DO 190 J = 1,N* (N-1)/2 HESL(J) = ZERO 190 CONTINUE DO 230 I = 1,29 D1 = DFLOAT(I)/C29 D2 = ONE S1 = ZERO S2 = X(1) DO 200 J = 2,N S1 = S1 + DFLOAT(J-1)*D2*X(J) D2 = D1*D2 S2 = S2 + D2*X(J) 200 CONTINUE T = TWO* (S1-S2**2-ONE)*D1**2 S3 = TWO*D1*S2 D2 = ONE/D1 M = 0 DO 220 J = 1,N T1 = DFLOAT(J-1) - S3 HESD(J) = HESD(J) + (T1**2-T)*D2**2 D3 = ONE/D1 DO 210 K = 1,J - 1 M = M + 1 HESL(M) = HESL(M) + (T1* (DFLOAT(K-1)-S3)-T)*D2*D3 D3 = D1*D3 210 CONTINUE D2 = D1*D2 220 CONTINUE 230 CONTINUE T3 = X(2) - X(1)**2 - ONE HESD(1) = HESD(1) + ONE - TWO* (T3-TWO*X(1)**2) HESD(2) = HESD(2) + ONE HESL(1) = HESL(1) - TWO*X(1) DO 240 J = 1,N HESD(J) = TWO*HESD(J) 240 CONTINUE DO 250 J = 1,N* (N-1)/2 HESL(J) = TWO*HESL(J) 250 CONTINUE RETURN C C PENALTY FUNCTION I. C 260 CONTINUE T1 = -CP25 DO 270 J = 1,N T1 = T1 + X(J)**2 270 CONTINUE D1 = TWO*AP TH = FOUR*BP*T1 M = 0 DO 290 J = 1,N HESD(J) = D1 + TH + EIGHT*X(J)**2 DO 280 K = 1,J - 1 M = M + 1 HESL(M) = EIGHT*X(J)*X(K) 280 CONTINUE 290 CONTINUE RETURN C C PENALTY FUNCTION II. C 300 CONTINUE T1 = -ONE DO 310 J = 1,N T1 = T1 + DFLOAT(N-J+1)*X(J)**2 310 CONTINUE D1 = EXP(CP1) D2 = ONE TH = FOUR*BP*T1 M = 0 DO 330 J = 1,N HESD(J) = EIGHT*BP* (DFLOAT(N-J+1)*X(J))**2 + DFLOAT(N-J+1)*TH S1 = EXP(X(J)/TEN) IF (J.GT.1) THEN S3 = S1 + S2 - D2* (D1+ONE) HESD(J) = HESD(J) + AP*S1* (S3+S1-ONE/D1+TWO*S1)/C50 HESD(J-1) = HESD(J-1) + AP*S2* (S2+S3)/C50 DO 320 K = 1,J - 1 M = M + 1 T1 = EXP(DFLOAT(K)/TEN) HESL(M) = EIGHT*DFLOAT(N-J+1)*DFLOAT(N-K+1)*X(J)*X(K) 320 CONTINUE HESL(M) = HESL(M) + AP*S1*S2/C50 END IF S2 = S1 D2 = D1*D2 330 CONTINUE HESD(1) = HESD(1) + TWO*BP RETURN C C BROWN BADLY SCALED FUNCTION. C 340 CONTINUE HESD(1) = TWO + TWO*X(2)**2 HESD(2) = TWO + TWO*X(1)**2 HESL(1) = FOUR*X(1)*X(2) - FOUR RETURN C C BROWN AND DENNIS FUNCTION. C 350 CONTINUE DO 360 I = 1,4 HESD(I) = ZERO 360 CONTINUE DO 370 I = 1,6 HESL(I) = ZERO 370 CONTINUE DO 380 I = 1,20 D1 = DFLOAT(I)/FIVE D2 = SIN(D1) T1 = X(1) + D1*X(2) - EXP(D1) T2 = X(3) + D2*X(4) - COS(D1) T = EIGHT*T1*T2 S1 = C12*T1**2 + FOUR*T2**2 S2 = C12*T2**2 + FOUR*T1**2 HESD(1) = HESD(1) + S1 HESD(2) = HESD(2) + S1*D1**2 HESD(3) = HESD(3) + S2 HESD(4) = HESD(4) + S2*D2**2 HESL(1) = HESL(1) + S1*D1 HESL(2) = HESL(2) + T HESL(4) = HESL(4) + T*D2 HESL(3) = HESL(3) + T*D1 HESL(5) = HESL(5) + T*D1*D2 HESL(6) = HESL(6) + S2*D2 380 CONTINUE RETURN C C GULF RESEARCH AND DEVELOPMENT FUNCTION. C 390 CONTINUE DO 400 I = 1,3 HESD(I) = ZERO HESL(I) = ZERO 400 CONTINUE D1 = TWO/THREE DO 410 I = 1,99 ARG = DFLOAT(I)/C100 R = (-FIFTY*LOG(ARG))**D1 + C25 - X(2) T1 = ABS(R)**X(3)/X(1) T2 = EXP(-T1) T3 = T1*T2* (T1*T2+ (T1-ONE)* (T2-ARG)) T = T1*T2* (T2-ARG) LOGR = LOG(ABS(R)) HESD(1) = HESD(1) + T3 - T HESD(2) = HESD(2) + (T+X(3)*T3)/R**2 HESD(3) = HESD(3) + T3*LOGR**2 HESL(1) = HESL(1) + T3/R HESL(2) = HESL(2) - T3*LOGR HESL(3) = HESL(3) + (T-X(3)*T3*LOGR)/R 410 CONTINUE HESD(1) = HESD(1)/X(1)**2 HESD(2) = HESD(2)*X(3) HESL(1) = HESL(1)*X(3)/X(1) HESL(2) = HESL(2)/X(1) DO 420 I = 1,3 HESD(I) = TWO*HESD(I) HESL(I) = TWO*HESL(I) 420 CONTINUE RETURN C C TRIGONOMETRIC FUNCTION. C 430 CONTINUE S1 = ZERO DO 440 J = 1,N HESD(J) = SIN(X(J)) S1 = S1 + COS(X(J)) 440 CONTINUE S2 = ZERO M = 0 DO 460 J = 1,N TH = COS(X(J)) T = DFLOAT(N+J) - HESD(J) - S1 - DFLOAT(J)*TH S2 = S2 + T DO 450 K = 1,J - 1 M = M + 1 HESL(M) = SIN(X(K))* (DFLOAT(N+J+K)*HESD(J)-TH) - + HESD(J)*COS(X(K)) HESL(M) = TWO*HESL(M) 450 CONTINUE HESD(J) = DFLOAT(J* (J+2)+N)*HESD(J)**2 + + TH* (TH-DFLOAT(2*J+2)*HESD(J)) + + T* (DFLOAT(J)*TH+HESD(J)) 460 CONTINUE DO 470 J = 1,N HESD(J) = TWO* (HESD(J)+COS(X(J))*S2) 470 CONTINUE RETURN C C EXTENDED ROSENBROCK FUNCTION. C 480 CONTINUE DO 490 J = 1,N* (N-1)/2 HESL(J) = ZERO 490 CONTINUE DO 500 J = 1,N,2 HESD(J+1) = C200 HESD(J) = C1200*X(J)**2 - C400*X(J+1) + TWO HESL(IX(J+1,J)) = -C400*X(J) 500 CONTINUE RETURN C C EXTENDED POWELL FUNCTION. C 510 CONTINUE DO 520 J = 1,N* (N-1)/2 HESL(J) = ZERO 520 CONTINUE DO 530 J = 1,N,4 T2 = X(J+1) - TWO*X(J+2) T3 = X(J) - X(J+3) S1 = C12*T2**2 S2 = C120*T3**2 HESD(J) = TWO + S2 HESD(J+1) = C200 + S1 HESD(J+2) = TEN + FOUR*S1 HESD(J+3) = TEN + S2 HESL(IX(J+1,J)) = TWO*TEN HESL(IX(J+2,J)) = ZERO HESL(IX(J+2,J+1)) = -TWO*S1 HESL(IX(J+3,J)) = -S2 HESL(IX(J+3,J+1)) = ZERO HESL(IX(J+3,J+2)) = -TEN 530 CONTINUE RETURN C C BEALE FUNCTION. C 540 CONTINUE S1 = ONE - X(2) T1 = C1P5 - X(1)*S1 S2 = ONE - X(2)**2 T2 = C2P25 - X(1)*S2 S3 = ONE - X(2)**3 T3 = C2P625 - X(1)*S3 HESD(1) = TWO* (S1**2+S2**2+S3**2) HESD(2) = TWO*X(1)* (X(1)+TWO*T2+FOUR*X(1)*X(2)**2+SIX*X(2)*T3+ + NINE*X(1)*X(2)**4) HESL(1) = TWO* (T1-X(1)*S1) + FOUR*X(2)* (T2-X(1)*S2) + + SIX* (T3-X(1)*S3)*X(2)**2 RETURN C C WOOD FUNCTION. C 550 CONTINUE HESD(1) = C1200*X(1)**2 - C400*X(2) + TWO HESD(2) = C220P2 HESD(3) = C1080*X(3)**2 - C360*X(4) + TWO HESD(4) = C200P2 HESL(1) = -C400*X(1) HESL(2) = ZERO HESL(3) = ZERO HESL(4) = ZERO HESL(5) = C19P8 HESL(6) = -C360*X(3) RETURN C C CHEBYQUAD FUNCTION. C 560 CONTINUE DO 570 I = 1,N FVEC(I) = ZERO 570 CONTINUE DO 590 J = 1,N T1 = ONE T2 = TWO*X(J) - ONE T = TWO*T2 DO 580 I = 1,N FVEC(I) = FVEC(I) + T2 TH = T*T2 - T1 T1 = T2 T2 = TH 580 CONTINUE 590 CONTINUE D1 = ONE/FLOAT(N) IEV = .FALSE. DO 600 I = 1,N FVEC(I) = D1*FVEC(I) IF (IEV) FVEC(I) = FVEC(I) + ONE/ (DFLOAT(I)**2-ONE) IEV = .NOT. IEV 600 CONTINUE D2 = TWO*D1 M = 0 DO 640 J = 1,N HESD(J) = FOUR*D1 T1 = ONE T2 = TWO*X(J) - ONE T = TWO*T2 S1 = ZERO S2 = TWO P1 = ZERO P2 = ZERO GVEC(1) = S2 DO 610 I = 2,N TH = FOUR*T2 + T*S2 - S1 S1 = S2 S2 = TH TH = T*T2 - T1 T1 = T2 T2 = TH TH = EIGHT*S1 + T*P2 - P1 P1 = P2 P2 = TH GVEC(I) = S2 HESD(J) = HESD(J) + FVEC(I)*TH + D1*S2**2 610 CONTINUE HESD(J) = D2*HESD(J) DO 630 K = 1,J - 1 M = M + 1 HESL(M) = ZERO TT1 = ONE TT2 = TWO*X(K) - ONE TT = TWO*TT2 SS1 = ZERO SS2 = TWO DO 620 I = 1,N HESL(M) = HESL(M) + SS2*GVEC(I) TTH = FOUR*TT2 + TT*SS2 - SS1 SS1 = SS2 SS2 = TTH TTH = TT*TT2 - TT1 TT1 = TT2 TT2 = TTH 620 CONTINUE HESL(M) = D2*D1*HESL(M) 630 CONTINUE 640 CONTINUE RETURN C C LAST CARD OF SUBROUTINE HESFCN. C END SHAR_EOF fi # end of overwriting check if test -f 'tn1.f' then echo shar: will not over-write existing file "'tn1.f'" else cat << SHAR_EOF > 'tn1.f' C ---------------------------------------------------------- C SEGMENT 1: C 1. Driver and Routines for Testing TNPACK-Minimization C for the Extended Rosenbrock Function C C********************************************************************* PROGRAM MINROS C ----------------------------------------------------------------- C TNPACK-Minimization Driver for the Extended Rosenbrock Function C of dimension N (N even) C ----------------------------------------------------------------- C .. Parameters .. INTEGER MAXN PARAMETER (MAXN=1000) INTEGER N,NZ,MP,LW,LIW PARAMETER (N=1000,NZ=N,MP=6,LW=10*N+5*NZ,LIW=7*N+5*NZ+1) C .. C .. Scalars in Common .. INTEGER NPROB C .. C .. Arrays in Common .. DOUBLE PRECISION HESD(MAXN),HESOD(MAXN) C .. C .. Local Scalars .. DOUBLE PRECISION F,RJ REAL TIMEED,TIMEST,TIMEUR INTEGER INFORM,J C .. C .. Local Arrays .. DOUBLE PRECISION G(N),PLIST(20),W(LW),X(N) REAL TARRAY(2) INTEGER IW(LIW),OPLIST(20) C .. C .. External Subroutines .. EXTERNAL ROSFUN,ROSHDP,ROSPAT,SETLIS,TNMIN C .. C .. Intrinsic Functions .. INTRINSIC COS C .. C .. Common blocks .. C .. CPU time COMMON NPROB COMMON /HESIAN/HESD,HESOD C .. C .. External Functions .. REAL ETIME EXTERNAL ETIME C .. NPROB = 1 IF (N.GT.MAXN) THEN WRITE (MP,FMT=9000) MAXN GO TO 20 END IF c C--- initial guess c DO 10 J = 1,N - 1,2 RJ = J X(J) = -1.2D0 - COS(RJ) X(J+1) = 1.D0 + COS(RJ) 10 CONTINUE c c--- Sets sample parameters for minimization c CALL SETLIS(N,OPLIST,PLIST,INFORM) C --------------------------------------------- C Modify OPLIST and PLIST elements as desired: C OPLIST(1) - Preconditioning option C OPLIST(5) - Max. PCG itns. C OPLIST(10) - Hd mult. option C PLIST(2) - Convergence parameter EPSG C PLIST(3) - Residual truncation parameter c Four new elements: c OPLIST(15) - Option selector for termination rule of PCG C 1: the modified negative curvature test; C 2: the descent direction test C Default value is 1 C OPLIST(16) - Option selector for line search stopping rule C 1: Criterion 1 C 2: Criterion 2 (more lenient stopping rule) C Default value is 1 C OPLIST(17) - option selector for modified Cholesky factorization C 1: the MC by Gill and Murray C 2: the UMC by Xie and Schlick C Default value is 1 C PLIST(8) - Parameter tau of UMC. Default value is 10.d0 C --------------------------------------------- OPLIST(1) = 1 OPLIST(5) = 100 OPLIST(10) = 0 OPLIST(12) = 1 PLIST(2) = 1.0D-08 PLIST(3) = 0.1 c--- xie OPLIST(15) = 1 OPLIST(16) = 2 OPLIST(17) = 2 PLIST(8) = 10.0D0 c------------------------- c c--- Truncated Newton Interface Routine c TIMEST = ETIME(TARRAY) TIMEUR = TARRAY(1) CALL TNMIN(N,X,F,G,OPLIST,PLIST,INFORM,NZ,W,LW,IW,LIW,ROSFUN, + ROSPAT,ROSHDP) TIMEED = ETIME(TARRAY) - TIMEST TIMEUR = TARRAY(1) - TIMEUR WRITE (*,FMT=*) ' The CPU time in seconds ',TIMEUR 20 CONTINUE STOP 9000 FORMAT (/,5X,' N > MAXN, MAXN = ',I6,/,6X, + 'reset PARAMETER statement for MAXN',/) END C********************************************************************* SUBROUTINE ROSFUN(N,X,F,G,A,IA,JA,NOUT) C ------------------------------------------------------- C ROSENBROCK'S FUNCTION OF DIMENSION N - ASSEMBLE F,G,H,M C ------------------------------------------------------- c c--- F: ROSENBROCK'S FUNCTION c C .. Parameters .. INTEGER MAXN PARAMETER (MAXN=1000) C .. C .. Scalar Arguments .. DOUBLE PRECISION F INTEGER N,NOUT C .. C .. Array Arguments .. DOUBLE PRECISION A(*),G(N),X(N) INTEGER IA(N+1),JA(*) C .. C .. Arrays in Common .. DOUBLE PRECISION HESD(MAXN),HESOD(MAXN) C .. C .. Local Scalars .. DOUBLE PRECISION T1,T2 INTEGER J C .. C .. Common blocks .. COMMON /HESIAN/HESD,HESOD C .. F = 0.D0 DO 10 J = 1,N - 1,2 T1 = 1.D0 - X(J) T2 = 10.D0* (X(J+1)-X(J)*X(J)) F = F + T1*T1 + T2*T2 10 CONTINUE IF (NOUT.EQ.0) RETURN c c--- Computing the componets G(j) of the gradient of F c DO 20 J = 1,N - 1,2 T1 = 1.D0 - X(J) G(J+1) = 200.D0* (X(J+1)-X(J)*X(J)) G(J) = -2.D0* (X(J)*G(J+1)+T1) 20 CONTINUE IF (NOUT.EQ.1) RETURN C ------------------------------------------------------ c Computing Hessian matrix C H is stored in 2 arrays: HESD for diagonals, HESOD for C off-diagonals, stored by rows for the upper triangle. C HESOD is zero every even entry since H has the 2x2 C block-diagonal pattern: C C * * C * * C * * C * * C Etc. C C M is set to the diagonal of H. c c H_jj = HESD( J ) H_{j,j+1} = HESOD(J) c H_{j+1,j+1} = HESD(J+1) C ------------------------------------------------------ DO 30 J = 1,N - 1,2 HESD(J+1) = 200.D0 HESD(J) = 2.D0 - 400.D0*X(J+1) + 1200.D0*X(J)*X(J) HESOD(J) = -400.D0*X(J) 30 CONTINUE IF (NOUT.EQ.2) RETURN c c--- Computing the preconditor M: the diagonal of H c DO 40 J = 1,N - 1,2 A(J) = -2.D0* (200.D0*X(J+1)-1.D0-600.D0* (X(J)**2)) A(J+1) = 200.0D0 40 CONTINUE RETURN END C*********************************************************************** SUBROUTINE ROSHDP(N,D,HD,Y1,Y2) C ------------------------------------------------- C COMPUTE c HD = H*D c where H is the Hessian matrix and D is a vector. C ------------------------------------------------- C .. Parameters .. INTEGER MAXN PARAMETER (MAXN=1000) C .. C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION D(N),HD(N),Y1(N),Y2(N) C .. C .. Arrays in Common .. DOUBLE PRECISION HESD(MAXN),HESOD(MAXN) C .. C .. Local Scalars .. DOUBLE PRECISION OFDIAG INTEGER J C .. C .. Common blocks .. COMMON /HESIAN/HESD,HESOD C .. DO 10 J = 1,N - 1,2 OFDIAG = HESOD(J) HD(J) = HESD(J)*D(J) + OFDIAG*D(J+1) HD(J+1) = OFDIAG*D(J) + HESD(J+1)*D(J+1) 10 CONTINUE RETURN END C*********************************************************************** SUBROUTINE ROSPAT(N,X,A,IA,JA) C ------------------------------------------------------------- C DETERMINE The PATTERN OF M C ------------------------------------------------------------- C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION A(*),X(N) INTEGER IA(N+1),JA(*) C .. C .. Local Scalars .. INTEGER I C .. DO 10 I = 1,N IA(I) = I JA(I) = I 10 CONTINUE IA(N+1) = N + 1 RETURN END SHAR_EOF fi # end of overwriting check if test -f 'tn2.f' then echo shar: will not over-write existing file "'tn2.f'" else cat << SHAR_EOF > 'tn2.f' * -------------------------------------------------------------------- * * SEGMENT 2: * 2. Driver and Routines for Testing TNPACK-Minimization * for the Trigonometric Function. * * To test the re-ordering option in TNPACK, a preconditioner * M = ( m_{ij} ), a matrix of order n, is defined by * * m_{i,j} = H_{ii} if i = j * 0.1 if i = 1, j = n-1 * -0.1 if i = 1, j = n * 0.1 if i = n-1, j = 1 * -0.1 if i = n, j = 1 * 0.0 otherwise. * * where H_{ii} is the diagonal element of Hessian matrix. * ************************************************************************ PROGRAM MINTRG * ----------------------------------------------------------------- * TNPACK-MinimizatiON dRIVer for the Trigonometric Function * ----------------------------------------------------------------- C .. Parameters .. INTEGER MAXM,MAXLH PARAMETER (MAXM=5000,MAXLH=MAXM* (MAXM-1)/2) INTEGER N,NZ,MP,LW,LIW PARAMETER (N=1000,NZ=3*N,MP=6,LW=9*N+5*NZ+7*N,LIW=7*N+5*NZ+1) C .. C .. Scalars in Common .. INTEGER NPROB C .. C .. Arrays in Common .. DOUBLE PRECISION HESD(MAXM),HESL(MAXLH) C .. C .. Local Scalars .. DOUBLE PRECISION F,RI,XZERO REAL TIMEED,TIMEST,TIMEUR INTEGER I,INFORM C .. C .. Local Arrays .. DOUBLE PRECISION G(N),PLIST(20),W(LW),X(N) REAL TARRAY(2) INTEGER IW(LIW),OPLIST(20) C .. C .. External Subroutines .. EXTERNAL SETLIS,TNMIN,TRGFUN,TRGHDP,TRGPAT C .. C .. Intrinsic Functions .. INTRINSIC COS,DBLE C .. C .. Common blocks .. COMMON NPROB COMMON /HESIAN/HESL,HESD C .. C .. External Functions .. REAL ETIME EXTERNAL ETIME C .. NPROB = 1 IF (N.GT.MAXM) THEN WRITE (MP,FMT=9000) MAXM GO TO 20 END IF XZERO = 1.D0/DBLE(N) DO 10 I = 1,N RI = I X(I) = XZERO + 0.2D0*COS(RI) 10 CONTINUE CALL SETLIS(N,OPLIST,PLIST,INFORM) * --------------------------------------------- * Modify OPLIST elements as desired: * OPLIST(1) - Preconditioning option * OPLIST(4) - Max. Newton itns. * OPLIST(5) - Max. PCG itns. * OPLIST(6) - Max. F&G evals. * OPLIST(7) - 1: reoder; 0: do not reorder * OPLIST(10) - Hd mult. option * Four new elements: * OPLIST(15) - 1: Test 1A' 2: Test 2A * Default value is 1 * OPLIST(16) - Option selector for line search stopping rule * 1: Criterion 1 * 2: Criterion 2 (more lenient stopping rule) * Default value is 1 * OPLIST(17) - Option selector for modified Cholesky factorization * 1: the MC by Gill and Murray * 2: the UMC by Xie and Schlick * Default value is 1 * PLIST(8) - UMC parameter tau * --------------------------------------------- OPLIST(1) = 1 OPLIST(5) = 40 OPLIST(6) = 10*N OPLIST(7) = 0 OPLIST(10) = 0 OPLIST(12) = 1 c--- xie OPLIST(15) = 2 OPLIST(16) = 2 OPLIST(17) = 2 PLIST(8) = 0.5D0 c------------------------- TIMEST = ETIME(TARRAY) TIMEUR = TARRAY(1) CALL TNMIN(N,X,F,G,OPLIST,PLIST,INFORM,NZ,W,LW,IW,LIW,TRGFUN, + TRGPAT,TRGHDP) TIMEED = ETIME(TARRAY) - TIMEST TIMEUR = TARRAY(1) - TIMEUR WRITE (*,FMT=*) ' The CPU time in seconds ',TIMEUR 20 CONTINUE STOP 9000 FORMAT (/,5X,' N > MAXM, MAXM = ',I6,/,6X, + 'reset PARAMETER statement for MAXM',/) END ************************************************************************ SUBROUTINE TRGFUN(N,X,F,G,A,IA,JA,NOUT) * -------------------------------------------------------- * TRIGONOMETRIC FUNCTION OF DIMENSION N - ASSEMBLE F,G,H,M * -------------------------------------------------------- C .. Parameters .. INTEGER MAXM,MAXLH PARAMETER (MAXM=5000,MAXLH=MAXM* (MAXM-1)/2) C .. C .. Scalar Arguments .. DOUBLE PRECISION F INTEGER N,NOUT C .. C .. Array Arguments .. DOUBLE PRECISION A(*),G(N),X(N) INTEGER IA(*),JA(*) C .. C .. Arrays in Common .. DOUBLE PRECISION HESD(MAXM),HESL(MAXLH) C .. C .. Local Scalars .. DOUBLE PRECISION SUMCOS,SUMTRM,TWO,ZERO INTEGER I,II,J,JJ,K C .. C .. Local Arrays .. DOUBLE PRECISION COSINE(MAXM),SINE(MAXM),TERM(MAXM) C .. C .. Intrinsic Functions .. INTRINSIC COS,DBLE,SIN C .. C .. Common blocks .. COMMON /HESIAN/HESL,HESD C .. C .. Statement Functions .. INTEGER IX C .. C .. Save statement .. SAVE /HESIAN/ C .. C .. Statement Function definitions .. IX(II,JJ) = ((II-1)* (II-2)/2) + JJ C .. ZERO = 0.D0 TWO = 2.D0 SUMCOS = ZERO DO 10 I = 1,N SINE(I) = SIN(X(I)) COSINE(I) = COS(X(I)) SUMCOS = SUMCOS + COSINE(I) 10 CONTINUE F = ZERO DO 20 I = 1,N TERM(I) = DBLE(N+I) - SUMCOS - DBLE(I)*COSINE(I) - SINE(I) F = F + TERM(I)**2 20 CONTINUE IF (NOUT.EQ.0) RETURN SUMTRM = ZERO DO 30 I = 1,N SUMTRM = SUMTRM + TERM(I) 30 CONTINUE DO 40 I = 1,N G(I) = SINE(I)*SUMTRM + (DBLE(I)*SINE(I)-COSINE(I))*TERM(I) G(I) = TWO*G(I) 40 CONTINUE IF (NOUT.EQ.1) RETURN DO 50 I = 1,N HESD(I) = DBLE(I* (I+2)+N)*SINE(I)**2 + + COSINE(I)* (SUMTRM+COSINE(I)- + (DBLE(2*I+2)*SINE(I))) + + TERM(I)* (DBLE(I)*COSINE(I)+SINE(I)) HESD(I) = TWO*HESD(I) 50 CONTINUE DO 70 I = 2,N DO 60 J = 1,I - 1 K = IX(I,J) HESL(K) = DBLE(N+I+J)*SINE(I)*SINE(J) - + SINE(I)*COSINE(J) - SINE(J)*COSINE(I) HESL(K) = TWO*HESL(K) 60 CONTINUE 70 CONTINUE IF (NOUT.EQ.2) RETURN DO I = 4,N + 2 A(I) = HESD(I-2) END DO A(1) = HESD(1) A(2) = 0.1D0 A(3) = -0.1D0 RETURN END ************************************************************************ SUBROUTINE TRGHDP(N,D,HD,XC,GC) * ------------------------------------------------------ * TRIGONOMETRIC FUNCTION OF DIMENSION N - COMPUTE HD * ------------------------------------------------------ C .. Parameters .. INTEGER MAXM,MAXLH PARAMETER (MAXM=5000,MAXLH=MAXM* (MAXM-1)/2) C .. C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION D(N),GC(N),HD(N),XC(N) C .. C .. Arrays in Common .. DOUBLE PRECISION HESD(MAXM),HESL(MAXLH) C .. C .. Local Scalars .. INTEGER I,II,J,JJ C .. C .. Common blocks .. COMMON /HESIAN/HESL,HESD C .. C .. Statement Functions .. INTEGER IX C .. C .. Save statement .. SAVE /HESIAN/ C .. C .. Statement Function definitions .. IX(II,JJ) = ((II-1)* (II-2)/2) + JJ C .. DO 10 I = 1,N HD(I) = HESD(I)*D(I) 10 CONTINUE DO 40 I = 1,N DO 20 J = 1,I - 1 HD(I) = HD(I) + HESL(IX(I,J))*D(J) 20 CONTINUE DO 30 J = I + 1,N HD(I) = HD(I) + HESL(IX(J,I))*D(J) 30 CONTINUE 40 CONTINUE RETURN END ************************************************************************ SUBROUTINE TRGPAT(N,X,A,IA,JA) * -------------------------------------------------------------- * TRIGONOMETRIC FUNCTION OF DIMENSION N - DETERMINE PATTERN OF M * -------------------------------------------------------------- C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION A(*),X(*) INTEGER IA(N+1),JA(*) C .. C .. Local Scalars .. INTEGER I C .. DO 10 I = 1,N IA(I) = I 10 CONTINUE DO I = 4,N + 2 JA(I) = I - 2 END DO JA(1) = 1 JA(2) = N - 1 JA(3) = N DO I = 2,N IA(I) = IA(I) + 2 END DO IA(N+1) = IA(N) + 1 RETURN END SHAR_EOF fi # end of overwriting check if test -f 'tn3.f' then echo shar: will not over-write existing file "'tn3.f'" else cat << SHAR_EOF > 'tn3.f' * ---------------------------------------------------------- * * SEGMENT 3: * 3. Driver and Routines for Testing TNPACK-Minimization * for the Functions defined in Alg. 566 * ************************************************************************ PROGRAM MINFUC566 * ----------------------------------------------------------------- * TNPACK-Minimization Driver for the Functions in Alg.566 * ----------------------------------------------------------------- C .. Parameters .. INTEGER MAXM,MAXLH PARAMETER (MAXM=5000,MAXLH=MAXM* (MAXM-1)/2) INTEGER NN,NZ,MP,LW,LIW PARAMETER (NN=1000,NZ=3*NN,MP=6,LW=10*NN+5*NZ,LIW=7*NN+5*NZ+1) C .. C .. Scalars in Common .. INTEGER NPROB C .. C .. Arrays in Common .. DOUBLE PRECISION HESD(MAXM),HESL(MAXLH) C .. C .. Local Scalars .. DOUBLE PRECISION F,FACTOR INTEGER INFORM,N C .. C .. Local Arrays .. DOUBLE PRECISION G(NN),PLIST(20),W(LW),X(NN) INTEGER IW(LIW),OPLIST(20) C .. C .. External Subroutines .. EXTERNAL INITPT,PROBSET,SETLIS,TNMIN,TRGFUN,TRGHDP,TRGPAT C .. C .. Common blocks .. COMMON NPROB COMMON /HESIAN/HESL,HESD C .. NPROB = 0 10 CONTINUE NPROB = NPROB + 1 WRITE (MP,FMT=9010) WRITE (MP,FMT=9000) NPROB WRITE (MP,FMT=9010) * *--- Enter the value of N * CALL PROBSET(NPROB,N) IF (N.GT.MAXM) THEN WRITE (MP,FMT=9020) MAXM GO TO 20 END IF * *--- Initial guess * FACTOR = 1.0 CALL INITPT(N,X,NPROB,FACTOR) CALL SETLIS(N,OPLIST,PLIST,INFORM) OPLIST(1) = 1 OPLIST(4) = 10000 OPLIST(5) = 100 OPLIST(6) = 4000 OPLIST(10) = 0 OPLIST(13) = 0 PLIST(2) = 1.D-8 PLIST(3) = 1.0D0 *-- xie OPLIST(15) = 2 OPLIST(16) = 2 OPLIST(17) = 2 PLIST(8) = 10.D0 *------------------------------ CALL TNMIN(N,X,F,G,OPLIST,PLIST,INFORM,NZ,W,LW,IW,LIW,TRGFUN, + TRGPAT,TRGHDP) 20 CONTINUE IF (NPROB.LT.18) GO TO 10 STOP 9000 FORMAT (/,20X,'Problem ',I2) 9010 FORMAT (/,14X,'=======================') 9020 FORMAT (/,5X,' N > MAXM, MAXM = ',I6,/,6X, + 'reset PARAMETER statement for MAXM',/) END ************************************************************************ SUBROUTINE TRGFUN(N,X,F,G,A,IA,JA,NOUT) * -------------------------------------------------------- * ASSEMBLE F,G,H,M * -------------------------------------------------------- C .. Parameters .. INTEGER MAXM,MAXLH PARAMETER (MAXM=5000,MAXLH=MAXM* (MAXM-1)/2) INTEGER LBAND PARAMETER (LBAND=1) C .. C .. Scalar Arguments .. DOUBLE PRECISION F INTEGER N,NOUT C .. C .. Array Arguments .. DOUBLE PRECISION A(*),G(N),X(N) INTEGER IA(*),JA(*) C .. C .. Scalars in Common .. INTEGER NPROB C .. C .. Arrays in Common .. DOUBLE PRECISION HESD(MAXM),HESL(MAXLH) C .. C .. Local Scalars .. INTEGER I,II,J,JJ,K,M C .. C .. External Subroutines .. EXTERNAL GRDFCN,HESFCN,OBJFCN C .. C .. Intrinsic Functions .. INTRINSIC MIN C .. C .. Common blocks .. COMMON NPROB COMMON /HESIAN/HESL,HESD C .. C .. Statement Functions .. INTEGER IX C .. C .. Save statement .. SAVE /HESIAN/ C .. C .. Statement Function definitions .. IX(II,JJ) = ((II-1)* (II-2)/2) + JJ C .. *--- THIS SUBROUTINE DEFINES THE OBJECTIVE FUNCTIONS OF EIGHTEEN * NONLINEAR UNCONSTRAINED MINIMIZATION PROBLEMS. * THE SUBROUTINE STATEMENT IS * * SUBROUTINE OBJFCN(N,X,F,NPROB) * * WHERE * * N IS A POSITIVE INTEGER INPUT VARIABLE. * * X IS AN INPUT ARRAY OF LENGTH N. * * F IS AN OUTPUT VARIABLE WHICH CONTAINS THE VALUE OF * THE NPROB OBJECTIVE FUNCTION EVALUATED AT X. * * NPROB IS A POSITIVE INTEGER INPUT VARIABLE WHICH DEFINES THE * NUMBER OF THE PROBLEM. NPROB MUST NOT EXCEED 18. CALL OBJFCN(N,X,F,NPROB) IF (NOUT.EQ.0) RETURN *--- THIS SUBROUTINE DEFINES THE GRADIENT VECTORS OF EIGHTEEN * NONLINEAR UNCONSTRAINED MINIMIZATION PROBLEMS. * * N IS A POSITIVE INTEGER INPUT VARIABLE. * * X IS AN INPUT ARRAY OF LENGTH N. * * G IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE COMPONENTS * OF THE GRADIENT VECTOR OF THE NPROB OBJECTIVE FUNCTION * EVALUATED AT X. * * NPROB IS A POSITIVE INTEGER INPUT VARIABLE WHICH DEFINES THE * NUMBER OF THE PROBLEM. NPROB MUST NOT EXCEED 18. *----------------------------------------------------------------------- CALL GRDFCN(N,X,G,NPROB) IF (NOUT.EQ.1) RETURN *--- THIS SUBROUTINE DEFINES THE HESSIAN MATRICES OF 18 * NONLINEAR UNCONSTRAINED MINIMIZATION PROBLEMS. * * N IS A POSITIVE INTEGER INPUT VARIABLE. * * X IS AN INPUT ARRAY OF LENGTH N. * * HESD IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE * DIAGONAL COMPONENTS OF THE HESSIAN MATRIX OF THE NPROB * OBJECTIVE FUNCTION EVALUATED AT X. * * HESL IS AN OUTPUT ARRAY OF LENGTH N*(N-1)/2 WHICH CONTAINS * THE LOWER TRIANGULAR PART OF THE HESSIAN MATRIX OF THE * NPROB OBJECTIVE FUNCTION EVALUATED AT X. * * NPROB IS A POSITIVE INTEGER INPUT VARIABLE WHICH DEFINES THE * NUMBER OF THE PROBLEM. NPROB MUST NOT EXCEED 18. *-------------------------------------------------------------------- CALL HESFCN(N,X,HESD,HESL,NPROB) IF (NOUT.EQ.2) RETURN K = 1 DO 20 I = 1,N A(K) = HESD(I) M = MIN(LBAND/2,N-I) IF (M.GT.0) THEN DO 10 J = 1,M A(K+J) = HESL(IX(I+J,I)) 10 CONTINUE END IF K = K + M + 1 20 CONTINUE RETURN END ************************************************************************ SUBROUTINE TRGHDP(N,D,HD,XC,GC) * ------------------------------------------------------ * COMPUTE HD * ------------------------------------------------------ C .. Parameters .. INTEGER MAXM,MAXLH PARAMETER (MAXM=5000,MAXLH=MAXM* (MAXM-1)/2) C .. C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION D(N),GC(N),HD(N),XC(N) C .. C .. Arrays in Common .. DOUBLE PRECISION HESD(MAXM),HESL(MAXLH) C .. C .. Local Scalars .. INTEGER I,II,J,JJ C .. C .. Common blocks .. COMMON /HESIAN/HESL,HESD C .. C .. Statement Functions .. INTEGER IX C .. C .. Save statement .. SAVE /HESIAN/ C .. C .. Statement Function definitions .. IX(II,JJ) = ((II-1)* (II-2)/2) + JJ C .. DO 10 I = 1,N HD(I) = HESD(I)*D(I) 10 CONTINUE DO 40 I = 1,N DO 20 J = 1,I - 1 HD(I) = HD(I) + HESL(IX(I,J))*D(J) 20 CONTINUE DO 30 J = I + 1,N HD(I) = HD(I) + HESL(IX(J,I))*D(J) 30 CONTINUE 40 CONTINUE RETURN END ************************************************************************ SUBROUTINE TRGPAT(N,X,A,IA,JA) * -------------------------------------------------------------- * DETERMINE PATTERN OF M * -------------------------------------------------------------- C .. Parameters .. INTEGER LBAND PARAMETER (LBAND=1) C .. C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION A(*),X(*) INTEGER IA(N+1),JA(*) C .. C .. Local Scalars .. INTEGER I,J,K,M C .. C .. Intrinsic Functions .. INTRINSIC MIN C .. K = 1 DO 20 I = 1,N IA(I) = K M = MIN(LBAND/2,N-I) DO 10 J = 0,M JA(K+J) = I + J 10 CONTINUE K = K + M + 1 20 CONTINUE IA(N+1) = K RETURN END SUBROUTINE PROBSET(NPROB,N) C .. Scalar Arguments .. INTEGER N,NPROB C .. C .. Local Arrays .. INTEGER NVARS(18) C .. C .. Intrinsic Functions .. INTRINSIC MAX,MIN,MOD C .. C .. Data statements .. DATA NVARS/3,6,3,2,3,0,0,0,0,2,4,3,0,0,0,2,4,0/ C .. IF (NPROB.EQ.0) STOP IF (NPROB.LT.1 .OR. NPROB.GT.18) THEN WRITE (*,FMT=9000) STOP END IF IF (NVARS(NPROB).NE.0) N = NVARS(NPROB) IF (NPROB.EQ.7) N = MAX(2,MIN(N,31)) IF (NPROB.EQ.14) N = MAX(2,N-MOD(N,2)) IF (NPROB.EQ.15) N = MAX(4,N-MOD(N,4)) IF (NPROB.EQ.18) N = MAX(2,MIN(N,50)) RETURN 9000 FORMAT (/,4X,'ERROR IN INPUT FILE',/) END SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'blas1.f' then echo shar: will not over-write existing file "'blas1.f'" else cat << SHAR_EOF > 'blas1.f' C*********************************************************************** C*********************************************************************** C BLAS LEVEL 1, DOUBLE PRECISION C FROM NETLIB, FRI OCT 12 17:03:00 EDT 1990 C*********************************************************** SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) C -------------------------------------------------------- C CONSTANT TIMES A VECTOR PLUS A VECTOR. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C -------------------------------------------------------- DOUBLE PRECISION DX(*),DY(*),DA INTEGER I,INCX,INCY,IX,IY,M,MP1,N IF (N. LE. 0) RETURN IF (DA .EQ. 0.0D0) RETURN IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 C -------------------------------------------------------- C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C -------------------------------------------------------- IX = 1 IY = 1 IF (INCX .LT. 0) IX = (-N+1)*INCX + 1 IF (INCY .LT. 0) IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DY(IY) + DA*DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C -------------------------------------------------------- C CODE FOR BOTH INCREMENTS EQUAL TO 1 C CLEAN-UP LOOP C -------------------------------------------------------- 20 M = MOD(N,4) IF (M .EQ. 0) GO TO 40 DO 30 I = 1,M DY(I) = DY(I) + DA*DX(I) 30 CONTINUE IF (N .LT. 4) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 DY(I) = DY(I) + DA*DX(I) DY(I + 1) = DY(I + 1) + DA*DX(I + 1) DY(I + 2) = DY(I + 2) + DA*DX(I + 2) DY(I + 3) = DY(I + 3) + DA*DX(I + 3) 50 CONTINUE RETURN END C*********************************************************************** C*********************************************************************** SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) C -------------------------------------------------------- C COPIES A VECTOR, X, TO A VECTOR, Y. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C -------------------------------------------------------- DOUBLE PRECISION DX(*),DY(*) INTEGER I,INCX,INCY,IX,IY,M,MP1,N IF (N .LE. 0) RETURN IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 C -------------------------------------------------------- C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C -------------------------------------------------------- IX = 1 IY = 1 IF (INCX .LT. 0) IX = (-N+1)*INCX + 1 IF (INCY .LT. 0) IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C -------------------------------------------------------- C CODE FOR BOTH INCREMENTS EQUAL TO 1 C CLEAN-UP LOOP C -------------------------------------------------------- 20 M = MOD(N,7) IF (M .EQ. 0) GO TO 40 DO 30 I = 1,M DY(I) = DX(I) 30 CONTINUE IF (N .LT. 7) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 DY(I) = DX(I) DY(I + 1) = DX(I + 1) DY(I + 2) = DX(I + 2) DY(I + 3) = DX(I + 3) DY(I + 4) = DX(I + 4) DY(I + 5) = DX(I + 5) DY(I + 6) = DX(I + 6) 50 CONTINUE RETURN END C*********************************************************************** C*********************************************************************** DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) C -------------------------------------------------------- C FORMS THE DOT PRODUCT OF TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C -------------------------------------------------------- DOUBLE PRECISION DX(*),DY(*),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N DDOT = 0.0D0 DTEMP = 0.0D0 IF (N .LE. 0) RETURN IF (INCX.EQ.1. AND. INCY.EQ.1) GO TO 20 C -------------------------------------------------------- C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C -------------------------------------------------------- IX = 1 IY = 1 IF (INCX .LT. 0) IX = (-N+1)*INCX + 1 IF (INCY .LT. 0) IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP = DTEMP + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE DDOT = DTEMP RETURN C -------------------------------------------------------- C CODE FOR BOTH INCREMENTS EQUAL TO 1 C CLEAN-UP LOOP C -------------------------------------------------------- 20 M = MOD(N,5) IF (M .EQ. 0) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DX(I)*DY(I) 30 CONTINUE IF (N .LT. 5) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,5 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4) 50 CONTINUE 60 DDOT = DTEMP RETURN END C*********************************************************************** C*********************************************************************** DOUBLE PRECISION FUNCTION DNRM2 (N,DX,INCX) INTEGER NEXT DOUBLE PRECISION DX(*),CUTLO,CUTHI,HITEST,SUM,XMAX,ZERO,ONE DATA ZERO,ONE /0.0D0,1.0D0/ C -------------------------------------------------------- C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() 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.L.LAWSON, 1978 JAN 08 C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE C HOPEFULLY APPLICABLE TO ALL MACHINES. C CUTLO = MAXIMUM OF DSQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF DSQRT(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 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 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 / INTEGER N, NN, INCX, i, J DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / IF (N .GT. 0) GO TO 10 DNRM2 = ZERO GO TO 300 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N * INCX C ------------------ C BEGIN MAIN LOOP C ------------------ I = 1 20 GO TO NEXT,(30, 50, 70, 110) 30 IF (DABS(DX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C ------------------ C PHASE 1. SUM IS ZERO C ------------------ 50 IF (DX(I) .EQ. ZERO) GO TO 200 IF (DABS(DX(I)) .GT. CUTLO) GO TO 85 C ------------------ C PREPARE FOR PHASE 2. C ------------------ ASSIGN 70 TO NEXT GO TO 105 C ------------------ C PREPARE FOR PHASE 4. C ------------------ 100 I = J ASSIGN 110 TO NEXT SUM = (SUM / DX(I)) / DX(I) 105 XMAX = DABS(DX(I)) GO TO 115 C ------------------ C SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C ------------------ 70 IF (DABS(DX(I)) .GT. CUTLO) GO TO 75 C ------------------ C COMMON CODE FOR PHASES 2 AND 4. C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. C ------------------ 110 IF (DABS(DX(I)) .LE. XMAX) GO TO 115 SUM = ONE + SUM * (XMAX / DX(I))**2 XMAX = DABS(DX(I)) GO TO 200 115 SUM = SUM + (DX(I)/XMAX)**2 GO TO 200 C ------------------ C PREPARE FOR PHASE 3. C ------------------ 75 SUM = (SUM * XMAX) * XMAX C ------------------ C FOR REAL OR D.P. SET HITEST = CUTHI/N C FOR COMPLEX SET HITEST = CUTHI/(2*N) C ------------------ 85 HITEST = CUTHI/FLOAT(N) C ------------------ C PHASE 3. SUM IS MID-RANGE. NO SCALING. C ------------------ DO 95 J =I,NN,INCX IF (DABS(DX(J)) .GE. HITEST) GO TO 100 95 SUM = SUM + DX(J)**2 DNRM2 = DSQRT(SUM) GO TO 300 200 CONTINUE I = I + INCX IF (I .LE. NN) GO TO 20 C ------------------ C END OF MAIN LOOP. C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. C ------------------ DNRM2 = XMAX * DSQRT(SUM) 300 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check if test -f 'tnpack.f' then echo shar: will not over-write existing file "'tnpack.f'" else cat << SHAR_EOF > 'tnpack.f' C ---------------------------------------------------------- C SEGMENT 4: C 4. TNPACK routines (new version) c---------------------------------------------------------- c C*********************************************************************** C*********************************************************************** C * C - "TNPACK" - * C * C A PACKAGE OF ROUTINES FOR A * C LARGE-SCALE TRUNCATED NEWTON ALGORITHM * C * C copyright (c) 1990 by * C * C Tamar Schlick & Aaron Fogelson * c * c Updated (November 1998) by * c * c Dexuan Xie & Tamar Schlick * C * C The changes made involve the following components: * C * C (a) Modified Cholesky factorization * C (b) Negative curvature test * C (c) Line search scheme * C * C For details, see * C * C (1) Xie, D. and Schlick, T. Remark on the Updated Truncated * C Newton Minimization Package, Algorithm 702, ACM Trans. Math * C Softw., in Press (1998) * C * C (2) Xie, D. and Schlick, T. Efficient implementation of the * C truncated-Newton algorithm for large-scale chemistry * C applications, SIAM J. Opt., in Press (1998) * C * C The manuscripts are also available from the authors. Please send * C requests by email to * C * C dexuan@cims.nyu.edu or schlick@nyu.edu * C * C*********************************************************************** C * C PLEASE NOTE: * C * C (i) Double Precision is used througout. * C (ii) Documentation for the program is given separately; here * C only brief comments are provided. * C (iii) It is essential that the user become familiar with all input * C options and parameters before using this package. Entry to * C TNPACK from the user's driver must first be done through a * C call to routine SETLIS. SETLIS sets sample input options and * C parameters, some of which the user may then modify to suit * C the problem. Routine TNMIN should then be called. * C (iv) For runs on Cray computers, the Yale Sparse Matrix Package * C parameter RATIO must be changed from 2 to 1 - routine SDRVMD * C (v) If running on VAX/VMS systems with built-in facilities * C for BLAS, you may need to compile this program with something * C like FOR/NOBLAS, since some standard BLAS are included here. * C * C*********************************************************************** C * C TNPACK COMPONENTS: * C * C (A) USER-SUPPLIED SUBROUTINES (generic names): * C ----------------------------------------- * C * C * CALFGH - calculates F,G,H and M (preconditioner) at X * C * CALHDP - calculates HD, a Hessian/vector product * C * CALPAT - determines the sparsity pattern of M * C * C (B) TRUNCATED NEWTON SUBROUTINES: * C ---------------------------- * C * C * TNMIN - the TN driver routine * C * OUTER - the outer Newton iteration * C * INNER - the inner PCG iteration * C * OURHDP - computes Hessian/vector products * C * SETLIS - sets sample parameters for minimization * C * CHKLIS - checks input options and parameters * C * DIVWRK - divides work space * C * MLINES, NEWSTP - perform the line search * C * C (C) YSMP (YALE SPARSE MATRIX PACKAGE) SUBROUTINES: * C --------------------------------------------- * C * C * SDRVMD - the driver for solving Mz=r * C * ODRV - the driver for finding M's reordering * C * MD,MDI,MDM,MDP,MDU - the minimum-degree reordering routines * C * SRO - prepares for symmetric reordering of M * C * SSF - performs symbolic factorization of M * C * SNFMOD - performs numerical factorization of M * C * SNS - performs numerical solution of Mz=r * C * C (D) OTHER SUBPROGRAMS: * C ---------------- * C * C * DCOPY, DAXPY (blas subroutines) * C * DNRM2, DDOT, DMACH (blas functions) * C * C*********************************************************************** C*********************************************************************** SUBROUTINE TNMIN(N,X,F,G,OPLIST,PLIST,INFORM,NZ,W,LW,IW,LIW, + CALFGH,CALPAT,CALHDP) C -------------------------------------------------------- C TNMIN: Truncated Newton Interface Routine C -------------------------------------------------------- C .. Scalar Arguments .. DOUBLE PRECISION F INTEGER INFORM,LIW,LW,N,NZ C .. C .. Array Arguments .. DOUBLE PRECISION G(N),PLIST(20),W(LW),X(N) INTEGER IW(LIW),OPLIST(20) C .. C .. Subroutine Arguments .. EXTERNAL CALFGH,CALHDP,CALPAT C .. C .. Scalars in Common .. DOUBLE PRECISION TAU INTEGER ICALLS,IXA,IXD,IXD1,IXHD,IXIA,IXIPER,IXISP,IXJA,IXP,IXPER, + IXRSP,IXY,IXZ,MC,MP,NPROB,SRLS C .. C .. Local Scalars .. INTEGER I,MARK,NSP C .. C .. External Subroutines .. EXTERNAL CHKLIS,DIVWRK,OUTER C .. C .. Intrinsic Functions .. INTRINSIC MAX C .. C .. Common blocks .. COMMON NPROB COMMON /TAUINPUT/TAU,SRLS,MC COMMON /TN001/MP,ICALLS COMMON /TN004/IXP,IXY,IXD,IXHD,IXZ,IXA,IXRSP,IXD1 COMMON /TN005/IXPER,IXIPER,IXIA,IXJA,IXISP C .. C .. Save statement .. C -------------------------------------------------------- C check parameters and divide work space C -------------------------------------------------------- c c--- INFORM: indicates the status of TNPACK c SAVE /TN001/,/TN004/,/TN005/ C .. IF (INFORM.NE.-10) THEN WRITE (*,FMT=9000) INFORM = -1 RETURN ELSE INFORM = 0 END IF MARK = 0 c c--- Check input options & parameters (OPLIST,PLIST) c c INFORM = -2: input errors in N, OPLIST,PLIST c CALL CHKLIS(N,OPLIST,PLIST,MARK) IF (MARK.LT.0) THEN INFORM = -2 RETURN END IF c c--- Specifies the unit number for printing c MP = OPLIST(11) IF (ICALLS.EQ.0) THEN c--- xie: NPROB is used only for testing Alg.566 to print out c the input information only once for running all 18 problems. C To do so, uncomment "if (NPROB .eq. 1) then" and "endif" below c Then compile with make test3'' c--------------------------------------------------------------------- IF (NPROB.EQ.1) THEN WRITE (MP,FMT=9010) N, (OPLIST(I),I=1,6) WRITE (MP,FMT=9020) (OPLIST(I),I=7,17) WRITE (MP,FMT=9030) (PLIST(I),I=1,8) END IF ELSE WRITE (MP,FMT=9040) ICALLS + 1 END IF ICALLS = ICALLS + 1 c c--- Compute starting indices for work vectors W,IW c MARK = 0 CALL DIVWRK(N,NZ,LW,LIW,MARK) IF (MARK.LT.0) THEN INFORM = -3 RETURN END IF c c--- Outer Newton iteration of the TN method c NSP = 3*N + 4*MAX(N,NZ) CALL OUTER(N,X,F,G,OPLIST,PLIST,INFORM,NZ,NSP,W(IXP),W(IXY), + W(IXD),W(IXHD),W(IXZ),W(IXA),W(IXRSP),IW(IXPER), + IW(IXIPER),IW(IXIA),IW(IXJA),IW(IXISP),CALFGH,CALPAT, + CALHDP,W(IXD1)) RETURN 9000 FORMAT (/,2X,'SETLIS MUST BE CALLED BEFORE TNMIN! ') 9010 FORMAT (/,' PARAMETER INFO, ENTERING TNPACK',/,8X,' N = ',I8,/, + 5X,'OPLIST:',/,8X,' 1) IPCG = ',I8,8X, + '(preconditioning option)',/,8X,' 2) IPRINT = ',I8,8X, + '(printing option)',/,8X,' 3) IPFREQ = ',I8,8X, + '(printing frequency)',/,8X,' 4) MXITN = ',I8,8X, + '(max Newton itns.)',/,8X,' 5) MXITCG = ',I8,8X, + '(max PCG itns.)',/,8X,' 6) MAXNF = ',I8,8X, + '(max F&G evals.)') 9020 FORMAT (8X,' 7) IORDER = ',I8,8X, + '(M-reordering option - calc. per.)',/,8X,' 8) IPERPR = ', + I8,8X,'(M-reordering printing option)',/,8X, + ' 9) IPKNOW = ',I8,8X,'(M-reordering option -', + ' per. known)',/,8X,'10) IHD = ',I8,8X, + '(Numeric HD option)',/,8X,'11) MP = ',I8,8X, + '(printing unit)',/,8X,'12) IEXIT = ',I8,8X, + '(descent option for neg.',' curvature)',/,8X, + '13) ITT = ',I8,8X,'(truncation test option)',/,8X, + '14) MAXFEV = ',I8,8X,'(max F evals. in line search)',/,8X, + '15) TA = ',I8,8X, + '(descent direction test option by PCG)',/,8X, + '16) SRLS = ',I8,8X,'(stopping rule for line search)',/, + 8X,'17) MC = ',I8,8X,'(modified Cholesky)') 9030 FORMAT (5X,'PLIST:',/,8X,' 1) EPSF = ',1P,E10.3,6X, + '(controls conv. test, F accuracy)',/,8X,' 2) EPSG = ', + 1P,E10.3,6X,'(controls conv. test, G accuracy)',/,8X, + ' 3) ETAR = ',1P,E10.3,6X, + '(controls res.-based truncation)',/,8X,' 4) ETAQ = ',1P, + E10.3,6X,'(controls quad.-model based ','truncation)',/,8X, + ' 5) PRODK = ',1P,E10.3,6X,'(controls PCG neg curv. test)' + ,/,8X,' 6) FTOL = ',1P,E10.3,6X, + '(controls F test in line search)',/,8X,' 7) GTOL = ',1P, + E10.3,6X,'(controls G test in line ','search)',/,8X, + ' 8) TAU = ',1P,E10.3,6X,'(parameter in UMC)') 9040 FORMAT (/,6X,'ICALLS =',I8) END C*********************************************************************** C*********************************************************************** SUBROUTINE OUTER(N,X,F,G,OPLIST,PLIST,INFORM,NZ,NSP,P,Y,D,HD,Z,A, + RSP,PER,IPER,IA,JA,ISP,CALFGH,CALPAT,CALHDP,D1) C -------------------------------------------------------- C OUTER: Outer Newton iteration of the TN method C -------------------------------------------------------- C -------------------------------------------------------- C subroutines and functions called: C calfgh, calpat, calhdp (user-supplied) C inner, mlines, odrv, sdrvmd, dnrm2 C -------------------------------------------------------- C -------------------------------------------------------- C Initialize: C c * Set the parameters c * Give an initial guess X_0 c * Compute E(X_0) and G_0 C -------------------------------------------------------- C .. Scalar Arguments .. DOUBLE PRECISION F INTEGER INFORM,N,NSP,NZ C .. C .. Array Arguments .. DOUBLE PRECISION A(N+NZ),D(N),D1(N),G(N),HD(N),P(N),PLIST(20), + RSP(NSP),X(N),Y(N),Z(N) INTEGER IA(N+1),IPER(N),ISP(NSP),JA(N+NZ),OPLIST(20),PER(N) C .. C .. Subroutine Arguments .. EXTERNAL CALFGH,CALHDP,CALPAT C .. C .. Scalars in Common .. DOUBLE PRECISION EPSMCH,ETAQ,ETAR,FTOL,GTOL,PRODK,SQEPS2,SQRTN,TAU INTEGER ICALLS,IEXIT,IHD,ITT,MAXFEV,MC,MP,SRLS,TA C .. C .. Local Scalars .. DOUBLE PRECISION EPSF,EPSF2,EPSF3,EPSG,FOLD,GNORM,LAMBDA,ONE,ONEF, + PERCNT,RNORM,XNORM,XSNORM INTEGER ESP,I,ILINE,IORDER,IPCG,IPERPR,IPFREQ,IPKNOW,IPRINT,ITR, + ITRCG,ITRMAJ,LENA,MAXNF,MODE,MXITCG,MXITN,NFEV,NFUN,NOUT, + OFLAG,OPATH,SFLAG,SPATH LOGICAL T1,T123,T2,T3,T4 C .. C .. External Functions .. DOUBLE PRECISION DNRM2 EXTERNAL DNRM2 C .. C .. External Subroutines .. EXTERNAL INNER,MLINES,ODRV,SDRVMD C .. C .. Intrinsic Functions .. INTRINSIC ABS,DBLE,MOD,SQRT C .. C .. Common blocks .. COMMON /TAINPUT/TA COMMON /TAUINPUT/TAU,SRLS,MC COMMON /TN001/MP,ICALLS COMMON /TN002/PRODK,ETAR,ETAQ,IEXIT,IHD,ITT COMMON /TN003/FTOL,GTOL,MAXFEV COMMON /TN006/EPSMCH,SQEPS2,SQRTN C .. C .. Save statement .. SAVE /TN001/,/TN002/,/TN003/,/TN006/ C .. C .. Local Arrays .. DOUBLE PRECISION RDUM(1),ZDUM(1) C .. IPCG = OPLIST(1) IPRINT = OPLIST(2) IPFREQ = OPLIST(3) MXITN = OPLIST(4) MXITCG = OPLIST(5) MAXNF = OPLIST(6) IORDER = OPLIST(7) IPERPR = OPLIST(8) IPKNOW = OPLIST(9) IHD = OPLIST(10) MP = OPLIST(11) IEXIT = OPLIST(12) ITT = OPLIST(13) MAXFEV = OPLIST(14) EPSF = PLIST(1) EPSG = PLIST(2) ETAR = PLIST(3) ETAQ = PLIST(4) PRODK = PLIST(5) FTOL = PLIST(6) GTOL = PLIST(7) c-- xie: new parameters TA = OPLIST(15) SRLS = OPLIST(16) MC = OPLIST(17) TAU = PLIST(8) c---------------------------- ITRMAJ = 0 ITRCG = 0 NFUN = 0 EPSF2 = SQRT(EPSF) EPSF3 = EPSF** (1./3.) ONE = 1.0D0 C -------------------------------------------------------- C determine the pattern of the preconditioner M C -------------------------------------------------------- IF (IPCG.EQ.1) THEN CALL CALPAT(N,X,A,IA,JA) LENA = IA(N+1) - 1 PERCNT = (DBLE(LENA)/DBLE(N* (N+1)/2))*100.0 IF (LENA.GT. (N+NZ)) THEN INFORM = -4 WRITE (MP,FMT=9000) LENA,N + NZ GO TO 30 END IF IF (ICALLS.LE.1) WRITE (MP,FMT=9010) LENA,PERCNT END IF C -------------------------------------------------------- C compute f, g, H, and M at x0 C -------------------------------------------------------- NOUT = 2 IF (IPCG.EQ.1) NOUT = 3 CALL CALFGH(N,X,F,G,A,IA,JA,NOUT) NFUN = NFUN + 1 IF (IPRINT.GE.1) THEN WRITE (MP,FMT=9020) WRITE (MP,FMT=9030) (X(I),I=1,N) IF (IPRINT.GE.2) THEN WRITE (MP,FMT=9040) WRITE (MP,FMT=9030) (G(I),I=1,N) END IF END IF C -------------------------------------------------------- C check for convergence at the first iteration C -------------------------------------------------------- XNORM = DNRM2(N,X,1)/SQRTN GNORM = DNRM2(N,G,1)/SQRTN IF (GNORM.LT. (EPSG* (ONE+ABS(F)))) THEN WRITE (MP,FMT=9050) ITRMAJ,F,GNORM INFORM = 0 WRITE (MP,FMT=9060) INFORM ITRCG = 0 GO TO 30 END IF C -------------------------------------------------------- C when reordering M, call ODRV to: i) compute the permutation C arrays and then prepare for symmetric reordering of M, or C ii) just prepare for the symmetric reordering of M when the C permutation arrays are known (IPKNOW=1). Then call SDRVMD C to factorize M symbolically. C -------------------------------------------------------- IF (IPCG.EQ.1) THEN IF (IPKNOW.EQ.1) THEN OPATH = 5 ELSE OPATH = 4 DO 10 I = 1,N PER(I) = I IPER(I) = I 10 CONTINUE END IF OFLAG = 0 IF (IORDER.EQ.1) THEN CALL ODRV(N,IA,JA,A,PER,IPER,NSP,ISP,OPATH,OFLAG) IF (IPERPR.EQ.1 .AND. IPKNOW.NE.1) THEN WRITE (MP,FMT=9070) WRITE (MP,FMT=9080) (PER(I),I=1,N) END IF END IF SPATH = 4 IF (OFLAG.EQ.0) CALL SDRVMD(N,PER,IPER,IA,JA,A,RDUM,ZDUM,NSP, + ISP,RSP,ESP,SPATH,SFLAG) IF (OFLAG.NE.0 .OR. SFLAG.NE.0) THEN WRITE (MP,FMT=9090) OFLAG,SFLAG GO TO 30 END IF END IF WRITE (MP,FMT=9050) ITRMAJ,F,GNORM C -------------------------------------------------------- C MAIN LOOP BEGINS C -------------------------------------------------------- 20 CONTINUE MODE = 0 ITRMAJ = ITRMAJ + 1 C -------------------------------------------------------- C print x and f if specified C -------------------------------------------------------- IF (IPFREQ.GT.0) THEN IF (MOD(ITRMAJ,IPFREQ).EQ.0) THEN WRITE (MP,FMT=9100) ITRMAJ WRITE (MP,FMT=9030) (X(I),I=1,N) END IF END IF C -------------------------------------------------------- C check if either the number of maximum-allowed Newton C iterations or function evaluations has been exceeded C -------------------------------------------------------- IF (ITRMAJ.GE.MXITN) THEN INFORM = -5 WRITE (MP,FMT=9110) ITRMAJ,MXITN GO TO 30 END IF IF (NFUN.GT.MAXNF) THEN INFORM = -6 WRITE (MP,FMT=9120) GO TO 30 END IF C -------------------------------------------------------- C call INNER to solve for the search vector p C -------------------------------------------------------- CALL INNER(N,MODE,MXITCG,ITRMAJ,ITR,IPCG,NZ,NSP,X,G,P,Y,PER,IPER, + IA,JA,A,ISP,RSP,D,HD,Z,XNORM,CALFGH,CALHDP,D1) RNORM = DNRM2(N,Y,1) C -------------------------------------------------------- C save old value of f for convergence tests C -------------------------------------------------------- FOLD = F ITRCG = ITRCG + ITR C -------------------------------------------------------- C call line search C -------------------------------------------------------- LAMBDA = ONE ILINE = 0 CALL MLINES(N,X,F,G,P,LAMBDA,ILINE,Y,NFEV,CALFGH) NFUN = NFUN + NFEV C -------------------------------------------------------- C summarize progress and prepare for convergence tests C -------------------------------------------------------- GNORM = DNRM2(N,G,1)/SQRTN XNORM = DNRM2(N,X,1)/SQRTN WRITE (MP,FMT=9130) ITRMAJ,F,GNORM,ITR,MODE,RNORM,LAMBDA,NFUN XSNORM = LAMBDA* (DNRM2(N,P,1)/SQRTN) C -------------------------------------------------------- C call a routine MONIT, as desired, to print further info C -------------------------------------------------------- C C IF (ICALLS .EQ. 1) CALL MONIT(ITRMAJ,ICALLS,F,GNORM) C -------------------------------------------------------- C four convergence tests are performed: C tests 1 & 2 check for convergence of the f and x sequences C tests 3 & 4 check that |g| is suff. small C -------------------------------------------------------- c EPSF2 = 1.E-8 ONEF = ONE + ABS(F) T1 = (FOLD-F) .LT. (EPSF*ONEF) T2 = XSNORM .LT. (EPSF2* (ONE+XNORM)) T3 = GNORM .LT. (EPSF3*ONEF) T4 = GNORM .LT. (EPSG*ONEF) T123 = T1 .AND. (T2 .AND. T3) C ------------------------------------------------------------------- C to check conv. tests for difficult problems, uncomment next 4 lines C C WRITE (MP,*) ' t1:',FOLD-F,EPSF*ONEF,' ',t1 C WRITE (MP,*) ' t2:',XSNORM,EPSF2*(ONE + XNORM),' ',t2 C WRITE (MP,*) ' t3:',GNORM,EPSF3*ONEF,' ',t3 C WRITE (MP,*) ' t4:',GNORM,EPSG*ONEF,' ',t4 C ------------------------------------------------------------------- IF (ILINE.EQ.1) THEN IF (T123 .OR. T4) INFORM = 1 ELSE IF (T3 .OR. T4) THEN INFORM = 2 ELSE INFORM = 3 END IF IF (ILINE.EQ.0) THEN WRITE (MP,FMT=9140) ELSE IF (ILINE.EQ.2) THEN WRITE (MP,FMT=9150) ELSE IF (ILINE.EQ.3) THEN WRITE (MP,FMT=9160) ELSE IF (ILINE.EQ.4) THEN WRITE (MP,FMT=9170) ELSE IF (ILINE.EQ.5) THEN WRITE (MP,FMT=9180) ELSE IF (ILINE.EQ.6) THEN WRITE (MP,FMT=9190) ELSE IF (ILINE.EQ.7) THEN WRITE (MP,FMT=9200) END IF END IF IF (INFORM.EQ.1 .OR. ILINE.NE.1) THEN WRITE (MP,FMT=9210) ITRCG WRITE (MP,FMT=9220) INFORM IF (ILINE.NE.1 .AND. INFORM.NE.1) WRITE (MP,FMT=9230) IF (T1) WRITE (MP,FMT=9240) IF (T2) WRITE (MP,FMT=9250) IF (T3) WRITE (MP,FMT=9260) IF (T4) WRITE (MP,FMT=9270) WRITE (MP,FMT=9280) IF (IPRINT.GE.1) THEN WRITE (MP,FMT=9290) WRITE (MP,FMT=9030) (X(I),I=1,N) IF (IPRINT.GE.2) THEN WRITE (MP,FMT=9300) WRITE (MP,FMT=9030) (G(I),I=1,N) END IF END IF GO TO 30 END IF C -------------------------------------------------------- C line search OK and conv. criteria are not satisfied; C compute H and M, prepare for symmetric reordering of M, C and continue main loop C -------------------------------------------------------- IF (IPCG.EQ.0 .AND. IHD.EQ.1) GO TO 20 IF (IPCG.EQ.0 .AND. IHD.EQ.0) THEN NOUT = 4 ELSE IF (IPCG.EQ.1 .AND. IHD.EQ.0) THEN NOUT = 5 ELSE NOUT = 6 END IF CALL CALFGH(N,X,F,G,A,IA,JA,NOUT) IF (IPCG.EQ.1 .AND. IORDER.EQ.1) THEN OFLAG = 0 OPATH = 5 c---xie: re-calculate IA and JA because they have been changed by c using reordering routine ODRV CALL CALPAT(N,X,A,IA,JA) c--------------------------------------------------- CALL ODRV(N,IA,JA,A,PER,IPER,NSP,ISP,OPATH,OFLAG) IF (OFLAG.NE.0) THEN WRITE (MP,FMT=9310) OFLAG GO TO 30 END IF END IF GO TO 20 C C -------------------------------------------------------- C MAIN LOOP ENDS C -------------------------------------------------------- 30 CONTINUE RETURN 9000 FORMAT (/,1X,' INSUFFICIENT STORAGE FOR M. LENGTH of M = ',I8,/, + 1X,' STORAGE AVAILABLE = N+NZ = ',I8,/) 9010 FORMAT (/,1X,'***** upper-M has ',I8,' nonzeros (',F6.2, + ' % ) ***** ',/) 9020 FORMAT (/,' INITIAL X: ',/) 9030 FORMAT (6 (3X,F10.6)) 9040 FORMAT (/,' INITIAL G: ',/) 9050 FORMAT (2X,'____________________________',/,' ITN.',1X, + ' F ',5X,' |G| ',2X,' #CG ',3X,'mode',2X, + ' |R| ',2X,' LAMBDA',3X,' NFG',/,I5,1P,E13.4,1P,E12.2) 9060 FORMAT (' INFORM = ',I5,/,' (|G(X0)| sufficiently small)',/,1X, + '____________________________',/) 9070 FORMAT (/,1X,' PERMUTATION ARRAY FROM ODRV: ',/) 9080 FORMAT (12 (1X,I5)) 9090 FORMAT (/,1X,' OFLAG = ',I6,' SFLAG = ',I6,/,1X, + 'INSUFF. STORAGE', + ' FOR YSMP (INCREASE NZ) OR DUPLICATE ENTRY IN A',/) 9100 FORMAT (/,' X AT ITRMAJ',I5,':',/) 9110 FORMAT (/,1X,' MXITN (Max Newton iterations) EXCEEDED!',2 (I8),/) 9120 FORMAT (/,1X,' MAXNF (Max function evaluations) EXCEEDED!',/) 9130 FORMAT (I5,1P,E13.4,1P,E12.2,1X,I5,4X,I2,1X,1P,E11.2,1P,E11.2,I5) 9140 FORMAT (/,1X,' IMPROPER INPUT PARAMETERS DURING THE LINE SEARCH.', + /) 9150 FORMAT (/,1X,' RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY', + ' IN THE LINE SEARCH',/,'IS AT MOST XTOL.',/) 9160 FORMAT (/,1X,' NUMBER OF CALLS TO FUNCTION IN THE', + ' LINE SEARCH HAS REACHED MAXFEV.') 9170 FORMAT (/,1X,' THE STEP IN THE LINE SEARCH IS TOO SMALL.',/) 9180 FORMAT (/,1X,' THE STEP IN THE LINE SEARCH IS TOO LARGE.',/) 9190 FORMAT (/,1X,' ROUNDING ERRORS PREVENT FURTHER PROGRESS IN ',/, + ' THE LINE SEARCH.',/) 9200 FORMAT (/,1X,' P IS NOT A DESCENT DIRECTION.',/) 9210 FORMAT (/,1X,'____________________________',2X,'#CG,tot',/,29X,I7) 9220 FORMAT (' INFORM = ',I5) 9230 FORMAT (/,' CHECK LS PARAMETERS, TRY RELAXING CONV. CRITERIA.',/, + ' AND TRY A DIFFERENT TRUNCATION PARAMETER).',/) 9240 FORMAT (' CONV. TEST 1 SATISFIED ') 9250 FORMAT (' CONV. TEST 2 SATISFIED ') 9260 FORMAT (' CONV. TEST 3 SATISFIED ') 9270 FORMAT (' CONV. TEST 4 SATISFIED ') 9280 FORMAT (1X,'____________________________') 9290 FORMAT (/,' FINAL X: ',/) 9300 FORMAT (/,' FINAL G: ',/) 9310 FORMAT (/,1X,' IN TNMIN, after ODRV for new M, OFLAG = ',I6,/,1X, + 'INSUFFICIENT STORAGE FOR YSMP. INCREASE NZ.',/) END C*********************************************************************** C*********************************************************************** SUBROUTINE INNER(N,MODE,MXITCG,ITRMAJ,ITR,IPCG,NZ,NSP,X,G,P,RES, + PER,IPER,IA,JA,A,ISP,RSP,D,HD,Z,XNORM,CALFGH, + CALHDP,D1) C -------------------------------------------------------- C INNER: Inner PCG Iteration C -------------------------------------------------------- C -------------------------------------------------------- C subroutines and functions called: C calhdp (user-supplied), ourhdp, C sdrvmd, dcopy, daxpy, dnrm2, ddot C -------------------------------------------------------- C -------------------------------------------------------- C STEP 1 - Initialization C C (a) Compute ||g||, set p=0, and res= -g C (b) If preconditioning, perform the numerical C factorization of M and solve Mz=res; C else (M=I) set z=res C (c) Set d=z and compute delta=d*d and z*res C -------------------------------------------------------- C .. Scalar Arguments .. DOUBLE PRECISION XNORM INTEGER IPCG,ITR,ITRMAJ,MODE,MXITCG,N,NSP,NZ C .. C .. Array Arguments .. DOUBLE PRECISION A(N+NZ),D(N),D1(N),G(N),HD(N),P(N),RES(N), + RSP(NSP),X(N),Z(N) INTEGER IA(N+1),IPER(N),ISP(NSP),JA(N+NZ),PER(N) C .. C .. Subroutine Arguments .. EXTERNAL CALFGH,CALHDP C .. C .. Scalars in Common .. DOUBLE PRECISION ETAQ,ETAR,PRODK,TAU INTEGER ICALLS,IEXIT,IHD,ITT,MC,MP,SRLS,TA C .. C .. Local Scalars .. DOUBLE PRECISION ALPHA,BETA,DELTA,ETA,GNORM,PRODCT,PTG,PTR,QNEW, + QOLD,RITR,RMAJ,RNORM,RTZ,RTZOLD,TEMP,TEMP_OLD, + ZETA INTEGER ESP,I,SFLAG,SPATH C .. C .. External Functions .. DOUBLE PRECISION DDOT,DNRM2 EXTERNAL DDOT,DNRM2 C .. C .. External Subroutines .. EXTERNAL DAXPY,DCOPY,OURHDP,SDRVMD C .. C .. Intrinsic Functions .. INTRINSIC ABS,MIN C .. C .. Common blocks .. COMMON /TAINPUT/TA COMMON /TAUINPUT/TAU,SRLS,MC COMMON /TN001/MP,ICALLS COMMON /TN002/PRODK,ETAR,ETAQ,IEXIT,IHD,ITT C .. C .. Save statement .. SAVE /TN001/,/TN002/ C .. TEMP_OLD = 0.D0 ITR = 0 GNORM = DNRM2(N,G,1) DO I = 1,N P(I) = 0.0D0 RES(I) = -G(I) END DO QOLD = 0.D0 IF (IPCG.EQ.1) THEN SPATH = 2 CALL SDRVMD(N,PER,IPER,IA,JA,A,RES,Z,NSP,ISP,RSP,ESP,SPATH, + SFLAG) IF (SFLAG.GT.0) THEN WRITE (MP,FMT=9000) SFLAG GO TO 30 END IF ELSE CALL DCOPY(N,RES,1,Z,1) END IF CALL DCOPY(N,Z,1,D,1) DELTA = DDOT(N,D,1,D,1) RTZ = DDOT(N,RES,1,Z,1) RNORM = DDOT(N,RES,1,RES,1) RMAJ = ITRMAJ C -------------------------------------------------------- C STEP 2 - MAIN LOOP BEGINS: C (a) compute the Hessian-vector product Hd C (b) if d*Hd is small, exit with a descent direction C chosen according to: On ITR 1, p=-g C on ITR>1, p= current p c Note: we use Test 1A' or Test 2A in (b). C -------------------------------------------------------- 10 CONTINUE c--- Test if ITR is larger than the maximum number of PCG iterations IF (ITR.GE.MXITCG) THEN MODE = 3 GO TO 30 END IF C--- Compute a Hessian/vector product Hd C IHD=1: Compute Hd by a finite-difference of gradient C 0: Compute Hd directly C C D - input vector C HD - resultant product vector C X - current coordinate C G - gradient C-------------------------------------- IF (IHD.EQ.1) THEN CALL OURHDP(N,D,HD,X,G,RSP(1),ITR,XNORM,CALFGH) ELSE CALL CALHDP(N,D,HD,X,G) END IF ITR = ITR + 1 RITR = ITR PRODCT = DDOT(N,D,1,HD,1) IF (MC.EQ.1) GO TO 20 C--- xie: The singularity test ZETA = 1.D-20 IF (ABS(RTZ).LE.ZETA*RNORM*RNORM .OR. ABS(PRODCT).LE.ZETA) THEN IF (ITR.EQ.1) THEN MODE = 1 DO I = 1,N P(I) = -G(I) END DO ELSE MODE = 5 END IF GO TO 30 END IF C--------------------------------------------- 20 CONTINUE c--- xie: Determine descent search directions c c TA = 1: use Test 1A' (the modified negative curvature test). c A default choice c TA = 2: use Test 2A (the strong negative curvature test) c---------------------------------------------------------- c---xie: Test 1A': IF (TA.EQ.1) THEN IF (PRODCT.LE. (DELTA*PRODK)) THEN IF (ITR.EQ.1) THEN MODE = 1 DO I = 1,N P(I) = -G(I) END DO ELSE MODE = 2 c--- xie: 8/20/98 To test D as search direction IF (IEXIT.EQ.2) THEN CALL DCOPY(N,D,1,P,1) END IF c---------------------------------------------- END IF GO TO 30 END IF END IF c--------------------- end of Test 1A' ALPHA = RTZ/PRODCT RTZOLD = RTZ c-- xie: keep the current direction P in D1 DO I = 1,N D1(I) = P(I) END DO c--- Update P: P = P + ALPHA * D CALL DAXPY(N,ALPHA,D,1,P,1) c--- xie: Test 2A c c MODE = 1: G*P > 0 at the first CG or PCG iteration c MODE = 4: G*P > 0 or G*P becomes increasing c c PRODK = 10E-10 as the default value c--------------------------------------------------------- IF (TA.EQ.2) THEN TEMP = DDOT(N,P,1,G,1) IF (TEMP.GE.TEMP_OLD) THEN IF (ITR.EQ.1) THEN MODE = 1 DO I = 1,N P(I) = -G(I) END DO ELSE MODE = 4 c--- exit with the previous PCG iterate (i.e., D1) as the C search direction DO I = 1,N P(I) = D1(I) END DO END IF GO TO 30 ELSE TEMP_OLD = TEMP - ZETA END IF END IF c-------------------------- End of Test 2A c c--- Update residual vector: RES = RES - ALPHA * HD c CALL DAXPY(N,-ALPHA,HD,1,RES,1) c c--- STEP 3 - Truncation Test c IF (ITT.EQ.1) THEN PTR = DDOT(N,P,1,RES,1) PTG = DDOT(N,P,1,G,1) QNEW = 0.5D0* (PTR+PTG) IF (RITR* (1.D0- (QOLD/QNEW)).LE.ETAQ) THEN MODE = 0 GO TO 30 END IF QOLD = QNEW ELSE RNORM = DNRM2(N,RES,1) ETA = MIN(ETAR/RMAJ,GNORM) IF (RNORM.LE. (GNORM*ETA)) THEN MODE = 0 GO TO 30 END IF END IF c---------------------- End the test C -------------------------------------------------------- C STEP 4 - Continuation of PCG. C (a) if preconditioning, solve Mz=res by re-using C the factors of M; else (M=I) set z=res, C (b) update d and delta C (c) go on to the next PCG iteration. C -------------------------------------------------------- IF (IPCG.EQ.1) THEN SPATH = 3 c c--- Solving Mz=r by the unusual modified Cholesky factorization c CALL SDRVMD(N,PER,IPER,IA,JA,A,RES,Z,NSP,ISP,RSP,ESP,SPATH, + SFLAG) IF (SFLAG.GT.0) THEN WRITE (MP,FMT=9000) SFLAG GO TO 30 END IF ELSE c c--- COPIES A VECTOR, RES, TO A VECTOR, Z c CALL DCOPY(N,RES,1,Z,1) END IF RTZ = DDOT(N,RES,1,Z,1) BETA = RTZ/RTZOLD DO I = 1,N D(I) = BETA*D(I) + Z(I) END DO DELTA = DDOT(N,D,1,D,1) GO TO 10 C -------------------------------------------------------- C MAIN LOOP ENDS C -------------------------------------------------------- 30 CONTINUE RETURN 9000 FORMAT (/,1X,' IN TNPCG, SFLAG = ',I6,/,1X, +'INSUFF. STORAGE FOR YSMP (INCREASE NZ) OR DUPLICATE ENTR +Y IN A',/) END C*********************************************************************** C*********************************************************************** SUBROUTINE OURHDP(N,D,HD,X,G,Y,ITR,XNORM,CALFGH) C ------------------------------------------------------ C OURHDP: Our numeric Hessian/vector multiplication C ------------------------------------------------------ C -------------------------------------------------------- C DELH is set to: 2*SQRT(EPS)*(1+|X|) / |D|, with division C precautions and upper/lower limits C -------------------------------------------------------- C .. Parameters .. DOUBLE PRECISION ONE,DMAX,DIMAX PARAMETER (ONE=1.D0,DMAX=0.1D0,DIMAX=ONE/DMAX) C .. C .. Scalar Arguments .. DOUBLE PRECISION XNORM INTEGER ITR,N C .. C .. Array Arguments .. DOUBLE PRECISION D(N),G(N),HD(N),X(N),Y(N) C .. C .. Subroutine Arguments .. EXTERNAL CALFGH C .. C .. Scalars in Common .. DOUBLE PRECISION EPSMCH,SQEPS2,SQRTN C .. C .. Local Scalars .. DOUBLE PRECISION DDEN,DELH,DELHI,DNORM,DNUM,FNEW INTEGER I C .. C .. External Functions .. DOUBLE PRECISION DNRM2 EXTERNAL DNRM2 C .. C .. Intrinsic Functions .. INTRINSIC MAX C .. C .. Common blocks .. COMMON /TN006/EPSMCH,SQEPS2,SQRTN C .. C .. Save statement .. SAVE /TN006/,DNUM C .. C .. Local Arrays .. DOUBLE PRECISION ADUM(1) INTEGER IADUM(1),JADUM(1) C .. IF (ITR.EQ.0) DNUM = SQEPS2* (ONE+XNORM*SQRTN) DNORM = DNRM2(N,D,1) DDEN = MAX(DNUM*DIMAX,DNORM) DELH = MAX(DNUM/DDEN,DNUM*0.01) DELHI = ONE/DELH DO 10 I = 1,N Y(I) = X(I) + DELH*D(I) 10 CONTINUE CALL CALFGH(N,Y,FNEW,HD,ADUM,IADUM,JADUM,1) DO 20 I = 1,N HD(I) = (HD(I)-G(I))*DELHI 20 CONTINUE RETURN END C*********************************************************************** C*********************************************************************** SUBROUTINE SETLIS(N,OPLIST,PLIST,INFORM) C -------------------------------------------------------- C SETLIS: Set sample values for the TNMIN-call list of C options and parameters (OPLIST & PLIST). C The user should then change some of these items C - to suit the problem - before the TNMIN call C -------------------------------------------------------- C -------------------------------------------------------- C function called: dmach (blas) C -------------------------------------------------------- C .. Scalar Arguments .. INTEGER INFORM,N C .. C .. Array Arguments .. DOUBLE PRECISION PLIST(20) INTEGER OPLIST(20) C .. C .. Scalars in Common .. DOUBLE PRECISION EPSMCH,SQEPS2,SQRTN INTEGER ICALLS,MP C .. C .. External Functions .. DOUBLE PRECISION DMACH EXTERNAL DMACH C .. C .. Intrinsic Functions .. INTRINSIC DBLE,SQRT C .. C .. Common blocks .. COMMON /TN001/MP,ICALLS COMMON /TN006/EPSMCH,SQEPS2,SQRTN C .. C .. Save statement .. SAVE /TN001/,/TN006/ C .. ICALLS = 0 C -------------------------------------------------------- C Brief Description of parameters (see manual for details; C M below refers to the preconditioner): C -------------------------------------------------------- C INFORM - diagnostic/status parameter C ICALLS - counter for number of TNPACK calls (multiple C calls may be made for a series of minimizations) C OPLIST(1) - preconditioning option (0 - No, 1 - Yes) C OPLIST(2) - general output controler (0,1,2) C OPLIST(3) - output controler, frequency of printing X (0,1,...) C OPLIST(4) - limit of total Newton iterations C OPLIST(5) - limit of PCG iterations in each Newton iteration C OPLIST(6) - limit of total function evaluations C OPLIST(7) - M-reordering option (0 - Do not reorder, 1 - Reorder) C OPLIST(8) - printing option of M's reordering (permutation array) C (0 - Do not print, 1 - Print) C OPLIST(9) - option for specifying whether the permutation array C for reordering M is known (0 - Not known, 1 - Known). C (This option is useful for multiple TNPACK minimiza- C tions, involving reorderings of M, if the sparsity C structure of M does not change. The user can reset C OPLIST(9) from 0 to 1 before the second TNMIN call). C OPLIST(10) - Hessian/vector multiplication option (0 - use a user- C supplied routine, 1 - use our default finite-difference C design routine); IHD in COMMON/TN002 C OPLIST(11) - unit number for printing; MP in COMMON/TN001 c c--- xie: modify OPLIST(12) C OPLIST(12) - option selector for combination of exit directions in C case of detection of negative curvature for PCG itns. c Only 1: -G/P is used. The other selectors (such as c 0: -G/-G, 2: -G/D, 3: -Y/-G, 4: -Y/P, C 5: -Y/D, where Y=-[M**(-1)]G) are deleted. c----------------------------------------------------------------- c C OPLIST(13) - option selector for PCG truncation test C (0: residual-based, 1: quadratic-model based); C ITT in COMMON/TN002 C OPLIST(14) - limit for number of function calls in the line search C (at each Newton iteration); MAXFEV in COMMON/TN003 c C PLIST(1) - desired F accuracy in min. conv. test C PLIST(2) - desired ||G|| accuracy in min. conv. test C PLIST(3) - truncation parameter for residual test; C ETAR in COMMON/TN002 C PLIST(4) - truncation parameter for quadratic-model test; C ETAQ in COMMON/TN002 C PLIST(5) - tolerance for 'negative curvature' test; C PRODK in COMMON/TN002 C PLIST(6) - line-search conv. parameter, controlling function C decrease; FTOL in COMMON/TN003 C PLIST(7) - line-search conv. parameter, controlling reduction C of derivative magnitude; GTOL in COMMON/TN003 c c--- xie: new parameters c OPLIST(15) - option selector for termination rule of PCG C 1: the modified negative curvature test 1A' C 2: the descent direction test 2A) C Default value is 1 C C OPLIST(16) - option selector for line search stopping rule C 1: Criterion 1 C 2: Criterion 2 (more lenient stopping rule) C Default value is 1 C OPLIST(17) - option selector for modified Cholesky factorization C 1: the MC by Gill and Murray C 2: the UMC by Xie and Schlick C Default value is 1 c PLIST(8) - unusual modified Cholesky factorization parameter TAU c Default value is 10.d0 C --------------------------------------------------------------- IF (N.LE.0) THEN WRITE (*,FMT=9000) N RETURN END IF INFORM = -10 OPLIST(1) = 0 OPLIST(2) = 0 OPLIST(3) = 0 OPLIST(4) = 100*N OPLIST(5) = 10*N OPLIST(6) = 100*N OPLIST(7) = 0 OPLIST(8) = 0 OPLIST(9) = 0 OPLIST(10) = 1 OPLIST(11) = 6 OPLIST(12) = 1 OPLIST(13) = 0 OPLIST(14) = 20 MP = OPLIST(11) EPSMCH = DMACH(1) SQEPS2 = 2.D0*SQRT(EPSMCH) SQRTN = SQRT(DBLE(N)) PLIST(1) = 1.0D-10 PLIST(2) = SQRT(EPSMCH) PLIST(3) = 0.5D0 PLIST(4) = 0.5D0 PLIST(5) = 1.0D-10 PLIST(6) = 1.0D-04 PLIST(7) = 0.9D0 c-- xie: defualt new parameters PLIST(8) = 1.0D1 OPLIST(15) = 1 OPLIST(16) = 1 OPLIST(17) = 1 c-------------------------------- RETURN 9000 FORMAT (/,2X,' N <= 0, N = ',I10) END C*********************************************************************** C*********************************************************************** SUBROUTINE CHKLIS(N,OPLIST,PLIST,MARK) C -------------------------------------------------------- C CHKLIS: Check input options & parameters (OPLIST,PLIST) C -------------------------------------------------------- C .. Scalar Arguments .. INTEGER MARK,N C .. C .. Array Arguments .. DOUBLE PRECISION PLIST(20) INTEGER OPLIST(20) C .. C .. Scalars in Common .. DOUBLE PRECISION EPSMCH,SQEPS2,SQRTN C .. C .. Local Scalars .. DOUBLE PRECISION ONE,SQEPS,ZERO C .. C .. Intrinsic Functions .. INTRINSIC SQRT C .. C .. Common blocks .. COMMON /TN006/EPSMCH,SQEPS2,SQRTN C .. C .. Save statement .. SAVE /TN006/ C .. ZERO = 0.D0 ONE = 1.D0 SQEPS = SQRT(EPSMCH) IF (N.LE.0) THEN WRITE (*,FMT=9000) N MARK = -1 RETURN END IF IF (OPLIST(1).NE.0 .AND. OPLIST(1).NE.1) THEN OPLIST(1) = 0 WRITE (*,FMT=9010) END IF IF (OPLIST(2).LT.0 .OR. OPLIST(2).GT.2) THEN OPLIST(2) = 0 WRITE (*,FMT=9020) END IF IF (OPLIST(3).LT.0) THEN OPLIST(3) = 0 WRITE (*,FMT=9030) END IF IF (OPLIST(4).LE.0 .OR. OPLIST(5).LE.0 .OR. OPLIST(6).LE.0) THEN WRITE (*,FMT=9040) MARK = -4 RETURN END IF IF (OPLIST(1).EQ.0) THEN OPLIST(7) = 0 OPLIST(8) = 0 OPLIST(9) = 0 ELSE IF (OPLIST(7).NE.0 .AND. OPLIST(7).NE.1) THEN OPLIST(7) = 0 WRITE (*,FMT=9050) END IF IF (OPLIST(7).EQ.0 .AND. (OPLIST(8).NE.0.OR. + OPLIST(9).NE.0)) THEN OPLIST(8) = 0 OPLIST(9) = 0 END IF END IF IF (OPLIST(10).NE.0 .AND. OPLIST(10).NE.1) THEN OPLIST(10) = 1 WRITE (*,FMT=9060) END IF IF (OPLIST(11).LE.0) THEN WRITE (*,FMT=9070) MARK = -11 RETURN END IF IF (OPLIST(12).LT.0 .OR. OPLIST(12).GT.5) THEN OPLIST(12) = 1 WRITE (*,FMT=9080) OPLIST(12) END IF IF (OPLIST(13).NE.0 .AND. OPLIST(13).NE.1) THEN OPLIST(13) = 0 WRITE (*,FMT=9090) END IF IF (OPLIST(14).LE.0 .OR. OPLIST(14).GT.40) THEN OPLIST(14) = 20 WRITE (*,FMT=9100) END IF IF (PLIST(1).LT.EPSMCH) THEN PLIST(1) = EPSMCH*100.D0 WRITE (*,FMT=9110) PLIST(1) END IF IF (PLIST(2).LT.EPSMCH) THEN PLIST(2) = SQEPS*10.D0 WRITE (*,FMT=9120) PLIST(2),SQEPS END IF IF (OPLIST(13).EQ.0) THEN IF (PLIST(3).LT.ZERO .OR. PLIST(3).GT.ONE) THEN PLIST(3) = 0.5D0 WRITE (*,FMT=9130) PLIST(3) ELSE IF (PLIST(3).LT.1.D-4) THEN WRITE (*,FMT=9150) END IF END IF IF (OPLIST(13).EQ.1) THEN IF (PLIST(4).LT.ZERO .OR. PLIST(4).GT.ONE) THEN PLIST(4) = 0.5D0 WRITE (*,FMT=9140) PLIST(4) ELSE IF (PLIST(4).LT.1.D-4) THEN WRITE (*,FMT=9160) END IF END IF IF (PLIST(5).LT.EPSMCH) THEN PLIST(5) = SQEPS*0.1D0 WRITE (*,FMT=9170) PLIST(5) END IF IF (PLIST(6).LT.ZERO .OR. PLIST(7).LT.ZERO .OR. + PLIST(6).GT.ONE .OR. PLIST(7).GT.ONE .OR. + PLIST(6).GE.PLIST(7)) THEN PLIST(6) = 1.0D-04 PLIST(7) = 0.9D0 WRITE (*,FMT=9180) END IF c-- xie: set default values of OPLIST(15), OPLIST(16), OPLIST(17), c-- and PLIST(8) IF (OPLIST(15).LE.0 .OR. OPLIST(15).GT.2) OPLIST(15) = 1 IF (OPLIST(16).LE.0 .OR. OPLIST(16).GT.2) OPLIST(16) = 1 IF (OPLIST(17).LE.0 .OR. OPLIST(17).GT.2) OPLIST(17) = 1 IF (PLIST(8).LT.ZERO) PLIST(8) = 1.0D1 c------------------------------------------------ RETURN 9000 FORMAT (/,2X,'N <= 0, N = ',I10) 9010 FORMAT (/,2X,'OPLIST(1) OUT-OF-RANGE, reset to 0 (no PCG)') 9020 FORMAT (/,2X,'OPLIST(2) OUT-OF-RANGE, reset to 0') 9030 FORMAT (/,2X,'OPLIST(3) OUT-OF-RANGE, reset to 0') 9040 FORMAT (/,2X,'OPLIST (4,5, and/or 6) <= 0') 9050 FORMAT (/,2X,'OPLIST(7) OUT-OF-RANGE, reset to 0 (no reordering)' + ) 9060 FORMAT (/,2X, + 'OPLIST(10) OUT-OF-RANGE, reset to 1 (our Hd routine)') 9070 FORMAT (/,2X,'OPLIST(11), PRINTING UNIT, <= 0') 9080 FORMAT (/,2X,'OPLIST(12) OUT-OF-RANGE, reset to',I10) 9090 FORMAT (/,2X, + 'OPLIST(13) OUT-OF-RANGE, reset to 0 (residual test)') 9100 FORMAT (/,2X,'OPLIST(14) <=0 or unreasonably high, reset to 20') 9110 FORMAT (/,2X,'PLIST(1) < EPSMCH, reset to',1P,E15.5) 9120 FORMAT (/,2X,'PLIST(2) < EPSMCH, reset to',1P,E15.5,/,2X, + 'However, around SQRT(EPSMCH) is recommended:',1P,E15.5) 9130 FORMAT (/,2X,'PLIST(3) OUT-OF-RANGE, reset to',1P,E15.5) 9140 FORMAT (/,2X,'PLIST(4) OUT-OF-RANGE, reset to',1P,E15.5) 9150 FORMAT (/,2X,'NOTE: A SMALL INPUT VALUE FOR PLIST(3) MAY NOT',/, + 2X,'EXPLOIT THE BENEFIT OF TRUNCATION') 9160 FORMAT (/,2X,'NOTE: A SMALL INPUT VALUE FOR PLIST(4) MAY NOT',/, + 2X,'EXPLOIT THE BENEFIT OF TRUNCATION') 9170 FORMAT (/,2X,'PLIST(5) < EPSMCH, reset to',1P,E15.5) 9180 FORMAT (/,2X,'PLIST(6 and/or 7) OUT-OF-RANGE, both were reset',/, + 2X,'to default values') END C*********************************************************************** C*********************************************************************** SUBROUTINE DIVWRK(N,NZ,LW,LIW,MARK) C -------------------------------------------------------- C DIVWRK: Compute starting indices for work vectors W,IW C -------------------------------------------------------- C -------------------------------------------------------- C W is partitioned to vectors: C P(n),Y(n),D(n),HD(n),Z(n),A(n+nz),RSP(nsp), C IW is partitioned to vectors: C PER(n),IPER(n),IA(n+1),JA(n+nz),ISP(nsp), C where nz = # of nonzeros in strict upper-triang. of M, C and nsp = 3*n + 4*max(n,nz). C -------------------------------------------------------- C .. Scalar Arguments .. INTEGER LIW,LW,MARK,N,NZ C .. C .. Scalars in Common .. INTEGER ICALLS,IXA,IXD,IXD1,IXHD,IXIA,IXIPER,IXISP,IXJA,IXP,IXPER, + IXRSP,IXY,IXZ,MP C .. C .. Local Scalars .. INTEGER NEEDIW,NEEDW,NSP C .. C .. Intrinsic Functions .. INTRINSIC MAX C .. C .. Common blocks .. COMMON /TN001/MP,ICALLS COMMON /TN004/IXP,IXY,IXD,IXHD,IXZ,IXA,IXRSP,IXD1 COMMON /TN005/IXPER,IXIPER,IXIA,IXJA,IXISP C .. C .. Save statement .. SAVE /TN001/,/TN004/,/TN005/ C .. NSP = 3*N + 4*MAX(N,NZ) NEEDW = 6*N + NZ + NSP NEEDIW = 4*N + NZ + NSP + 1 IF (N.LE.0 .OR. NZ.LT.0 .OR. LW.LT.NEEDW .OR. LIW.LT.NEEDIW) THEN WRITE (*,FMT=9000) N,NZ,LW,NEEDW,LIW,NEEDIW MARK = -1 RETURN END IF IXP = 1 IXY = IXP + N IXD = IXY + N IXHD = IXD + N IXZ = IXHD + N IXA = IXZ + N IXRSP = IXA + N + NZ IXD1 = IXRSP + 3*N + 4*NZ IXPER = 1 IXIPER = 1 + N IXIA = IXIPER + N IXJA = IXIA + N + 1 IXISP = IXJA + N + NZ RETURN 9000 FORMAT (/,2X,'ERROR in DIMENSIONS of N,NZ,LW,or LIW',/,2X,'N = ', + I10,'NZ = ',I10,/,2X,'LW = ',I10,'(NEED ',I10,')','LIW = ', + I10,'(NEED ',I10,')') END C*********************************************************************** C*********************************************************************** SUBROUTINE MLINES(N,X,F,G,S,STP,INFO,WA,NFEV,CALFGH) C -------------------------------------------------------- C MLINES: Line-search algorithm of More' and Thuente C -------------------------------------------------------- C -------------------------------------------------------- C subroutines and functions called: C calfgh (user-supplied) C newstp C daxpy, dcopy, ddot (blas) C -------------------------------------------------------- C .. Scalar Arguments .. DOUBLE PRECISION F,STP INTEGER INFO,N,NFEV C .. C .. Array Arguments .. DOUBLE PRECISION G(N),S(N),WA(N),X(N) C .. C .. Subroutine Arguments .. EXTERNAL CALFGH C .. C .. Scalars in Common .. DOUBLE PRECISION FTOL,GTOL,TAU INTEGER MAXFEV,MC,SRLS C .. C .. Local Scalars .. DOUBLE PRECISION DG,DGINIT,DGM,DGTEST,DGX,DGXM,DGY,DGYM,FINIT,FM, + FTEST1,FX,FXM,FY,FYM,P5,P66,STMAX,STMIN,STPMAX, + STPMIN,STX,STY,WIDTH,WIDTH1,XTOL,XTRAPF,ZERO INTEGER INFOC LOGICAL BRACKT,STAGE1 C .. C .. External Functions .. DOUBLE PRECISION DDOT EXTERNAL DDOT C .. C .. External Subroutines .. EXTERNAL DAXPY,DCOPY,NEWSTP C .. C .. Intrinsic Functions .. INTRINSIC ABS,MAX,MIN C .. C .. Common blocks .. COMMON /TAUINPUT/TAU,SRLS,MC COMMON /TN003/FTOL,GTOL,MAXFEV C .. C .. Save statement .. SAVE C .. C .. Local Arrays .. DOUBLE PRECISION ADUM(1) INTEGER IDUM(1),JDUM(1) C .. C .. Data statements .. DATA P5,P66,XTRAPF,ZERO/0.5D0,0.66D0,4.0D0,0.0D0/ DATA XTOL,STPMIN,STPMAX/1.0D-17,1.0D-20,1.0D+20/ C .. INFOC = 1 C -------------------------------------------------------- C CHECK THE INPUT PARAMETERS FOR ERRORS. C -------------------------------------------------------- IF (N.LE.0 .OR. STP.LE.ZERO .OR. FTOL.LT.ZERO .OR. + GTOL.LT.ZERO .OR. XTOL.LT.ZERO .OR. STPMIN.LT.ZERO .OR. + STPMAX.LT.STPMIN .OR. MAXFEV.LE.0) RETURN C -------------------------------------------------------- C COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION C AND CHECK THAT S IS A DESCENT DIRECTION. C -------------------------------------------------------- DGINIT = DDOT(N,G,1,S,1) IF (DGINIT.GE.ZERO) THEN WRITE (*,FMT=*) ' GTP in MLINES = ',DGINIT INFO = 7 RETURN END IF C -------------------------------------------------------- C INITIALIZE LOCAL VARIABLES. C -------------------------------------------------------- BRACKT = .FALSE. STAGE1 = .TRUE. NFEV = 0 FINIT = F DGTEST = FTOL*DGINIT WIDTH = STPMAX - STPMIN WIDTH1 = WIDTH/P5 CALL DCOPY(N,X,1,WA,1) C -------------------------------------------------------- C THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP, C FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP. C THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP, C FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF C THE INTERVAL OF UNCERTAINTY. C THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP, C FUNCTION, AND DERIVATIVE AT THE CURRENT STEP. C -------------------------------------------------------- STX = ZERO FX = FINIT DGX = DGINIT STY = ZERO FY = FINIT DGY = DGINIT C -------------------------------------------------------- C START OF ITERATION. C -------------------------------------------------------- 10 CONTINUE C -------------------------------------------------------- C SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND C TO THE PRESENT INTERVAL OF UNCERTAINTY. C -------------------------------------------------------- IF (BRACKT) THEN STMIN = MIN(STX,STY) STMAX = MAX(STX,STY) ELSE STMIN = STX STMAX = STP + XTRAPF* (STP-STX) END IF C -------------------------------------------------------- C FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN. C -------------------------------------------------------- STP = MAX(STP,STPMIN) STP = MIN(STP,STPMAX) C -------------------------------------------------------- C IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET C STP BE THE LOWEST POINT OBTAINED SO FAR. C -------------------------------------------------------- IF ((BRACKT.AND. (STP.LE.STMIN.OR.STP.GE.STMAX)) .OR. + NFEV.GE.MAXFEV-1 .OR. INFOC.EQ.0 .OR. + (BRACKT.AND.STMAX-STMIN.LE.XTOL*STMAX)) STP = STX C -------------------------------------------------------- C EVALUATE THE FUNCTION AND GRADIENT AT STP C AND COMPUTE THE DIRECTIONAL DERIVATIVE. C -------------------------------------------------------- CALL DCOPY(N,WA,1,X,1) CALL DAXPY(N,STP,S,1,X,1) C CALL CALFGH(N,X,F,G,ADUM,IDUM,JDUM,1) C INFO = 0 NFEV = NFEV + 1 DG = DDOT(N,G,1,S,1) FTEST1 = FINIT + STP*DGTEST C -------------------------------------------------------- C TEST FOR CONVERGENCE. C -------------------------------------------------------- IF ((BRACKT.AND. (STP.LE.STMIN.OR.STP.GE.STMAX)) .OR. + INFOC.EQ.0) INFO = 6 IF (STP.EQ.STPMAX .AND. F.LE.FTEST1 .AND. DG.LE.DGTEST) INFO = 5 IF (STP.EQ.STPMIN .AND. (F.GT.FTEST1.OR.DG.GE.DGTEST)) INFO = 4 IF (NFEV.GE.MAXFEV) INFO = 3 IF (BRACKT .AND. STMAX-STMIN.LE.XTOL*STMAX) INFO = 2 c--- xie: the stopping rule of the line search c The original one: IF (SRLS.EQ.1) THEN IF (F.LE.FTEST1 .AND. ABS(DG).LE.GTOL* (-DGINIT)) INFO = 1 ELSE c-- xie: the lenient stopping rule of the line search IF (F.LE.FTEST1 .AND. (DG.GE.GTOL*DGINIT.OR. + DG.LE. (2.0D0-GTOL)*DGINIT)) INFO = 1 END IF C -------------------------------------------------------- C CHECK FOR TERMINATION. C -------------------------------------------------------- IF (INFO.NE.0) RETURN C -------------------------------------------------------- C IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED C FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE. C -------------------------------------------------------- IF (STAGE1 .AND. F.LE.FTEST1 .AND. + DG.GE.MIN(FTOL,GTOL)*DGINIT) STAGE1 = .FALSE. C -------------------------------------------------------- C A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF C WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED C FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE C DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN C OBTAINED BUT THE DECREASE IS NOT SUFFICIENT. C -------------------------------------------------------- IF (STAGE1 .AND. F.LE.FX .AND. F.GT.FTEST1) THEN C -------------------------------------------------------- C DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES. C -------------------------------------------------------- FM = F - STP*DGTEST FXM = FX - STX*DGTEST FYM = FY - STY*DGTEST DGM = DG - DGTEST DGXM = DGX - DGTEST DGYM = DGY - DGTEST C -------------------------------------------------------- C CALL NEWSTP TO UPDATE THE INTERVAL OF UNCERTAINTY C AND TO COMPUTE THE NEW STEP. C -------------------------------------------------------- CALL NEWSTP(STX,FXM,DGXM,STY,FYM,DGYM,STP,FM,DGM,BRACKT,STMIN, + STMAX,INFOC) C -------------------------------------------------------- C RESET THE FUNCTION AND GRADIENT VALUES FOR F. C -------------------------------------------------------- FX = FXM + STX*DGTEST FY = FYM + STY*DGTEST DGX = DGXM + DGTEST DGY = DGYM + DGTEST ELSE C -------------------------------------------------------- C CALL NEWSTP TO UPDATE THE INTERVAL OF UNCERTAINTY C AND TO COMPUTE THE NEW STEP. C -------------------------------------------------------- CALL NEWSTP(STX,FX,DGX,STY,FY,DGY,STP,F,DG,BRACKT,STMIN,STMAX, + INFOC) END IF C -------------------------------------------------------- C FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE C INTERVAL OF UNCERTAINTY. C -------------------------------------------------------- IF (BRACKT) THEN IF (ABS(STY-STX).GE.P66*WIDTH1) STP = STX + P5* (STY-STX) WIDTH1 = WIDTH WIDTH = ABS(STY-STX) END IF C -------------------------------------------------------- C END OF ITERATION. C -------------------------------------------------------- GO TO 10 END C C*********************************************************************** C*********************************************************************** SUBROUTINE NEWSTP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,STPMIN, + STPMAX,INFO) C -------------------------------------------------------- C NEWSTP: Compute safeguarded step for linesearch; update C interval of uncertainty C -------------------------------------------------------- C argument-list variables: C .. Scalar Arguments .. DOUBLE PRECISION DP,DX,DY,FP,FX,FY,STP,STPMAX,STPMIN,STX,STY INTEGER INFO LOGICAL BRACKT C .. C .. Local Scalars .. DOUBLE PRECISION GAMMA,P,Q,R,S,SGND,STPC,STPF,STPQ,THETA,XSAFE LOGICAL BOUND C .. C .. Intrinsic Functions .. INTRINSIC ABS,MAX,MIN,SQRT C .. C-- xie C other variables: C .. Scalars in Common .. DOUBLE PRECISION TAU INTEGER MC,SRLS C .. C .. Common blocks .. COMMON /TAUINPUT/TAU,SRLS,MC C .. XSAFE = 1.0D-3 c----------------------- INFO = 0 C -------------------------------------------------------- C CHECK THE INPUT PARAMETERS FOR ERRORS. C -------------------------------------------------------- IF ((BRACKT.AND. (STP.LE.MIN(STX,STY).OR.STP.GE.MAX(STX, + STY))) .OR. DX* (STP-STX).GE.0.0 .OR. STPMAX.LT.STPMIN) RETURN C -------------------------------------------------------- C DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN. C -------------------------------------------------------- SGND = DP* (DX/ABS(DX)) C -------------------------------------------------------- C FIRST CASE. A HIGHER FUNCTION VALUE. C THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER C TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN, C ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN. C -------------------------------------------------------- IF (FP.GT.FX) THEN INFO = 1 BOUND = .TRUE. THETA = 3* (FX-FP)/ (STP-STX) + DX + DP S = MAX(ABS(THETA),ABS(DX),ABS(DP)) GAMMA = S*SQRT((THETA/S)**2- (DX/S)* (DP/S)) IF (STP.LT.STX) GAMMA = -GAMMA P = (GAMMA-DX) + THETA Q = ((GAMMA-DX)+GAMMA) + DP R = P/Q STPC = STX + R* (STP-STX) STPQ = STX + ((DX/ ((FX-FP)/ (STP-STX)+DX))/2)* (STP-STX) IF (ABS(STPC-STX).LT.ABS(STPQ-STX)) THEN STPF = STPC ELSE STPF = STPC + (STPQ-STPC)/2 END IF c--- xie: to avoid a too small STPE IF (SRLS.EQ.1) GO TO 10 IF (STP.GT.STX) THEN STPF = MAX(STX+XSAFE* (STP-STX),STPF) ELSE STPF = MIN(STX+XSAFE* (STP-STX),STPF) END IF c--------------------------------------------------------------- 10 BRACKT = .TRUE. C -------------------------------------------------------- C SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF C OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC C STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP, C THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN. C -------------------------------------------------------- ELSE IF (SGND.LT.0.0) THEN INFO = 2 BOUND = .FALSE. THETA = 3* (FX-FP)/ (STP-STX) + DX + DP S = MAX(ABS(THETA),ABS(DX),ABS(DP)) GAMMA = S*SQRT((THETA/S)**2- (DX/S)* (DP/S)) IF (STP.GT.STX) GAMMA = -GAMMA P = (GAMMA-DP) + THETA Q = ((GAMMA-DP)+GAMMA) + DX R = P/Q STPC = STP + R* (STX-STP) STPQ = STP + (DP/ (DP-DX))* (STX-STP) IF (ABS(STPC-STP).GT.ABS(STPQ-STP)) THEN STPF = STPC ELSE STPF = STPQ END IF BRACKT = .TRUE. C -------------------------------------------------------- C THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE C SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES. C THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY C IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC C IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE C EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO C COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP C CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN. C -------------------------------------------------------- ELSE IF (ABS(DP).LT.ABS(DX)) THEN INFO = 3 BOUND = .TRUE. THETA = 3* (FX-FP)/ (STP-STX) + DX + DP S = MAX(ABS(THETA),ABS(DX),ABS(DP)) C -------------------------------------------------------- C THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND C TO INFINITY IN THE DIRECTION OF THE STEP. C -------------------------------------------------------- GAMMA = S*SQRT(MAX(0.0D0, (THETA/S)**2- (DX/S)* (DP/S))) IF (STP.GT.STX) GAMMA = -GAMMA P = (GAMMA-DP) + THETA Q = (GAMMA+ (DX-DP)) + GAMMA R = P/Q IF (R.LT.0.0 .AND. GAMMA.NE.0.0) THEN STPC = STP + R* (STX-STP) ELSE IF (STP.GT.STX) THEN STPC = STPMAX ELSE STPC = STPMIN END IF STPQ = STP + (DP/ (DP-DX))* (STX-STP) IF (BRACKT) THEN IF (ABS(STP-STPC).LT.ABS(STP-STPQ)) THEN STPF = STPC ELSE STPF = STPQ END IF ELSE IF (ABS(STP-STPC).GT.ABS(STP-STPQ)) THEN STPF = STPC ELSE STPF = STPQ END IF END IF C -------------------------------------------------------- C FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE C SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES C NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP C IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN. C -------------------------------------------------------- ELSE INFO = 4 BOUND = .FALSE. IF (BRACKT) THEN THETA = 3* (FP-FY)/ (STY-STP) + DY + DP S = MAX(ABS(THETA),ABS(DY),ABS(DP)) GAMMA = S*SQRT((THETA/S)**2- (DY/S)* (DP/S)) IF (STP.GT.STY) GAMMA = -GAMMA P = (GAMMA-DP) + THETA Q = ((GAMMA-DP)+GAMMA) + DY R = P/Q STPC = STP + R* (STY-STP) STPF = STPC ELSE IF (STP.GT.STX) THEN STPF = STPMAX ELSE STPF = STPMIN END IF END IF C -------------------------------------------------------- C UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT C DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE. C -------------------------------------------------------- c--- xie: according the scheme, SGND is defined by c SGND = DP*(STX - STP) c But in the original code it was missed, leading to c SGND = DP*(DX/ABS(DX)) c So we add it as blow SGND = DP* (STX-STP) c-------------------------------- IF (FP.GT.FX) THEN STY = STP FY = FP DY = DP ELSE IF (SGND.LT.0.0) THEN STY = STX FY = FX DY = DX END IF STX = STP FX = FP DX = DP END IF C -------------------------------------------------------- C COMPUTE THE NEW STEP AND SAFEGUARD IT. C -------------------------------------------------------- STPF = MIN(STPMAX,STPF) STPF = MAX(STPMIN,STPF) STP = STPF IF (BRACKT .AND. BOUND) THEN IF (STY.GT.STX) THEN STP = MIN(STX+0.66* (STY-STX),STP) ELSE STP = MAX(STX+0.66* (STY-STX),STP) END IF END IF RETURN END C*********************************************************************** C*********************************************************************** DOUBLE PRECISION FUNCTION DMACH(JOB) C -------------------------------------------------------- C SMACH COMPUTES MACHINE PARAMETERS OF FLOATING POINT C ARITHMETIC FOR USE IN TESTING ONLY. NOT REQUIRED BY C LINPACK PROPER. C IF TROUBLE WITH AUTOMATIC COMPUTATION OF THESE QUANTITIES, C THEY CAN BE SET BY DIRECT ASSIGNMENT STATEMENTS. C ASSUME THE COMPUTER HAS C B = BASE OF ARITHMETIC C T = NUMBER OF BASE B DIGITS C L = SMALLEST POSSIBLE EXPONENT C U = LARGEST POSSIBLE EXPONENT C THEN C EPS = B**(1-T) C TINY = 100.0*B**(-L+T) C HUGE = 0.01*B**(U-T) C DMACH SAME AS SMACH EXCEPT T, L, U APPLY TO C DOUBLE PRECISION. C CMACH SAME AS SMACH EXCEPT IF COMPLEX DIVISION C IS DONE BY C 1/(X+I*Y) = (X-I*Y)/(X**2+Y**2) C THEN C TINY = SQRT(TINY) C HUGE = SQRT(HUGE) C JOB IS 1, 2 OR 3 FOR EPSILON, TINY AND HUGE, RESPECTIVELY. C -------------------------------------------------------- C .. Scalar Arguments .. INTEGER JOB C .. C .. Local Scalars .. DOUBLE PRECISION EPS,HUGE,S,TINY C .. EPS = 1.0D0 10 EPS = EPS/2.0D0 S = 1.0D0 + EPS IF (S.GT.1.0D0) GO TO 10 EPS = 2.0D0*EPS S = 1.0D0 20 TINY = S S = S/16.0D0 IF (S*1.0.NE.0.0D0) GO TO 20 TINY = (TINY/EPS)*100.0 HUGE = 1.0D0/TINY IF (JOB.EQ.1) DMACH = EPS IF (JOB.EQ.2) DMACH = TINY IF (JOB.EQ.3) DMACH = HUGE RETURN END C*********************************************************************** C 1/15/81 C*********************************************************************** C ODRV -- DRIVER FOR SPARSE MATRIX REORDERING ROUTINES C*********************************************************************** SUBROUTINE ODRV(N,IA,JA,A,P,IP,NSP,ISP,PATH,FLAG) C C DESCRIPTION C C ODRV FINDS A MINIMUM DEGREE ORDERING OF THE ROWS AND COLUMNS OF A C SYMMETRIC MATRIX M STORED IN (IA,JA,A) FORMAT (SEE BELOW). FOR THE C REORDERED MATRIX, THE WORK AND STORAGE REQUIRED TO PERFORM GAUSSIAN C ELIMINATION IS (USUALLY) SIGNIFICANTLY LESS. C C IF ONLY THE NONZERO ENTRIES IN THE UPPER TRIANGLE OF M ARE BEING C STORED, THEN ODRV SYMMETRICALLY REORDERS (IA,JA,A), (OPTIONALLY) C WITH THE DIAGONAL ENTRIES PLACED FIRST IN EACH ROW. THIS IS TO C ENSURE THAT IF M(I,J) WILL BE IN THE UPPER TRIANGLE OF M WITH C RESPECT TO THE NEW ORDERING, THEN M(I,J) IS STORED IN ROW I (AND C THUS M(J,I) IS NOT STORED); WHEREAS IF M(I,J) WILL BE IN THE C STRICT LOWER TRIANGLE OF M, THEN M(J,I) IS STORED IN ROW J (AND C THUS M(I,J) IS NOT STORED). C C C STORAGE OF SPARSE MATRICES C C THE NONZERO ENTRIES OF THE MATRIX M ARE STORED ROW-BY-ROW IN THE C ARRAY A. TO IDENTIFY THE INDIVIDUAL NONZERO ENTRIES IN EACH ROW, C WE NEED TO KNOW IN WHICH COLUMN EACH ENTRY LIES. THESE COLUMN C INDICES ARE STORED IN THE ARRAY JA; I.E., IF A(K) = M(I,J), THEN C JA(K) = J. TO IDENTIFY THE INDIVIDUAL ROWS, WE NEED TO KNOW WHERE C EACH ROW STARTS. THESE ROW POINTERS ARE STORED IN THE ARRAY IA; C I.E., IF M(I,J) IS THE FIRST NONZERO ENTRY (STORED) IN THE I-TH ROW C AND A(K) = M(I,J), THEN IA(I) = K. MOREOVER, IA(N+1) POINTS TO C THE FIRST LOCATION FOLLOWING THE LAST ELEMENT IN THE LAST ROW. C THUS, THE NUMBER OF ENTRIES IN THE I-TH ROW IS IA(I+1) - IA(I), C THE NONZERO ENTRIES IN THE I-TH ROW ARE STORED CONSECUTIVELY IN C C A(IA(I)), A(IA(I)+1), ..., A(IA(I+1)-1), C C AND THE CORRESPONDING COLUMN INDICES ARE STORED CONSECUTIVELY IN C C JA(IA(I)), JA(IA(I)+1), ..., JA(IA(I+1)-1). C C SINCE THE COEFFICIENT MATRIX IS SYMMETRIC, ONLY THE NONZERO ENTRIES C IN THE UPPER TRIANGLE NEED BE STORED. FOR EXAMPLE, THE MATRIX C C ( 1 0 2 3 0 ) C ( 0 4 0 0 0 ) C M = ( 2 0 5 6 0 ) C ( 3 0 6 7 8 ) C ( 0 0 0 8 9 ) C C COULD BE STORED AS C C \ 1 2 3 4 5 6 7 8 9 10 11 12 13 C ---+-------------------------------------- C IA \ 1 4 5 8 12 14 C JA \ 1 3 4 2 1 3 4 1 3 4 5 4 5 C A \ 1 2 3 4 2 5 6 3 6 7 8 8 9 C C OR (SYMMETRICALLY) AS C C \ 1 2 3 4 5 6 7 8 9 C ---+-------------------------- C IA \ 1 4 5 7 9 10 C JA \ 1 3 4 2 3 4 4 5 5 C A \ 1 2 3 4 5 6 7 8 9 . C C C PARAMETERS C C N - ORDER OF THE MATRIX C C IA - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING POINTERS TO DELIMIT C ROWS IN JA AND A; DIMENSION = N+1 C C JA - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING THE COLUMN INDICES C CORRESPONDING TO THE ELEMENTS OF A; DIMENSION = NUMBER OF C NONZERO ENTRIES IN (THE UPPER TRIANGLE OF) M C C A - REAL ONE-DIMENSIONAL ARRAY CONTAINING THE NONZERO ENTRIES IN C (THE UPPER TRIANGLE OF) M, STORED BY ROWS; DIMENSION = C NUMBER OF NONZERO ENTRIES IN (THE UPPER TRIANGLE OF) M C C P - INTEGER ONE-DIMENSIONAL ARRAY USED TO RETURN THE PERMUTATION C OF THE ROWS AND COLUMNS OF M CORRESPONDING TO THE MINIMUM C DEGREE ORDERING; DIMENSION = N C C IP - INTEGER ONE-DIMENSIONAL ARRAY USED TO RETURN THE INVERSE OF C THE PERMUTATION RETURNED IN P; DIMENSION = N C C NSP - DECLARED DIMENSION OF THE ONE-DIMENSIONAL ARRAY ISP; NSP C MUST BE AT LEAST 3N+4K, WHERE K IS THE NUMBER OF NONZEROES C IN THE STRICT UPPER TRIANGLE OF M C C ISP - INTEGER ONE-DIMENSIONAL ARRAY USED FOR WORKING STORAGE; C DIMENSION = NSP C C PATH - INTEGER PATH SPECIFICATION; VALUES AND THEIR MEANINGS ARE - C 1 FIND MINIMUM DEGREE ORDERING ONLY C 2 FIND MINIMUM DEGREE ORDERING AND REORDER SYMMETRICALLY C STORED MATRIX (USED WHEN ONLY THE NONZERO ENTRIES IN C THE UPPER TRIANGLE OF M ARE BEING STORED) C 3 REORDER SYMMETRICALLY STORED MATRIX AS SPECIFIED BY C INPUT PERMUTATION (USED WHEN AN ORDERING HAS ALREADY C BEEN DETERMINED AND ONLY THE NONZERO ENTRIES IN THE C UPPER TRIANGLE OF M ARE BEING STORED) C 4 SAME AS 2 BUT PUT DIAGONAL ENTRIES AT START OF EACH ROW C 5 SAME AS 3 BUT PUT DIAGONAL ENTRIES AT START OF EACH ROW C C FLAG - INTEGER ERROR FLAG; VALUES AND THEIR MEANINGS ARE - C 0 NO ERRORS DETECTED C 9N+K INSUFFICIENT STORAGE IN MD C 10N+1 INSUFFICIENT STORAGE IN ODRV C 11N+1 ILLEGAL PATH SPECIFICATION C C C CONVERSION FROM REAL TO DOUBLE PRECISION C C CHANGE THE REAL DECLARATIONS IN ODRV AND SRO TO DOUBLE PRECISION C DECLARATIONS. C C----------------------------------------------------------------------- C C.... REAL A(1) C C----INITIALIZE ERROR FLAG AND VALIDATE PATH SPECIFICATION C .. Scalar Arguments .. INTEGER FLAG,N,NSP,PATH C .. C .. Array Arguments .. DOUBLE PRECISION A(*) INTEGER IA(*),IP(*),ISP(*),JA(*),P(*) C .. C .. Local Scalars .. INTEGER HEAD,L,MAX,NEXT,Q,TMP,V LOGICAL DFLAG C .. C .. External Subroutines .. EXTERNAL MD,SRO C .. FLAG = 0 IF (PATH.LT.1 .OR. 5.LT.PATH) GO TO 50 C C----ALLOCATE STORAGE AND FIND MINIMUM DEGREE ORDERING IF ((PATH-1)* (PATH-2)* (PATH-4).NE.0) GO TO 10 MAX = (NSP-N)/2 V = 1 L = V + MAX HEAD = L + MAX NEXT = HEAD + N IF (MAX.LT.N) GO TO 40 C CALL MD(N,IA,JA,MAX,ISP(V),ISP(L),ISP(HEAD),P,IP,ISP(V),FLAG) IF (FLAG.NE.0) GO TO 30 C C----ALLOCATE STORAGE AND SYMMETRICALLY REORDER MATRIX 10 IF ((PATH-2)* (PATH-3)* (PATH-4)* (PATH-5).NE.0) GO TO 20 TMP = (NSP+1) - N Q = TMP - (IA(N+1)-1) IF (Q.LT.1) GO TO 40 C DFLAG = PATH .EQ. 4 .OR. PATH .EQ. 5 CALL SRO(N,IP,IA,JA,A,ISP(TMP),ISP(Q),DFLAG) C 20 RETURN C C ** ERROR -- ERROR DETECTED IN MD 30 RETURN C ** ERROR -- INSUFFICIENT STORAGE 40 FLAG = 10*N + 1 RETURN C ** ERROR -- ILLEGAL PATH SPECIFIED 50 FLAG = 11*N + 1 RETURN END C*********************************************************************** C SSF -- SYMBOLIC UT-D-U FACTORIZATION OF SPARSE SYMMETRIC MATRIX C*********************************************************************** SUBROUTINE SSF(N,P,IP,IA,JA,IJU,JU,IU,JUMAX,Q,MARK,JL,FLAG) C C ADDITIONAL PARAMETERS C C Q - INTEGER ONE-DIMENSIONAL WORK ARRAY; DIMENSION = N C C MARK - INTEGER ONE-DIMENSIONAL WORK ARRAY; DIMENSION = N C C JL - INTEGER ONE-DIMENSIONAL WORK ARRAY; DIMENSION = N C C C DEFINITIONS OF INTERNAL PARAMETERS (DURING K-TH STAGE OF ELIMINATION) C C Q CONTAINS AN ORDERED LINKED LIST REPRESENTATION OF THE NONZERO C STRUCTURE OF THE K-TH ROW OF U -- C Q(K) IS THE FIRST COLUMN WITH A NONZERO ENTRY C Q(I) IS THE NEXT COLUMN WITH A NONZERO ENTRY AFTER COLUMN I C IN EITHER CASE, Q(I) = N+1 INDICATES THE END OF THE LIST C C JL CONTAINS LISTS OF ROWS TO BE MERGED INTO UNELIMINATED ROWS -- C I GE K => JL(I) IS THE FIRST ROW TO BE MERGED INTO ROW I C I LT K => JL(I) IS THE ROW FOLLOWING ROW I IN SOME LIST OF ROWS C IN EITHER CASE, JL(I) = 0 INDICATES THE END OF A LIST C C MARK(I) IS THE LAST ROW STORED IN JU FOR WHICH U(MARK(I),I) NE 0 C C JUMIN AND JUPTR ARE THE INDICES IN JU OF THE FIRST AND LAST C ELEMENTS IN THE LAST ROW SAVED IN JU C C LUK IS THE NUMBER OF NONZERO ENTRIES IN THE K-TH ROW C C----------------------------------------------------------------------- C C C----INITIALIZATION C .. Scalar Arguments .. INTEGER FLAG,JUMAX,N C .. C .. Array Arguments .. INTEGER IA(*),IJU(*),IP(*),IU(*),JA(*),JL(*),JU(*),MARK(*),P(*), + Q(*) C .. C .. Local Scalars .. INTEGER I,J,JMAX,JMIN,JUMIN,JUPTR,K,LMAX,LUI,LUK,M,QM,TAG,VJ LOGICAL CLIQUE C .. JUMIN = 1 JUPTR = 0 IU(1) = 1 DO 10 K = 1,N MARK(K) = 0 JL(K) = 0 10 CONTINUE C C----FOR EACH ROW K DO 190 K = 1,N LUK = 0 Q(K) = N + 1 C TAG = MARK(K) CLIQUE = .FALSE. IF (JL(K).NE.0) CLIQUE = JL(JL(K)) .EQ. 0 C C------INITIALIZE NONZERO STRUCTURE OF K-TH ROW TO ROW P(K) OF M JMIN = IA(P(K)) JMAX = IA(P(K)+1) - 1 IF (JMIN.GT.JMAX) GO TO 40 DO 30 J = JMIN,JMAX VJ = IP(JA(J)) IF (VJ.LE.K) GO TO 30 C QM = K 20 M = QM QM = Q(M) IF (QM.LT.VJ) GO TO 20 IF (QM.EQ.VJ) GO TO 200 LUK = LUK + 1 Q(M) = VJ Q(VJ) = QM IF (MARK(VJ).NE.TAG) CLIQUE = .FALSE. C 30 CONTINUE C C------IF EXACTLY ONE ROW IS TO BE MERGED INTO THE K-TH ROW AND THERE IS C------A NONZERO ENTRY IN EVERY COLUMN IN THAT ROW IN WHICH THERE IS A C------NONZERO ENTRY IN ROW P(K) OF M, THEN DO NOT COMPUTE FILL-IN, JUST C------USE THE COLUMN INDICES FOR THE ROW WHICH WAS TO HAVE BEEN MERGED 40 IF (.NOT.CLIQUE) GO TO 50 IJU(K) = IJU(JL(K)) + 1 LUK = IU(JL(K)+1) - (IU(JL(K))+1) GO TO 170 C C------MODIFY NONZERO STRUCTURE OF K-TH ROW BY COMPUTING FILL-IN C------FOR EACH ROW I TO BE MERGED IN 50 LMAX = 0 IJU(K) = JUPTR C I = K 60 I = JL(I) IF (I.EQ.0) GO TO 100 C C--------MERGE ROW I INTO K-TH ROW LUI = IU(I+1) - (IU(I)+1) JMIN = IJU(I) + 1 JMAX = IJU(I) + LUI QM = K C DO 80 J = JMIN,JMAX VJ = JU(J) 70 M = QM QM = Q(M) IF (QM.LT.VJ) GO TO 70 IF (QM.EQ.VJ) GO TO 80 LUK = LUK + 1 Q(M) = VJ Q(VJ) = QM QM = VJ 80 CONTINUE C C--------REMEMBER LENGTH AND POSITION IN JU OF LONGEST ROW MERGED IF (LUI.LE.LMAX) GO TO 90 LMAX = LUI IJU(K) = JMIN C 90 GO TO 60 C C------IF THE K-TH ROW IS THE SAME LENGTH AS THE LONGEST ROW MERGED, C------THEN USE THE COLUMN INDICES FOR THAT ROW 100 IF (LUK.EQ.LMAX) GO TO 170 C C------IF THE TAIL OF THE LAST ROW SAVED IN JU IS THE SAME AS THE HEAD C------OF THE K-TH ROW, THEN OVERLAP THE TWO SETS OF COLUMN INDICES -- C--------SEARCH LAST ROW SAVED FOR FIRST NONZERO ENTRY IN K-TH ROW ... I = Q(K) IF (JUMIN.GT.JUPTR) GO TO 120 DO 110 JMIN = JUMIN,JUPTR IF (JU(JMIN)-I) 110,130,120 110 CONTINUE 120 GO TO 150 C C--------... AND THEN TEST WHETHER TAIL MATCHES HEAD OF K-TH ROW 130 IJU(K) = JMIN DO 140 J = JMIN,JUPTR IF (JU(J).NE.I) GO TO 150 I = Q(I) IF (I.GT.N) GO TO 170 140 CONTINUE JUPTR = JMIN - 1 C C------SAVE NONZERO STRUCTURE OF K-TH ROW IN JU 150 I = K JUMIN = JUPTR + 1 JUPTR = JUPTR + LUK IF (JUPTR.GT.JUMAX) GO TO 210 DO 160 J = JUMIN,JUPTR I = Q(I) JU(J) = I MARK(I) = K 160 CONTINUE IJU(K) = JUMIN C C------ADD K TO ROW LIST FOR FIRST NONZERO ELEMENT IN K-TH ROW 170 IF (LUK.LE.1) GO TO 180 I = JU(IJU(K)) JL(K) = JL(I) JL(I) = K C 180 IU(K+1) = IU(K) + LUK 190 CONTINUE C FLAG = 0 RETURN C C ** ERROR -- DUPLICATE ENTRY IN A 200 FLAG = 2*N + P(K) RETURN C ** ERROR -- INSUFFICIENT STORAGE FOR JU 210 FLAG = 6*N + K RETURN END C C*********************************************************************** C SNS -- SOLUTION OF SPARSE SYMMETRIC POSITIVE DEFINITE SYSTEM OF C LINEAR EQUATIONS MX = B GIVEN UT-D-U FACTORIZATION OF M C*********************************************************************** SUBROUTINE SNS(N,P,D,IJU,JU,IU,U,Z,B,TMP) C REAL D(*), U(*), Z(*), B(*), TMP(*), TMPK, SUM C C ADDITIONAL PARAMETERS C C TMP - REAL ONE-DIMENSIONAL WORK ARRAY; DIMENSION = N C C----------------------------------------------------------------------- C C----SET TMP TO PERMUTED B C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION B(*),D(*),TMP(*),U(*),Z(*) INTEGER IJU(*),IU(*),JU(*),P(*) C .. C .. Local Scalars .. DOUBLE PRECISION SUM,TMPK INTEGER I,J,JMAX,JMIN,K,MU C .. DO 10 K = 1,N TMP(K) = B(P(K)) 10 CONTINUE C C----SOLVE UT D Y = B BY FORWARD SUBSTITUTION DO 40 K = 1,N TMPK = TMP(K) JMIN = IU(K) JMAX = IU(K+1) - 1 IF (JMIN.GT.JMAX) GO TO 30 MU = IJU(K) - JMIN DO 20 J = JMIN,JMAX TMP(JU(MU+J)) = TMP(JU(MU+J)) + U(J)*TMPK 20 CONTINUE 30 TMP(K) = TMPK*D(K) 40 CONTINUE C C----SOLVE U X = Y BY BACK SUBSTITUTION K = N DO 70 I = 1,N SUM = TMP(K) JMIN = IU(K) JMAX = IU(K+1) - 1 IF (JMIN.GT.JMAX) GO TO 60 MU = IJU(K) - JMIN DO 50 J = JMIN,JMAX SUM = SUM + U(J)*TMP(JU(MU+J)) 50 CONTINUE 60 TMP(K) = SUM Z(P(K)) = SUM K = K - 1 70 CONTINUE C RETURN END C*********************************************************************** C*********************************************************************** C*********************************************************************** C SRO -- SYMMETRIC REORDERING OF SPARSE SYMMETRIC MATRIX C*********************************************************************** SUBROUTINE SRO(N,IP,IA,JA,A,Q,R,DFLAG) C C DESCRIPTION C C THE NONZERO ENTRIES OF THE MATRIX M ARE ASSUMED TO BE STORED C SYMMETRICALLY IN (IA,JA,A) FORMAT (I.E., NOT BOTH M(I,J) AND M(J,I) C ARE STORED IF I NE J). C C SRO DOES NOT REARRANGE THE ORDER OF THE ROWS, BUT DOES MOVE C NONZEROES FROM ONE ROW TO ANOTHER TO ENSURE THAT IF M(I,J) WILL BE C IN THE UPPER TRIANGLE OF M WITH RESPECT TO THE NEW ORDERING, THEN C M(I,J) IS STORED IN ROW I (AND THUS M(J,I) IS NOT STORED); WHEREAS C IF M(I,J) WILL BE IN THE STRICT LOWER TRIANGLE OF M, THEN M(J,I) IS C STORED IN ROW J (AND THUS M(I,J) IS NOT STORED). C C C ADDITIONAL PARAMETERS C C Q - INTEGER ONE-DIMENSIONAL WORK ARRAY; DIMENSION = N C C R - INTEGER ONE-DIMENSIONAL WORK ARRAY; DIMENSION = NUMBER OF C NONZERO ENTRIES IN THE UPPER TRIANGLE OF M C C DFLAG - LOGICAL VARIABLE; IF DFLAG = .TRUE., THEN STORE NONZERO C DIAGONAL ELEMENTS AT THE BEGINNING OF THE ROW C C----------------------------------------------------------------------- C C REAL A(*), AK C C C--PHASE 1 -- FIND ROW IN WHICH TO STORE EACH NONZERO C----INITIALIZE COUNT OF NONZEROES TO BE STORED IN EACH ROW C .. Scalar Arguments .. INTEGER N LOGICAL DFLAG C .. C .. Array Arguments .. DOUBLE PRECISION A(*) INTEGER IA(*),IP(*),JA(*),Q(*),R(*) C .. C .. Local Scalars .. DOUBLE PRECISION AK INTEGER I,ILAST,J,JAK,JDUMMY,JMAX,JMIN,K C .. DO 10 I = 1,N Q(I) = 0 10 CONTINUE C C----FOR EACH NONZERO ELEMENT A(J) DO 30 I = 1,N JMIN = IA(I) JMAX = IA(I+1) - 1 IF (JMIN.GT.JMAX) GO TO 30 DO 20 J = JMIN,JMAX C C--------FIND ROW (=R(J)) AND COLUMN (=JA(J)) IN WHICH TO STORE A(J) ... K = JA(J) IF (IP(K).LT.IP(I)) JA(J) = I IF (IP(K).GE.IP(I)) K = I R(J) = K C C--------... AND INCREMENT COUNT OF NONZEROES (=Q(R(J)) IN THAT ROW Q(K) = Q(K) + 1 20 CONTINUE 30 CONTINUE C C C--PHASE 2 -- FIND NEW IA AND PERMUTATION TO APPLY TO (JA,A) C----DETERMINE POINTERS TO DELIMIT ROWS IN PERMUTED (JA,A) DO 40 I = 1,N IA(I+1) = IA(I) + Q(I) Q(I) = IA(I+1) 40 CONTINUE C C----DETERMINE WHERE EACH (JA(J),A(J)) IS STORED IN PERMUTED (JA,A) C----FOR EACH NONZERO ELEMENT (IN REVERSE ORDER) ILAST = 0 JMIN = IA(1) JMAX = IA(N+1) - 1 J = JMAX DO 70 JDUMMY = JMIN,JMAX I = R(J) IF (.NOT.DFLAG .OR. JA(J).NE.I .OR. I.EQ.ILAST) GO TO 50 C C------IF DFLAG, THEN PUT DIAGONAL NONZERO AT BEGINNING OF ROW R(J) = IA(I) ILAST = I GO TO 60 C C------PUT (OFF-DIAGONAL) NONZERO IN LAST UNUSED LOCATION IN ROW 50 Q(I) = Q(I) - 1 R(J) = Q(I) C 60 J = J - 1 70 CONTINUE C C C--PHASE 3 -- PERMUTE (JA,A) TO UPPER TRIANGULAR FORM (WRT NEW ORDERING) DO 90 J = JMIN,JMAX 80 IF (R(J).EQ.J) GO TO 90 K = R(J) R(J) = R(K) R(K) = K JAK = JA(K) JA(K) = JA(J) JA(J) = JAK AK = A(K) A(K) = A(J) A(J) = AK GO TO 80 90 CONTINUE C RETURN END C C*********************************************************************** C*********************************************************************** C MD -- MINIMUM DEGREE ALGORITHM (BASED ON ELEMENT MODEL) C*********************************************************************** SUBROUTINE MD(N,IA,JA,MAX,V,L,HEAD,LAST,NEXT,MARK,FLAG) C C DESCRIPTION C C MD FINDS A MINIMUM DEGREE ORDERING OF THE ROWS AND COLUMNS OF A C SYMMETRIC MATRIX M STORED IN (IA,JA,A) FORMAT. C C C ADDITIONAL PARAMETERS C C MAX - DECLARED DIMENSION OF THE ONE-DIMENSIONAL ARRAYS V AND L; C MAX MUST BE AT LEAST N+2K, WHERE K IS THE NUMBER OF C NONZEROES IN THE STRICT UPPER TRIANGLE OF M C C V - INTEGER ONE-DIMENSIONAL WORK ARRAY; DIMENSION = MAX C C L - INTEGER ONE-DIMENSIONAL WORK ARRAY; DIMENSION = MAX C C HEAD - INTEGER ONE-DIMENSIONAL WORK ARRAY; DIMENSION = N C C LAST - INTEGER ONE-DIMENSIONAL ARRAY USED TO RETURN THE PERMUTATION C OF THE ROWS AND COLUMNS OF M CORRESPONDING TO THE MINIMUM C DEGREE ORDERING; DIMENSION = N C C NEXT - INTEGER ONE-DIMENSIONAL ARRAY USED TO RETURN THE INVERSE OF C THE PERMUTATION RETURNED IN LAST; DIMENSION = N C C MARK - INTEGER ONE-DIMENSIONAL WORK ARRAY (MAY BE THE SAME AS V); C DIMENSION = N C C FLAG - INTEGER ERROR FLAG; VALUES AND THEIR MEANINGS ARE - C 0 NO ERRORS DETECTED C 11N+1 INSUFFICIENT STORAGE IN MD C C C DEFINITIONS OF INTERNAL PARAMETERS C C ---------+--------------------------------------------------------- C V(S) \ VALUE FIELD OF LIST ENTRY C ---------+--------------------------------------------------------- C L(S) \ LINK FIELD OF LIST ENTRY (0 => END OF LIST) C ---------+--------------------------------------------------------- C L(VI) \ POINTER TO ELEMENT LIST OF UNELIMINATED VERTEX VI C ---------+--------------------------------------------------------- C L(EJ) \ POINTER TO BOUNDARY LIST OF ACTIVE ELEMENT EJ C ---------+--------------------------------------------------------- C HEAD(D) \ VJ => VJ HEAD OF D-LIST D C \ 0 => NO VERTEX IN D-LIST D C C C \ VI UNELIMINATED VERTEX C \ VI IN EK \ VI NOT IN EK C ---------+-----------------------------+--------------------------- C NEXT(VI) \ UNDEFINED BUT NONNEGATIVE \ VJ => VJ NEXT IN D-LIST C \ \ 0 => VI TAIL OF D-LIST C ---------+-----------------------------+--------------------------- C LAST(VI) \ (NOT SET UNTIL MDP) \ -D => VI HEAD OF D-LIST D C \-VK => COMPUTE DEGREE \ VJ => VJ LAST IN D-LIST C \ EJ => VI PROTOTYPE OF EJ \ 0 => VI NOT IN ANY D-LIST C \ 0 => DO NOT COMPUTE DEGREE C ---------+-----------------------------+--------------------------- C MARK(VI) \ MARK(VK) \ NONNEGATIVE TAG < MARK(VK) C C C \ VI ELIMINATED VERTEX C \ EI ACTIVE ELEMENT \ OTHERWISE C ---------+-----------------------------+--------------------------- C NEXT(VI) \ -J => VI WAS J-TH VERTEX \ -J => VI WAS J-TH VERTEX C \ TO BE ELIMINATED \ TO BE ELIMINATED C ---------+-----------------------------+--------------------------- C LAST(VI) \ M => SIZE OF EI = M \ UNDEFINED C ---------+-----------------------------+--------------------------- C MARK(VI) \ -M => OVERLAP COUNT OF EI \ UNDEFINED C \ WITH EK = M C \ OTHERWISE NONNEGATIVE TAG C \ < MARK(VK) C C----------------------------------------------------------------------- C C C----INITIALIZATION C .. Scalar Arguments .. INTEGER FLAG,MAX,N C .. C .. Array Arguments .. INTEGER HEAD(*),IA(*),JA(*),L(*),LAST(*),MARK(*),NEXT(*),V(*) C .. C .. Local Scalars .. INTEGER DMIN,EK,K,TAG,TAIL,VK C .. C .. External Subroutines .. EXTERNAL MDI,MDM,MDP,MDU C .. C .. Equivalences .. EQUIVALENCE (VK,EK) C .. TAG = 0 CALL MDI(N,IA,JA,MAX,V,L,HEAD,LAST,NEXT,MARK,TAG,FLAG) IF (FLAG.NE.0) RETURN C K = 0 DMIN = 1 C C----WHILE K < N DO 10 IF (K.GE.N) GO TO 40 C C------SEARCH FOR VERTEX OF MINIMUM DEGREE 20 IF (HEAD(DMIN).GT.0) GO TO 30 DMIN = DMIN + 1 GO TO 20 C C------REMOVE VERTEX VK OF MINIMUM DEGREE FROM DEGREE LIST 30 VK = HEAD(DMIN) HEAD(DMIN) = NEXT(VK) IF (HEAD(DMIN).GT.0) LAST(HEAD(DMIN)) = -DMIN C C------NUMBER VERTEX VK, ADJUST TAG, AND TAG VK K = K + 1 NEXT(VK) = -K LAST(EK) = DMIN - 1 TAG = TAG + LAST(EK) MARK(VK) = TAG C C------FORM ELEMENT EK FROM UNELIMINATED NEIGHBORS OF VK CALL MDM(VK,TAIL,V,L,LAST,NEXT,MARK) C C------PURGE INACTIVE ELEMENTS AND DO MASS ELIMINATION CALL MDP(K,EK,TAIL,V,L,HEAD,LAST,NEXT,MARK) C C------UPDATE DEGREES OF UNELIMINATED VERTICES IN EK CALL MDU(EK,DMIN,V,L,HEAD,LAST,NEXT,MARK) C GO TO 10 C C----GENERATE INVERSE PERMUTATION FROM PERMUTATION 40 DO 50 K = 1,N NEXT(K) = -NEXT(K) LAST(NEXT(K)) = K 50 CONTINUE C RETURN END C C*********************************************************************** C MDI -- INITIALIZATION C*********************************************************************** SUBROUTINE MDI(N,IA,JA,MAX,V,L,HEAD,LAST,NEXT,MARK,TAG,FLAG) C C----INITIALIZE DEGREES, ELEMENT LISTS, AND DEGREE LISTS C .. Scalar Arguments .. INTEGER FLAG,MAX,N,TAG C .. C .. Array Arguments .. INTEGER HEAD(*),IA(*),JA(*),L(*),LAST(*),MARK(*),NEXT(*),V(*) C .. C .. Local Scalars .. INTEGER DVI,J,JMAX,JMIN,SFS,VI,VJ C .. DO 10 VI = 1,N MARK(VI) = 1 L(VI) = 0 HEAD(VI) = 0 10 CONTINUE SFS = N + 1 C C----CREATE NONZERO STRUCTURE C----FOR EACH NONZERO ENTRY A(VI,VJ) IN STRICT UPPER TRIANGLE DO 30 VI = 1,N JMIN = IA(VI) JMAX = IA(VI+1) - 1 IF (JMIN.GT.JMAX) GO TO 30 DO 20 J = JMIN,JMAX VJ = JA(J) IF (VI.GE.VJ) GO TO 20 IF (SFS.GE.MAX) GO TO 50 C C------ENTER VJ IN ELEMENT LIST FOR VI MARK(VI) = MARK(VI) + 1 V(SFS) = VJ L(SFS) = L(VI) L(VI) = SFS SFS = SFS + 1 C C------ENTER VI IN ELEMENT LIST FOR VJ MARK(VJ) = MARK(VJ) + 1 V(SFS) = VI L(SFS) = L(VJ) L(VJ) = SFS SFS = SFS + 1 20 CONTINUE 30 CONTINUE C C----CREATE DEGREE LISTS AND INITIALIZE MARK VECTOR DO 40 VI = 1,N DVI = MARK(VI) NEXT(VI) = HEAD(DVI) HEAD(DVI) = VI LAST(VI) = -DVI IF (NEXT(VI).GT.0) LAST(NEXT(VI)) = VI MARK(VI) = TAG 40 CONTINUE C RETURN C C ** ERROR -- INSUFFICIENT STORAGE 50 FLAG = 9*N + VI RETURN END C C*********************************************************************** C MDM -- FORM ELEMENT FROM UNELIMINATED NEIGHBORS OF VK C*********************************************************************** SUBROUTINE MDM(VK,TAIL,V,L,LAST,NEXT,MARK) C C----INITIALIZE TAG AND LIST OF UNELIMINATED NEIGHBORS C .. Scalar Arguments .. INTEGER TAIL,VK C .. C .. Array Arguments .. INTEGER L(*),LAST(*),MARK(*),NEXT(*),V(*) C .. C .. Local Scalars .. INTEGER B,BLP,BLPMAX,ES,LB,LS,S,TAG,VB,VS C .. C .. Equivalences .. EQUIVALENCE (VS,ES) C .. TAG = MARK(VK) TAIL = VK C C----FOR EACH VERTEX/ELEMENT VS/ES IN ELEMENT LIST OF VK LS = L(VK) 10 S = LS IF (S.EQ.0) GO TO 50 LS = L(S) VS = V(S) IF (NEXT(VS).LT.0) GO TO 20 C C------IF VS IS UNELIMINATED VERTEX, THEN TAG AND APPEND TO LIST OF C------UNELIMINATED NEIGHBORS MARK(VS) = TAG L(TAIL) = S TAIL = S GO TO 40 C C------IF ES IS ACTIVE ELEMENT, THEN ... C--------FOR EACH VERTEX VB IN BOUNDARY LIST OF ELEMENT ES 20 LB = L(ES) BLPMAX = LAST(ES) DO 30 BLP = 1,BLPMAX B = LB LB = L(B) VB = V(B) C C----------IF VB IS UNTAGGED VERTEX, THEN TAG AND APPEND TO LIST OF C----------UNELIMINATED NEIGHBORS IF (MARK(VB).GE.TAG) GO TO 30 MARK(VB) = TAG L(TAIL) = B TAIL = B 30 CONTINUE C C--------MARK ES INACTIVE MARK(ES) = TAG C 40 GO TO 10 C C----TERMINATE LIST OF UNELIMINATED NEIGHBORS 50 L(TAIL) = 0 C RETURN END C C*********************************************************************** C MDP -- PURGE INACTIVE ELEMENTS AND DO MASS ELIMINATION C*********************************************************************** SUBROUTINE MDP(K,EK,TAIL,V,L,HEAD,LAST,NEXT,MARK) C C----INITIALIZE TAG C .. Scalar Arguments .. INTEGER EK,K,TAIL C .. C .. Array Arguments .. INTEGER HEAD(*),L(*),LAST(*),MARK(*),NEXT(*),V(*) C .. C .. Local Scalars .. INTEGER ES,EVI,FREE,I,ILP,ILPMAX,LI,LS,LVI,S,TAG,VI C .. TAG = MARK(EK) C C----FOR EACH VERTEX VI IN EK LI = EK ILPMAX = LAST(EK) IF (ILPMAX.LE.0) GO TO 120 DO 110 ILP = 1,ILPMAX I = LI LI = L(I) VI = V(LI) C C------REMOVE VI FROM DEGREE LIST IF (LAST(VI).EQ.0) GO TO 30 IF (LAST(VI).GT.0) GO TO 10 HEAD(-LAST(VI)) = NEXT(VI) GO TO 20 10 NEXT(LAST(VI)) = NEXT(VI) 20 IF (NEXT(VI).GT.0) LAST(NEXT(VI)) = LAST(VI) C C------REMOVE INACTIVE ITEMS FROM ELEMENT LIST OF VI 30 LS = VI 40 S = LS LS = L(S) IF (LS.EQ.0) GO TO 60 ES = V(LS) IF (MARK(ES).LT.TAG) GO TO 50 FREE = LS L(S) = L(LS) LS = S 50 GO TO 40 C C------IF VI IS INTERIOR VERTEX, THEN REMOVE FROM LIST AND ELIMINATE 60 LVI = L(VI) IF (LVI.NE.0) GO TO 70 L(I) = L(LI) LI = I C K = K + 1 NEXT(VI) = -K LAST(EK) = LAST(EK) - 1 GO TO 110 C C------ELSE ... C--------CLASSIFY VERTEX VI 70 IF (L(LVI).NE.0) GO TO 90 EVI = V(LVI) IF (NEXT(EVI).GE.0) GO TO 90 IF (MARK(EVI).LT.0) GO TO 80 C C----------IF VI IS PROTOTYPE VERTEX, THEN MARK AS SUCH, INITIALIZE C----------OVERLAP COUNT FOR CORRESPONDING ELEMENT, AND MOVE VI TO END C----------OF BOUNDARY LIST LAST(VI) = EVI MARK(EVI) = -1 L(TAIL) = LI TAIL = LI L(I) = L(LI) LI = I GO TO 100 C C----------ELSE IF VI IS DUPLICATE VERTEX, THEN MARK AS SUCH AND ADJUST C----------OVERLAP COUNT FOR CORRESPONDING ELEMENT 80 LAST(VI) = 0 MARK(EVI) = MARK(EVI) - 1 GO TO 100 C C----------ELSE MARK VI TO COMPUTE DEGREE 90 LAST(VI) = -EK C C--------INSERT EK IN ELEMENT LIST OF VI 100 V(FREE) = EK L(FREE) = L(VI) L(VI) = FREE 110 CONTINUE C C----TERMINATE BOUNDARY LIST 120 L(TAIL) = 0 C RETURN END C C*********************************************************************** C MDU -- UPDATE DEGREES OF UNELIMINATED VERTICES IN EK C*********************************************************************** SUBROUTINE MDU(EK,DMIN,V,L,HEAD,LAST,NEXT,MARK) C C----INITIALIZE TAG C .. Scalar Arguments .. INTEGER DMIN,EK C .. C .. Array Arguments .. INTEGER HEAD(*),L(*),LAST(*),MARK(*),NEXT(*),V(*) C .. C .. Local Scalars .. INTEGER B,BLP,BLPMAX,DVI,ES,EVI,I,ILP,ILPMAX,S,TAG,VB,VI,VS C .. C .. Equivalences .. EQUIVALENCE (VS,ES) C .. TAG = MARK(EK) - LAST(EK) C C----FOR EACH VERTEX VI IN EK I = EK ILPMAX = LAST(EK) IF (ILPMAX.LE.0) GO TO 110 DO 100 ILP = 1,ILPMAX I = L(I) VI = V(I) IF (LAST(VI)) 10,100,80 C C------IF VI NEITHER PROTOTYPE NOR DUPLICATE VERTEX, THEN MERGE ELEMENTS C------TO COMPUTE DEGREE 10 TAG = TAG + 1 DVI = LAST(EK) C C--------FOR EACH VERTEX/ELEMENT VS/ES IN ELEMENT LIST OF VI S = L(VI) 20 S = L(S) IF (S.EQ.0) GO TO 90 VS = V(S) IF (NEXT(VS).LT.0) GO TO 30 C C----------IF VS IS UNELIMINATED VERTEX, THEN TAG AND ADJUST DEGREE MARK(VS) = TAG DVI = DVI + 1 GO TO 50 C C----------IF ES IS ACTIVE ELEMENT, THEN EXPAND C------------CHECK FOR OUTMATCHED VERTEX 30 IF (MARK(ES).LT.0) GO TO 60 C C------------FOR EACH VERTEX VB IN ES B = ES BLPMAX = LAST(ES) DO 40 BLP = 1,BLPMAX B = L(B) VB = V(B) C C--------------IF VB IS UNTAGGED, THEN TAG AND ADJUST DEGREE IF (MARK(VB).GE.TAG) GO TO 40 MARK(VB) = TAG DVI = DVI + 1 40 CONTINUE C 50 GO TO 20 C C------ELSE IF VI IS OUTMATCHED VERTEX, THEN ADJUST OVERLAPS BUT DO NOT C------COMPUTE DEGREE 60 LAST(VI) = 0 MARK(ES) = MARK(ES) - 1 70 S = L(S) IF (S.EQ.0) GO TO 100 ES = V(S) IF (MARK(ES).LT.0) MARK(ES) = MARK(ES) - 1 GO TO 70 C C------ELSE IF VI IS PROTOTYPE VERTEX, THEN CALCULATE DEGREE BY C------INCLUSION/EXCLUSION AND RESET OVERLAP COUNT 80 EVI = LAST(VI) DVI = LAST(EK) + LAST(EVI) + MARK(EVI) MARK(EVI) = 0 C C------INSERT VI IN APPROPRIATE DEGREE LIST 90 NEXT(VI) = HEAD(DVI) HEAD(DVI) = VI LAST(VI) = -DVI IF (NEXT(VI).GT.0) LAST(NEXT(VI)) = VI IF (DVI.LT.DMIN) DMIN = DVI C 100 CONTINUE C 110 RETURN END C*********************************************************************** C*********************************************************************** C Modified Routines of the Yale Sparse Matrix Package (YSMP): C C SDRVMD - a modified form of the SDRV driver. C SNFMOD - a modified form of the SNF routine. C C copyright (c) 1990 by Tamar Schlick C C*********************************************************************** C Yale's SDRV solves a linear system for (symmetric) positive- C definite matrices. It calls the following routines: C SSF (for symbolic factorization) C SNF (for numerical factorization) and C SNS (for numerical solution). C Our goal is to solve large sparse symmetric linear systems for C matrices that are not necessarily pos-def. Thus, we C replace SNF by SNFMOD so that a modified-Cholesky (MCF), rather than C a Cholesky, factorization is performed. In SDRV, we replace C the statement "CALL SNF" by "CALL SNFMOD". C In Yale's SDRV, the diagnostic parameter FLAG is set to zero in C the non pos-def case (and control is returned to the main program). C Here, instead, we set FLAG in SNFMOD to a negative integer if the C matrix is not sufficiently pos-def. C Specifically, FLAG is set to minus C the position of the diagonal element in the original matrix whose C modification was of the largest magnitude. Recall that in MCF we C produce matrices E,D, and U so that C M + E = UT-D-U where E and D are C diagonal and U is unit upper-triangular. FLAG records the index k C for the largest modification in E: ( E(P(k)) = max over i {E(i)} ). C C All modifications to the original YSMP code are indicated. C 1/15/81 C*********************************************************************** C SDRV -- DRIVER FOR SPARSE SYMMETRIC POSITIVE DEFINITE MATRIX ROUTINES C*********************************************************************** C ==================== change #1 (replacement) =====================1 C WAS: SUBROUTINE SDRV C ===================================================================== C ==================================================================end SUBROUTINE SDRVMD(N,P,IP,IA,JA,A,B,Z,NSP,ISP,RSP,ESP,PATH,FLAG) C C DESCRIPTION C C ==================== change #2 (replacement) =====================2 C WAS: SDRV SOLVES SPARSE SYMMETRIC POSITIVE DEFINITE SYSTEMS OF LINEAR C ===================================================================== C SDRVMD SOLVES SPARSE SYMMETRIC SYSTEMS OF LINEAR C ==================================================================end C EQUATIONS. THE SOLUTION PROCESS IS DIVIDED INTO THREE STAGES -- C C SSF - THE COEFFICIENT MATRIX M IS FACTORED SYMBOLICALLY TO C DETERMINE WHERE FILLIN WILL OCCUR DURING THE NUMERIC C FACTORIZATION. C C ==================== change #3 (replacement) =====================3 C WAS: SNF - M IS FACTORED NUMERICALLY INTO THE PRODUCT UT-D-U, WHERE C ===================================================================== C SNFMOD - M+E IS FACTORED NUMERICALLY BY THE GILL/MURRAY/WRIGHT C MODIFIED CHOLESKY FACTORIZATION: M + E = UT-D-U, WHERE C E IS DIAGONAL, C ==================================================================end C D IS DIAGONAL AND U IS UNIT UPPER TRIANGULAR. C C SNS - THE LINEAR SYSTEM MX = B IS SOLVED USING THE UT-D-U C FACTORIZATION FROM SNF. C C FOR SEVERAL SYSTEMS WITH THE SAME COEFFICIENT MATRIX, SSF AND SNF C NEED BE DONE ONLY ONCE (FOR THE FIRST SYSTEM); THEN SNS IS DONE C ONCE FOR EACH ADDITIONAL RIGHT-HAND SIDE. FOR SEVERAL SYSTEMS C WHOSE COEFFICIENT MATRICES HAVE THE SAME NONZERO STRUCTURE, SSF C NEED BE DONE ONLY ONCE (FOR THE FIRST SYSTEM); THEN SNF AND SNS C ARE DONE ONCE FOR EACH ADDITIONAL SYSTEM. C C C STORAGE OF SPARSE MATRICES C C THE NONZERO ENTRIES OF THE MATRIX M ARE STORED ROW-BY-ROW IN THE C ARRAY A. TO IDENTIFY THE INDIVIDUAL NONZERO ENTRIES IN EACH ROW, C WE NEED TO KNOW IN WHICH COLUMN EACH ENTRY LIES. THESE COLUMN C INDICES ARE STORED IN THE ARRAY JA; I.E., IF A(K) = M(I,J), THEN C JA(K) = J. TO IDENTIFY THE INDIVIDUAL ROWS, WE NEED TO KNOW WHERE C EACH ROW STARTS. THESE ROW POINTERS ARE STORED IN THE ARRAY IA; C I.E., IF M(I,J) IS THE FIRST NONZERO ENTRY (STORED) IN THE I-TH ROW C AND A(K) = M(I,J), THEN IA(I) = K. MOREOVER, IA(N+1) POINTS TO C THE FIRST LOCATION FOLLOWING THE LAST ELEMENT IN THE LAST ROW. C THUS, THE NUMBER OF ENTRIES IN THE I-TH ROW IS IA(I+1) - IA(I), C THE NONZERO ENTRIES IN THE I-TH ROW ARE STORED CONSECUTIVELY IN C C A(IA(I)), A(IA(I)+1), ..., A(IA(I+1)-1), C C AND THE CORRESPONDING COLUMN INDICES ARE STORED CONSECUTIVELY IN C C JA(IA(I)), JA(IA(I)+1), ..., JA(IA(I+1)-1). C C SINCE THE COEFFICIENT MATRIX IS SYMMETRIC, ONLY THE NONZERO ENTRIES C IN THE UPPER TRIANGLE NEED BE STORED, FOR EXAMPLE, THE MATRIX C C ( 1 0 2 3 0 ) C ( 0 4 0 0 0 ) C M = ( 2 0 5 6 0 ) C ( 3 0 6 7 8 ) C ( 0 0 0 8 9 ) C C COULD BE STORED AS C C \ 1 2 3 4 5 6 7 8 9 10 11 12 13 C ---+-------------------------------------- C IA \ 1 4 5 8 12 14 C JA \ 1 3 4 2 1 3 4 1 3 4 5 4 5 C A \ 1 2 3 4 2 5 6 3 6 7 8 8 9 C C OR (SYMMETRICALLY) AS C C \ 1 2 3 4 5 6 7 8 9 C ---+-------------------------- C IA \ 1 4 5 7 9 10 C JA \ 1 3 4 2 3 4 4 5 5 C A \ 1 2 3 4 5 6 7 8 9 . C C C REORDERING THE ROWS AND COLUMNS OF M C C A SYMMETRIC PERMUTATION OF THE ROWS AND COLUMNS OF THE COEFFICIENT C MATRIX M (E.G., WHICH REDUCES FILLIN OR ENHANCES NUMERICAL C STABILITY) MUST BE SPECIFIED. THE SOLUTION Z IS RETURNED IN THE C ORIGINAL ORDER. C C TO SPECIFY THE TRIVIAL ORDERING (I.E., THE IDENTITY PERMUTATION), C SET P(I) = IP(I) = I, I=1,...,N. IN THIS CASE, P AND IP CAN BE C THE SAME ARRAY. C C IF A NONTRIVIAL ORDERING (I.E., NOT THE IDENTITY PERMUTATION) IS C SPECIFIED AND M IS STORED SYMMETRICALLY (I.E., NOT BOTH M(I,J) AND C M(J,I) ARE STORED FOR I NE J), THEN ODRV SHOULD BE CALLED (WITH C PATH = 3 OR 5) TO SYMMETRICALLY REORDER (IA,JA,A) BEFORE CALLING C SDRV. THIS IS TO ENSURE THAT IF M(I,J) WILL BE IN THE UPPER C TRIANGLE OF M WITH RESPECT TO THE NEW ORDERING, THEN M(I,J) IS C STORED IN ROW I (AND THUS M(J,I) IS NOT STORED); WHEREAS IF M(I,J) C WILL BE IN THE STRICT LOWER TRIANGLE OF M, THEN M(J,I) IS STORED IN C ROW J (AND THUS M(I,J) IS NOT STORED). C C C PARAMETERS C C N - NUMBER OF VARIABLES/EQUATIONS C C P - INTEGER ONE-DIMENSIONAL ARRAY SPECIFYING A PERMUTATION OF C THE ROWS AND COLUMNS OF M; DIMENSION = N C C IP - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING THE INVERSE OF THE C PERMUTATION SPECIFIED IN P; I.E., IP(P(I)) = I, I=1,...,N; C DIMENSION = N C C IA - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING POINTERS TO DELIMIT C ROWS IN JA AND A; DIMENSION = N+1 C C JA - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING THE COLUMN INDICES C CORRESPONDING TO THE ELEMENTS OF A; DIMENSION = NUMBER OF C NONZERO ENTRIES IN M STORED C C A - REAL ONE-DIMENSIONAL ARRAY CONTAINING THE NONZERO ENTRIES IN C THE COEFFICIENT MATRIX M, STORED BY ROWS; DIMENSION = C NUMBER OF NONZERO ENTRIES IN M STORED C C B - REAL ONE-DIMENSIONAL ARRAY CONTAINING THE RIGHT-HAND SIDE B; C B AND Z CAN BE THE SAME ARRAY; DIMENSION = N C C Z - REAL ONE-DIMENSIONAL ARRAY CONTAINING THE SOLUTION X; Z AND C B CAN BE THE SAME ARRAY; DIMENSION = N C C NSP - DECLARED DIMENSION OF THE ONE-DIMENSIONAL ARRAYS ISP AND C RSP; NSP MUST BE (SUBSTANTIALLY) LARGER THAN 3N+2K, WHERE C K = NUMBER OF NONZERO ENTRIES IN THE UPPER TRIANGLE OF M C C ISP - INTEGER ONE-DIMENSIONAL ARRAY USED FOR WORKING STORAGE; ISP C AND RSP SHOULD BE EQUIVALENCED; DIMENSION = NSP C C RSP - REAL ONE-DIMENSIONAL ARRAY USED FOR WORKING STORAGE; RSP C AND ISP SHOULD BE EQUIVALENCED; DIMENSION = NSP C C ESP - INTEGER VARIABLE; IF SUFFICIENT STORAGE WAS AVAILABLE TO C PERFORM THE SYMBOLIC FACTORIZATION (SSF), THEN ESP IS SET TO C THE AMOUNT OF EXCESS STORAGE PROVIDED (NEGATIVE IF C INSUFFICIENT STORAGE WAS AVAILABLE TO PERFORM THE NUMERIC C FACTORIZATION (SNF)) C C PATH - INTEGER PATH SPECIFICATION; VALUES AND THEIR MEANINGS ARE - C 1 PERFORM SSF, SNF, AND SNS C 2 PERFORM SNF AND SNS (ISP/RSP IS ASSUMED TO HAVE BEEN C SET UP IN AN EARLIER CALL TO SDRV (FOR SSF)) C 3 PERFORM SNS ONLY (ISP/RSP IS ASSUMED TO HAVE BEEN SET C UP IN AN EARLIER CALL TO SDRV (FOR SSF AND SNF)) C 4 PERFORM SSF C 5 PERFORM SSF AND SNF C 6 PERFORM SNF ONLY (ISP/RSP IS ASSUMED TO HAVE BEEN SET C UP IN AN EARLIER CALL TO SDRV (FOR SSF)) C C FLAG - INTEGER ERROR FLAG; VALUES AND THEIR MEANINGS ARE - C C 0 NO ERRORS DETECTED C 2N+K DUPLICATE ENTRY IN A -- ROW = K C 6N+K INSUFFICIENT STORAGE IN SSF -- ROW = K C 7N+1 INSUFFICIENT STORAGE IN SNF C 8N+K ZERO PIVOT -- ROW = K C 10N+1 INSUFFICIENT STORAGE IN SDRV C 11N+1 ILLEGAL PATH SPECIFICATION C C ==================== change #4 (insertion) =======================4 C <0 MATRIX NOT SUFF. POS-DEF (detected in SNFMOD) C FLAG IS SET TO MINUS THE INDEX (IN THE ORIGINAL C MATRIX) OF THE LARGEST ADDITION in E C ==================================================================end C C C CONVERSION FROM REAL TO DOUBLE PRECISION C C CHANGE THE REAL DECLARATIONS IN SDRV, SNF, AND SNS TO DOUBLE C PRECISION DECLARATIONS; AND CHANGE THE VALUE IN THE DATA STATEMENT C FOR THE INTEGER VARIABLE RATIO (IN SDRV) FROM 1 TO 2. C C NOTE: FOR CRAY, SET RATIO to 1! C----------------------------------------------------------------------- C C REAL A(*), B(*), Z(*), RSP(*) C DATA RATIO/1/ C .. Scalar Arguments .. INTEGER ESP,FLAG,N,NSP,PATH C .. C .. Array Arguments .. DOUBLE PRECISION A(*),B(*),RSP(*),Z(*) INTEGER IA(*),IP(*),ISP(*),JA(*),P(*) C .. C .. Local Scalars .. INTEGER D,IJU,IL,IU,JL,JU,JUMAX,MARK,Q,RATIO,TMP,U,UMAX C .. C .. External Subroutines .. EXTERNAL SNFMOD,SNS,SSF C .. C .. Data statements .. DATA RATIO/2/ C .. C C----VALIDATE PATH SPECIFICATION IF (PATH.LT.1 .OR. 6.LT.PATH) GO TO 60 C C----ALLOCATE STORAGE AND FACTOR M SYMBOLICALLY TO DETERMINE FILL-IN IJU = 1 IU = IJU + N JL = IU + N + 1 JU = JL + N Q = (NSP+1) - N MARK = Q - N JUMAX = MARK - JU C IF ((PATH-1)* (PATH-4)* (PATH-5).NE.0) GO TO 10 IF (JUMAX.LE.0) GO TO 50 CALL SSF(N,P,IP,IA,JA,ISP(IJU),ISP(JU),ISP(IU),JUMAX,ISP(Q), + ISP(MARK),ISP(JL),FLAG) IF (FLAG.NE.0) GO TO 40 C C----ALLOCATE STORAGE AND FACTOR M NUMERICALLY 10 IL = JU + ISP(IJU+ (N-1)) TMP = ((IL-1)+ (RATIO-1))/RATIO + 1 D = TMP + N U = D + N UMAX = (NSP+1) - U ESP = UMAX - (ISP(IU+N)-1) C IF ((PATH-1)* (PATH-2)* (PATH-5)* (PATH-6).NE.0) GO TO 20 IF (UMAX.LE.0) GO TO 50 C ==================== change #5 (replacement) =====================5 C WAS: CALL SNF C ==================================================================end CALL SNFMOD(N,P,IP,IA,JA,A,RSP(D),ISP(IJU),ISP(JU),ISP(IU),RSP(U), + UMAX,ISP(IL),ISP(JL),FLAG) C ==================== change #6 (replacement) =====================6 C WAS: IF (FLAG.NE.0) GO TO 100 C ==================================================================end IF (FLAG.GT.0) GO TO 40 C C----SOLVE SYSTEM OF LINEAR EQUATIONS MX = B 20 IF ((PATH-1)* (PATH-2)* (PATH-3).NE.0) GO TO 30 IF (UMAX.LE.0) GO TO 50 CALL SNS(N,P,RSP(D),ISP(IJU),ISP(JU),ISP(IU),RSP(U),Z,B,RSP(TMP)) C 30 RETURN C C ** ERROR -- ERROR DETECTED IN SSF, SNF, OR SNS 40 RETURN C ** ERROR -- INSUFFICIENT STORAGE 50 FLAG = 10*N + 1 RETURN C ** ERROR -- ILLEGAL PATH SPECIFICATION 60 FLAG = 11*N + 1 RETURN END C C*********************************************************************** C*********************************************************************** C NUMERICAL FACTORIZATION OF SYMMETRIC MATRICES C*********************************************************************** C C ==================== change #1 (replacement) =====================1 C WAS: C SNF -- NUMERICAL UT-D-U FACTORIZATION OF SPARSE SYMMETRIC POSITIVE C DEFINITE MATRIX C SUBROUTINE SNF C ===================================================================== C C SNFMOD -- NUMERICAL FACTORIZATION OF SPARSE SYMMETRIC MATRICES M BY C THE GILL/MURRAY/WRIGHT MODIFIED CHOLESKY FACTORIZATION (GMW C MCF) WITHOUT PIVOTING. THE FACTORIZATION PRODUCES U,D, AND C E SO THAT M + E = UT-D-U, WHERE E AND D ARE DIAGONAL C MATRICES. THIS ROUTINE IS A MODIFICATION OF THE YSMP C routine SNF. ALL CHANGES ARE INDICATED. C C ==================================================================end SUBROUTINE SNFMOD(N,P,IP,IA,JA,A,D,IJU,JU,IU,U,UMAX,IL,JL,FLAG) C C ADDITIONAL PARAMETERS C C IL - INTEGER ONE-DIMENSIONAL WORK ARRAY; DIMENSION = N C C JL - INTEGER ONE-DIMENSIONAL WORK ARRAY; DIMENSION = N C C C DEFINITIONS OF INTERNAL PARAMETERS (DURING K-TH STAGE OF ELIMINATION) C C (D(I),I=K,N) CONTAINS THE K-TH ROW OF U (EXPANDED) C C IL(I) POINTS TO THE FIRST NONZERO ELEMENT IN COLUMNS K,...,N OF C ROW I OF U C C JL CONTAINS LISTS OF ROWS TO BE ADDED TO UNELIMINATED ROWS -- C I GE K => JL(I) IS THE FIRST ROW TO BE ADDED TO ROW I C I LT K => JL(I) IS THE ROW FOLLOWING ROW I IN SOME LIST OF ROWS C IN EITHER CASE, JL(I) = 0 INDICATES THE END OF A LIST C C----------------------------------------------------------------------- C .. Scalar Arguments .. INTEGER FLAG,N,UMAX C .. C .. Array Arguments .. DOUBLE PRECISION A(*),D(*),U(*) INTEGER IA(*),IJU(*),IL(*),IP(*),IU(*),JA(*),JL(*),JU(*),P(*) C .. C .. Scalars in Common .. DOUBLE PRECISION TAU INTEGER MC,SRLS C .. C .. Local Scalars .. DOUBLE PRECISION BOUND,DEL,DK,EK,ELT,ELT2,ELTNEW,EMAX,EPS,EPS1, + GAMMA,UKIDI,W,WW,XI,XIN,ZERO INTEGER I,ILI,IROW,J,JMAX,JMIN,JUMUJ,K,KK,KKMAX,KKMIN,MU,NEXTI,VJ C .. C .. Intrinsic Functions .. INTRINSIC ABS,MAX,MIN C .. C .. Common blocks .. COMMON /TAUINPUT/TAU,SRLS,MC C .. FLAG = 0 ZERO = 0.0D0 EMAX = ZERO C ==================================================================end C C----CHECK FOR SUFFICIENT STORAGE FOR U IF (IU(N+1)-1.GT.UMAX) GO TO 140 C C----INITIALIZATION DO 10 K = 1,N D(K) = 0 JL(K) = 0 10 CONTINUE C ==================== change #3 (insertion) =======================3 C Calculate GAMMA and XI, the largest magnitudes of the diag. and off- C diag. elements, respectively. When the diag. elts. are stored first C in A in each row (PATH = 4 or 5 in ODRV), GAMMA=max(GAMMA,A(IA(i))), C i=1,...,IA(n+1)-1. We assume that this IS the case. (If this were C later changed, then for each row I we would have to loop through KK C = KKMIN,..., KKMAX where KKMIN = IA(I), KKMAX = IA(I+1)-1, and test C whether I = JA(KK), ie. row index = column index. If this equality C holds, the element is a diagonal). Then calculate DEL and BOUND: C DEL = max ( max(XI,GAMMA)*EPS, EPS) where EPS is a small given C number, and BOUND = max ( XI/N, GAMMA, EPS). C ===================================================================== EPS = 1.0D-06 GAMMA = ZERO XI = ZERO DO 30 IROW = 1,N GAMMA = MAX(GAMMA,ABS(A(IA(IROW)))) KKMIN = IA(IROW) + 1 KKMAX = IA(IROW+1) - 1 IF (KKMIN.GT.KKMAX) GO TO 30 DO 20 KK = KKMIN,KKMAX XI = MAX(XI,ABS(A(KK))) 20 CONTINUE 30 CONTINUE EPS1 = MAX(GAMMA,XI)*EPS DEL = MAX(EPS,EPS1) XIN = N XIN = XI/XIN BOUND = MAX(GAMMA,XIN,EPS) C ==================================================================end C C----FOR EACH ROW K DO 130 K = 1,N C C------INITIALIZE K-TH ROW WITH ELEMENTS NONZERO IN ROW P(K) OF M JMIN = IA(P(K)) JMAX = IA(P(K)+1) - 1 IF (JMIN.GT.JMAX) GO TO 50 DO 40 J = JMIN,JMAX VJ = IP(JA(J)) IF (K.LE.VJ) D(VJ) = A(J) 40 CONTINUE C C------MODIFY K-TH ROW BY ADDING IN THOSE ROWS I WITH U(I,K) NE 0 C------FOR EACH ROW I TO BE ADDED IN 50 DK = D(K) I = JL(K) 60 IF (I.EQ.0) GO TO 90 NEXTI = JL(I) C C--------COMPUTE MULTIPLIER AND UPDATE DIAGONAL ELEMENT ILI = IL(I) UKIDI = -U(ILI)*D(I) DK = DK + UKIDI*U(ILI) U(ILI) = UKIDI C C--------ADD MULTIPLE OF ROW I TO K-TH ROW ... JMIN = ILI + 1 JMAX = IU(I+1) - 1 IF (JMIN.GT.JMAX) GO TO 80 MU = IJU(I) - IU(I) DO 70 J = JMIN,JMAX D(JU(MU+J)) = D(JU(MU+J)) + UKIDI*U(J) 70 CONTINUE C C--------... AND ADD I TO ROW LIST FOR NEXT NONZERO ENTRY IL(I) = JMIN J = JU(MU+JMIN) JL(I) = JL(J) JL(J) = I C 80 I = NEXTI GO TO 60 C C ==================== change #4 (replacement) =====================4 C WAS: C------CHECK FOR ZERO PIVOT C 9 IF (DK.EQ.0) GO TO 108 C ===================================================================== C STATEMENT 9 ABOVE WILL BE MODIFIED TO RESET Dk IN THE EVENT THE C THE MATRIX IS NOT SUFF. POSITIVE-DEFINITE. NOTE THAT EVEN WHEN Dk>0, C IT MAY BE MODIFIED IF THE MATRIX IS NOT POS. DEF! C C Dk is set as: Dk = MAX ( ABS(Dk), DEL, (ELT**2)/BOUND), where C ELT is the largest magnitude among the elements in the Kth row of U. C This restriction guarantees that all elts. of D are strictly positive C and that the elts. of the factors satisfy a uniform bound. C [ Recall that we work with the auxiliary quantities Vik = Uik * Dk. C The bound we want to impose on the elts. of U, C ( max(Uik)**2 ) * Dk <= BOUND, is equivalent to C ( max(Vik)**2 ) / Dk <= BOUND, or C Dk >= (max(Vik)**2) / BOUND.) C The value for ELT = max(Vik), max over i for fixed k, is found by C looping through the appropriate elements of U. These elements C are currently stored in part of D. ] C C ===================================================================== C 90 W = DK C ===================================================================== ELT = ZERO JMIN = IU(K) JMAX = IU(K+1) - 1 MU = IJU(K) - JMIN IF (JMIN.GT.JMAX) GO TO 110 DO 100 J = JMIN,JMAX ELTNEW = ABS(D(JU(MU+J))) ELT = MAX(ELT,ELTNEW) 100 CONTINUE 110 CONTINUE C ELT2 = ELT*ELT c--- xie: the modified Cholesky (MC) factorization c MC = 1: the original MC by Gill and Murray c MC = 2: the UMC by Xie & Schlick 6/12/96 c---------------------------------------------------- IF (MC.EQ.1) THEN DK = MAX(ABS(DK),DEL,ELT2/BOUND) EK = DK - W ELSE WW = DK + TAU IF (WW.GT.DEL) THEN DK = MAX(WW,ELT2/BOUND) ELSE IF (ABS(WW).LE.DEL) THEN DK = DEL ELSE DK = MIN(WW,-ELT2/BOUND) END IF EK = DK - W END IF c-------------------------- End of UMC IF (EK.GT.EMAX) THEN EMAX = EK FLAG = -P(K) END IF C------SAVE DIAGONAL ELEMENT D(K) = 1/DK C------SAVE NONZERO ENTRIES IN K-TH ROW OF U ... JMIN = IU(K) JMAX = IU(K+1) - 1 IF (JMIN.GT.JMAX) GO TO 130 MU = IJU(K) - JMIN DO 120 J = JMIN,JMAX JUMUJ = JU(MU+J) U(J) = D(JUMUJ) D(JUMUJ) = 0 120 CONTINUE C------ AND ADD K TO ROW LIST FOR FIRST NONZERO ENTRY IN K-TH ROW IL(K) = JMIN I = JU(MU+JMIN) JL(K) = JL(I) JL(I) = K 130 CONTINUE C uncomment next line to check for the largest diagonal modification C IF (FLAG .LT. 0) WRITE (6,*) ' NMAX, EMAX', FLAG, EMAX C C ==================== change #5 (deletion) ========================5 C WAS: FLAG = 0 C ==================================================================end RETURN C C ** ERROR -- INSUFFICIENT STORAGE FOR U 140 FLAG = 7*N + 1 RETURN END SHAR_EOF fi # end of overwriting check cd .. # End of shell archive exit 0