2 C Copyright 2002 Free Software Foundation, Inc.
4 C This file is part of GNU Radio
6 C GNU Radio is free software; you can redistribute it and/or modify
7 C it under the terms of the GNU General Public License as published by
8 C the Free Software Foundation; either version 3, or (at your option)
11 C GNU Radio is distributed in the hope that it will be useful,
12 C but WITHOUT ANY WARRANTY; without even the implied warranty of
13 C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 C GNU General Public License for more details.
16 C You should have received a copy of the GNU General Public License
17 C along with GNU Radio; see the file COPYING. If not, write to
18 C the Free Software Foundation, Inc., 51 Franklin Street,
19 C Boston, MA 02110-1301, USA.
21 DOUBLE PRECISION FUNCTION PRAX2(F,INITV,NDIM,OUT)
22 DOUBLE PRECISION INITV(128),OUT(128), F
26 DOUBLE PRECISION V,X,D,Q0,Q1,DMIN,EPSMCH,FX,H,QD0,QD1,QF1,
27 * SMALL,T,XLDT,XM2,XM4,DSEED,SCBD
29 COMMON /CPRAX/ V(128,128),X(128),D(128),Q0(128),Q1(128),
30 * DMIN,EPSMCH,FX,H,QD0,QD1,QF1,SMALL,T,XLDT,XM2,XM4,DSEED,SCBD,
31 * N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
41 C -1 produces no diagnostic output
59 C PRASET 1.0 JUNE 1995
61 C SET INITIAL VALUES FOR SOME QUANTITIES USED IN SUBROUTINE PRAXIS.
62 C THE USER CAN RESET THESE, IF DESIRED,
63 C AFTER CALLING PRASET AND BEFORE CALLING PRAXIS.
65 C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT,
66 C OKLAHOMA STATE UNIVERSITY
68 C ON MANY MACHINES, SUBROUTINE PRAXIS WILL CAUSE UNDERFLOW AND/OR
69 C DIVIDE CHECK WHEN COMPUTING EPSMCH**4 AND EPSMCH**(-4).
70 C IN THAT CASE, SET EPSMCH=1.0D-9 (OR POSSIBLY EPSMCH=1.0D-8)
71 C AFTER CALLING SUBROUTINE PRASET.
73 INTEGER N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
76 DOUBLE PRECISION V,X,D,Q0,Q1,DMIN,EPSMCH,FX,H,QD0,QD1,QF1,
77 * SMALL,T,XLDT,XM2,XM4,DSEED,SCBD
78 DOUBLE PRECISION A,B,XMID,XPLUS,RZERO,UNITR,RTWO
80 COMMON /CPRAX/ V(128,128),X(128),D(128),Q0(128),Q1(128),
81 * DMIN,EPSMCH,FX,H,QD0,QD1,QF1,SMALL,T,XLDT,XM2,XM4,DSEED,SCBD,
82 * N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
88 C NMAX IS THE DIMENSION OF THE ARRAYS V(*,*), X(*), D(*),
93 C NFMAX IS THE MAXIMUM NUMBER OF FUNCTION EVALUATIONS PERMITTED.
97 C LP IS THE LOGICAL UNIT NUMBER FOR PRINTED OUTPUT.
101 C T IS A CONVERGENCE TOLERANCE USED IN SUBROUTINE PRAXIS.
105 C JPRINT CONTROLS PRINTED OUTPUT IN PRAXIS.
109 C H IS AN ESTIMATE OF THE DISTANCE FROM THE INITIAL POINT
114 C USE BISECTION TO COMPUTE THE VALUE OF EPSMCH, "MACHINE EPSILON".
115 C EPSMCH IS THE SMALLEST FLOATING POINT (REAL OR DOUBLE PRECISION)
116 C NUMBER WHICH, WHEN ADDED TO ONE, GIVES A RESULT GREATER THAN ONE.
121 IF(XMID.LE.A .OR. XMID.GE.B) GO TO 20
123 IF(XPLUS.GT.UNITR) THEN
136 C JRANCH = 1 TO USE BRENT'S RANDOM,
137 C JRANCH = 2 TO USE FUNCTION DRANDM.
143 C DSEED IS AN INITIAL SEED FOR DRANDM,
144 C A SUBROUTINE THAT GENERATES PSEUDORANDOM NUMBERS
145 C UNIFORMLY DISTRIBUTED ON (0,1).
149 C SCBD IS AN UPPER BOUND ON THE SCALE FACTORS IN PRAXIS.
150 C IF THE AXES MAY BE BADLY SCALED (WHICH IS TO BE AVOIDED IF
151 C POSSIBLE) THEN SET SCBD = 10, OTHERWISE 1.
155 C ILLCIN IS THE INITIAL VALUE OF ILLC,
156 C THE FLAG THAT SIGNALS AN ILL-CONDITIONED PROBLEM.
157 C IF THE PROBLEM IS KNOWN TO BE ILL-CONDITIONED SET ILLCIN=1,
162 C KTM IS A CONVERGENCE SWITCH USED IN PRAXIS.
163 C KTM+1 IS THE NUMBER OF ITERATIONS WITHOUT IMPROVEMENT
164 C BEFORE THE ALGORITHM TERMINATES.
165 C KTM=4 IS VERY CAUTIOUS.
166 C USUALLY KTM=1 IS SATISFACTORY.
177 C PRAXIS 2.0 JUNE 1995
179 C THE PRAXIS PACKAGE MINIMIZES THE FUNCTION F(X,N) OF N
180 C VARIABLES X(1),...,X(N), USING THE PRINCIPAL AXIS METHOD.
181 C F MUST BE A SMOOTH (CONTINUOUSLY DIFFERENTIABLE) FUNCTION.
183 C "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
184 C RICHARD P. BRENT, PRENTICE-HALL 1973 (ISBN 0-13-022335-2),
187 C TRANSLATED FROM ALGOL W TO A.N.S.I. 1966 STANDARD BASIC FORTRAN
188 C BY ROSALEE TAYLOR AND SUE PINSKI, COMPUTER SCIENCE DEPARTMENT,
189 C OKLAHOMA STATE UNIVERSITY (DECEMBER 1973).
191 C UPDATED TO A.N.S.I. STANDARD FORTRAN 77 BY J. P. CHANDLER
192 C COMPUTER SCIENCE DEPARTMENT, OKLAHOMA STATE UNIVERSITY
195 C SUBROUTINE PRAXIS CALLS SUBPROGRAMS
196 C F, MINX, RANDOM (OR DRANDM), QUAD, MINFIT, SORT.
198 C SUBROUTINE QUAD CALLS MINX.
200 C SUBROUTINE MINX CALLS FLIN.
202 C SUBROUTINE FLIN CALLS F.
205 C INPUT QUANTITIES (SET IN THE CALLING PROGRAM)...
207 C F FUNCTION F(X,N) TO BE MINIMIZED
209 C X(*) INITIAL GUESS OF MINIMUM
211 C N DIMENSION OF X (NOTE... N MUST BE .GE. 2)
213 C H MAXIMUM STEP SIZE
217 C EPSMCH MACHINE PRECISION
219 C JPRINT PRINT SWITCH
222 C OUTPUT QUANTITIES...
224 C X(*) ESTIMATED POINT OF MINIMUM
228 C NL NUMBER OF LINEAR SEARCHES
230 C NF NUMBER OF FUNCTION EVALUATIONS
232 C V(*,*) EIGENVECTORS OF A
235 C D(*) EIGENVALUES OF A
241 C ON ENTRY X(*) HOLDS A GUESS. ON RETURN IT HOLDS THE ESTIMATED
242 C POINT OF MINIMUM, WITH (HOPEFULLY)
243 C ABS(ERROR) LESS THAN SQRT(EPSMCH)*ABS(X) + T, WHERE
244 C EPSMCH IS THE MACHINE PRECISION, THE SMALLEST NUMBER SUCH THAT
245 C (1 + EPSMCH) IS GREATER THAN 1.
249 C H IS THE MAXIMUM STEP SIZE, SET TO ABOUT THE MAXIMUM EXPECTED
250 C DISTANCE FROM THE GUESS TO THE MINIMUM. IF H IS SET TOO
251 C SMALL OR TOO LARGE THEN THE INITIAL RATE OF CONVERGENCE WILL
254 C THE USER SHOULD OBSERVE THE COMMENT ON HEURISTIC NUMBERS
255 C AT THE BEGINNING OF THE SUBROUTINE.
257 C JPRINT CONTROLS THE PRINTING OF INTERMEDIATE RESULTS.
258 C IT USES SUBROUTINES FLIN, MINX, QUAD, SORT, AND MINFIT.
259 C IF JPRINT = 1, F IS PRINTED AFTER EVERY N+1 OR N+2 LINEAR
260 C MINIMIZATIONS, AND FINAL X IS PRINTED, BUT INTERMEDIATE
261 C X ONLY IF N IS LESS THAN OR EQUAL TO 4.
262 C IF JPRINT = 2, EIGENVALUES OF A AND SCALE FACTORS ARE ALSO PRINTED.
263 C IF JPRINT = 3, F AND X ARE PRINTED AFTER EVERY FEW LINEAR
265 C IF JPRINT = 4, EIGENVECTORS ARE ALSO PRINTED.
266 C IF JPRINT = 5, ADDITIONAL DEBUGGING INFORMATION IS ALSO PRINTED.
268 C RANDOM RETURNS A RANDOM NUMBER UNIFORMLY DISTRIBUTED IN (0, 1).
270 C THIS SUBROUTINE IS MACHINE-INDEPENDENT, APART FROM THE
271 C SPECIFICATION OF EPSMCH. WE ASSUME THAT EPSMCH**(-4) DOES NOT
272 C OVERFLOW (IF IT DOES THEN EPSMCH MUST BE INCREASED), AND THAT ON
273 C FLOATING-POINT UNDERFLOW THE RESULT IS SET TO ZERO.
275 INTEGER N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
276 INTEGER ILLC,I,IK,IM,IMU,J,K,KL,KM1,KT,K2
278 DOUBLE PRECISION V,X,D,Q0,Q1,DMIN,EPSMCH,FX,H,QD0,QD1,QF1,
279 * SMALL,T,XLDT,XM2,XM4,DSEED,SCBD
280 DOUBLE PRECISION F, Y,Z,E, DABS,DSQRT,ZABS,ZSQRT,DRANDM,
281 * HUNDRD,HUNDTH,ONE,PT9,RHALF,TEN,TENTH,TWO,ZERO,
282 * DF,DLDFAC,DN,F1,XF,XL,T2,RANVAL,ARG,
283 * VLARGE,VSMALL,XLARGE,XLDS,FXVALU,F1VALU,S,SF,SL
287 DIMENSION Y(128),Z(128),E(128)
289 COMMON /CPRAX/ V(128,128),X(128),D(128),Q0(128),Q1(128),
290 * DMIN,EPSMCH,FX,H,QD0,QD1,QF1,SMALL,T,XLDT,XM2,XM4,DSEED,SCBD,
291 * N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
294 ZSQRT(ARG)=DSQRT(ARG)
308 C MACHINE DEPENDENT NUMBERS...
310 C ON MANY COMPUTERS, VSMALL WILL UNDERFLOW,
311 C AND COMPUTING XLARGE MAY CAUSE A DIVISION BY ZERO.
312 C IN THAT CASE, EPSMCH SHOULD BE SET EQUAL TO 1.0D-9
313 C (OR POSSIBLY 1.0D-8) BEFORE CALLING PRAXIS.
322 C HEURISTIC NUMBERS...
324 C IF THE AXES MAY BE BADLY SCALED (WHICH IS TO BE AVOIDED IF
325 C POSSIBLE) THEN SET SCBD = 10, OTHERWISE 1.
327 C IF THE PROBLEM IS KNOWN TO BE ILL-CONDITIONED SET ILLC = 1,
330 C KTM+1 IS THE NUMBER OF ITERATIONS WITHOUT IMPROVEMENT
331 C BEFORE THE ALGORITHM TERMINATES.
332 C KTM=4 IS VERY CAUTIOUS.
333 C USUALLY KTM=1 IS SATISFACTORY.
335 C BRENT RECOMMENDED THE FOLLOWING VALUES FOR MOST PROBLEMS...
341 C SCBD, ILLCIN, AND KTM ARE NOW IN COMMON.
342 C THEY ARE INITIALIZED IN SUBROUTINE PRASET,
343 C AND CAN BE RESET BY THE USER AFTER CALLING PRASET.
361 IF(H.LT.HUNDRD*T) H=HUNDRD*T
374 C Q0(*) AND Q1(*) ARE PREVIOUS X(*) POINTS,
375 C INITIALIZED IN PRAXIS, USED IN FLIN,
376 C AND CHANGED IN QUAD.
381 C Q0(*) WAS NOT INITIALIZED IN BRENT'S ALGOL PROCEDURE.
388 40 FORMAT(/' NL =',I10,5X,'NF =',I10/5X,'FX =',1PG15.7)
390 IF(N.LE.4 .OR. JPRINT.GT.2) THEN
391 WRITE(LP,50)(X(I),I=1,N)
392 50 FORMAT(/8X,'X'/(1X,1PG15.7,4G15.7))
403 C MINIMIZE ALONG THE FIRST DIRECTION.
405 IF(JPRINT.GE.5) WRITE(LP,70)D(1),S,FX
406 70 FORMAT(/' CALL NO. 1 TO MINX.'/
407 * 5X,'D(1) =',1PG15.7,5X,'S =',G15.7,5X,'FX =',G15.7)
410 CALL MINX(1,2,D(1),S,FXVALU,0,F)
418 IF(SF.LE.PT9*D(1) .OR. PT9*SF.GE.D(1)) THEN
445 C TAKE A RANDOM STEP TO GET OUT OF A RESOLUTION VALLEY.
447 C PRAXIS ASSUMES THAT RANDOM (OR DRANDM) RETURNS
448 C A PSEUDORANDOM NUMBER UNIFORMLY DISTRIBUTED IN (0,1),
449 C AND THAT ANY INITIALIZATION OF THE RANDOM NUMBER GENERATOR
450 C HAS ALREADY BEEN DONE.
460 S=(TENTH*XLDT+T2*TEN**KT)*(RANVAL-RHALF)
471 IF(JPRINT.GE.1) WRITE(LP,140)NF,SF,FX
472 140 FORMAT(/' ***** RANDOM STEP IN PRAXIS. NF =',I11/
473 * 5X,'SF =',1PG15.7,5X,'FX =',G15.7)
481 C MINIMIZE ALONG NON-CONJUGATE DIRECTIONS.
483 IF(JPRINT.GE.5) WRITE(LP,150)K2,D(K2),S,FX
484 150 FORMAT(/' CALL NO. 2 TO MINX.'/
485 * 5X,'K2 =',I4,5X,'D(K2) =',1PG15.7,5X,
486 * 'S =',G15.7/5X,'FX =',G15.7)
489 CALL MINX(K2,2,D(K2),S,FXVALU,0,F)
503 170 IF(ILLC.EQ.0 .AND. DF.LT.ZABS(HUNDRD*EPSMCH*FX)) THEN
505 C NO SUCCESS WITH ILLC=0, SO TRY ONCE WITH ILLC=1 .
514 IF(K.EQ.2 .AND. JPRINT.GT.1) THEN
515 WRITE(LP,180)(D(I),I=1,N)
516 180 FORMAT(/' NEW D'/(1X,1PG15.7,4G15.7))
520 IF(KM1.LT.1) GO TO 210
523 C MINIMIZE ALONG CONJUGATE DIRECTIONS.
525 IF(JPRINT.GE.5) WRITE(LP,190)K2,D(K2),S,FX
526 190 FORMAT(/' CALL NO. 3 TO MINX.'/
527 * 5X,'K2 =',I4,5X,'D(K2) =',1PG15.7,5X,
528 * 'S =',G15.7/5X,'FX =',G15.7)
532 CALL MINX(K2,2,D(K2),S,FXVALU,0,F)
548 IF(XLDS.GT.SMALL) THEN
550 C THROW AWAY THE DIRECTION KL AND MINIMIZE ALONG
551 C THE NEW CONJUGATE DIRECTION.
554 IF(K.GT.IK) GO TO 250
571 IF(JPRINT.GE.5) WRITE(LP,270)K,D(K),XLDS,F1
572 270 FORMAT(/' CALL NO. 4 TO MINX.'/
573 * 5X,'K =',I4,5X,'D(K) =',1PG15.7,5X,
574 * 'XLDS =',G15.7/5X,'F1 =',G15.7)
577 CALL MINX(K,4,D(K),XLDS,F1VALU,1,F)
579 IF(XLDS.LE.ZERO) THEN
589 IF(XLDT.LT.XLDS) XLDT=XLDS
593 IF(N.LE.4 .OR. JPRINT.GT.2) THEN
594 WRITE(LP,50)(X(I),I=1,N)
604 C SEE IF THE STEP LENGTH EXCEEDS HALF THE TOLERANCE.
606 IF(XLDT.GT.RHALF*T2) THEN
614 IF(KT.GT.KTM) GO TO 550
618 300 FORMAT(/' IN PRAXIS, NF REACHED THE LIMIT NFMAX =',I11/
619 * 5X,'THIS IS AN ABNORMAL TERMINATION.'/
620 * 5X,'THE FUNCTION HAS NOT BEEN MINIMIZED AND',
621 * ' THE RESULTING X(*) VECTOR SHOULD NOT BE USED.')
627 C TRY QUADRATIC EXTRAPOLATION IN CASE WE ARE STUCK IN A CURVED VALLEY.
634 IF(DN.LT.D(I)) DN=D(I)
640 340 FORMAT(/' NEW DIRECTIONS')
643 WRITE(LP,350)I,(V(I,J),J=1,N)
644 350 FORMAT(1X,I5,4X,1PG15.7,4G15.7/(10X,5G15.7))
658 C SCALE THE AXES TO TRY TO REDUCE THE CONDITION NUMBER.
669 IF(Z(I).LT.XM4) Z(I)=XM4
677 IF(Z(I).GT.SCBD) THEN
682 C IT APPEARS THAT THERE ARE TWO MISSING END; STATEMENTS
683 C AT THIS POINT IN BRENT'S LISTING.
688 C TRANSPOSE V FOR MINFIT.
701 C FIND THE SINGULAR VALUE DECOMPOSITION OF V.
702 C THIS GIVES THE EIGENVALUES AND PRINCIPAL AXES
703 C OF THE APPROXIMATING QUADRATIC FORM
704 C WITHOUT SQUARING THE CONDITION NUMBER.
706 440 CALL MINFIT(N,EPSMCH,VSMALL,V,D,E,NMAX,LP)
739 IF(DN*D(I).GT.XLARGE) THEN
741 ELSE IF(DN*D(I).LT.SMALL) THEN
744 D(I)=ONE/(DN*D(I))**2
748 C SORT THE NEW EIGENVALUES AND EIGENVECTORS.
753 IF(DMIN.LT.SMALL) DMIN=SMALL
755 IF(XM2*D(1).GT.DMIN) THEN
761 IF(JPRINT.GT.1 .AND. SCBD.GT.ONE) THEN
762 WRITE(LP,510)(Z(I),I=1,N)
763 510 FORMAT(/' SCALE FACTORS'/(1X,1PG15.7,4G15.7))
767 WRITE(LP,520)(D(I),I=1,N)
768 520 FORMAT(/' EIGENVALUES OF A'/(1X,1PG15.7,4G15.7))
774 530 FORMAT(/' EIGENVECTORS OF A')
777 WRITE(LP,350)I,(V(I,J),J=1,N)
781 C GO BACK TO THE MAIN LOOP.
784 C HANDLE THE CASE N .EQ. 1 IN AN AD HOC WAY.
785 C (BRENT DID NOT PROVIDE FOR THIS CASE.)
791 550 IF(JPRINT.GT.0) THEN
792 WRITE(LP,560)(X(I),I=1,N)
793 560 FORMAT(//7X,'X'/(1X,1PG15.7,4G15.7))
798 IF(JPRINT.GE.0) WRITE(LP,570)FX,NL,NF
799 570 FORMAT(/' EXIT PRAXIS. FX =',1PG25.17,5X,'NL =',I8,
809 C THIS SUBROUTINE LOOKS FOR THE MINIMUM ALONG
810 C A CURVE DEFINED BY Q0, Q1, AND X.
812 C "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
813 C RICHARD P. BRENT, PRENTICE-HALL 1973, PAGE 161
815 C SUBROUTINE QUAD IS CALLED BY SUBROUTINE PRAXIS.
817 INTEGER N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
820 DOUBLE PRECISION V,X,D,Q0,Q1,DMIN,EPSMCH,FX,H,QD0,QD1,QF1,
821 * SMALL,T,XLDT,XM2,XM4,DSEED,SCBD
822 DOUBLE PRECISION F, DSQRT,ZSQRT,ARG,
823 * ONE,QA,QB,QC,S,TWO,XL,ZERO,QF1VAL
827 COMMON /CPRAX/ V(128,128),X(128),D(128),Q0(128),Q1(128),
828 * DMIN,EPSMCH,FX,H,QD0,QD1,QF1,SMALL,T,XLDT,XM2,XM4,DSEED,SCBD,
829 * N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
831 ZSQRT(ARG)=DSQRT(ARG)
853 IF(QD0.GT.ZERO .AND. QD1.GT.ZERO .AND. NL.GE.3*N*N) THEN
855 IF(JPRINT.GE.1) WRITE(LP,20)NF,QD0,QD1,FX,QF1
856 20 FORMAT(/' ***** CALL MINX FROM QUAD. NF =',I11/
857 * 5X,'QD0 =',1PG15.7,5X,'QD1 =',G15.7/
858 * 5X,'FX =',G15.7,5X,'QF1 =',G15.7)
861 CALL MINX(0,2,S,XL,QF1VAL,1,F)
862 QA=XL*(XL-QD1)/(QD0*(QD0+QD1))
863 QB=(XL+QD0)*(QD1-XL)/(QD0*QD1)
864 QC=XL*(XL+QD0)/(QD1*(QD0+QD1))
877 X(I)=QA*S+QB*X(I)+QC*Q1(I)
885 SUBROUTINE MINX(J,NITS,D2,X1,F1,IFK,F)
887 C SUBROUTINE MINX MINIMIZES F FROM X IN THE DIRECTION V(*,J)
888 C UNLESS J IS LESS THAN 1, WHEN A QUADRATIC SEARCH IS DONE IN
889 C THE PLANE DEFINED BY Q0, Q1, AND X.
891 C "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
892 C RICHARD P. BRENT, PRENTICE-HALL 1973, PAGES 159-160
894 C SUBROUTINE MINX IS CALLED BY SUBROUTINES PRAXIS AND QUAD.
896 C D2 AND X1 RETURN RESULTS.
897 C J, NITS, F1 AND IFK ARE VALUE PARAMETERS THAT RETURN NOTHING.
898 C DO NOT SEND A COMMON VARIABLE TO MINX FOR PARAMETER F1.
901 C D2 IS AN APPROXIMATION TO HALF OF
902 C THE SECOND DERIVATIVE OF F (OR ZERO).
904 C X1 IS AN ESTIMATE OF DISTANCE TO MINIMUM,
905 C RETURNED AS THE DISTANCE FOUND.
907 C IF IFK = 1 THEN F1 IS FLIN(X1), OTHERWISE X1 AND F1 ARE
908 C IGNORED ON ENTRY UNLESS FINAL FX IS GREATER THAN F1.
910 C NITS CONTROLS THE NUMBER OF TIMES AN ATTEMPT IS MADE TO
911 C HALVE THE INTERVAL.
915 INTEGER N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
916 INTEGER IFK,J,NITS, I,IDZ,K
918 DOUBLE PRECISION V,X,D,Q0,Q1,DMIN,EPSMCH,FX,H,QD0,QD1,QF1,
919 * SMALL,T,XLDT,XM2,XM4,DSEED,SCBD
920 DOUBLE PRECISION D2,X1,
921 * DABS,DSQRT,ZABS,ZSQRT,ARG,
922 * HUNDTH,RHALF,TWO,ZERO,
923 * DENOM,D1,FM,F0,F1,F2,S,SF1,SX1,T2,XM,X2
925 COMMON /CPRAX/ V(128,128),X(128),D(128),Q0(128),Q1(128),
926 * DMIN,EPSMCH,FX,H,QD0,QD1,QF1,SMALL,T,XLDT,XM2,XM4,DSEED,SCBD,
927 * N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
929 ZSQRT(ARG)=DSQRT(ARG)
944 IF(D2.LT.EPSMCH) THEN
950 C FIND THE STEP SIZE.
964 T2=XM4*ZSQRT(ZABS(FX)/DENOM+S*XLDT)+XM2*XLDT
966 IF(IDZ.EQ.1 .AND. T2.GT.S) T2=S
967 IF(T2.LT.SMALL) T2=SMALL
968 IF(T2.GT.HUNDTH*H) T2=HUNDTH*H
970 IF(IFK.EQ.1 .AND. F1.LE.FM) THEN
975 IF(IFK.EQ.0 .OR. ZABS(X1).LT.T2) THEN
995 C EVALUATE FLIN AT ANOTHER POINT,
996 C AND ESTIMATE THE SECOND DERIVATIVE.
1004 CALL FLIN(X2,J,F,F2)
1011 D2=(X2*(F1-F0)-X1*(F2-F0))/(X1*X2*(X1-X2))
1013 IF(JPRINT.GE.5) WRITE(LP,30)X1,X2,F0,F1,F2,D2
1014 30 FORMAT(/' COMPUTE D2 IN SUBROUTINE MINX.'/
1015 * 5X,'X1 =',1PG15.7,5X,'X2 =',G15.7/
1016 * 5X,'F0 =',G15.7,5X,'F1 =',G15.7,5X,'F2 =',G15.7/
1020 C ESTIMATE THE FIRST DERIVATIVE AT 0.
1025 C PREDICT THE MINIMUM.
1027 IF(D2.LE.SMALL) THEN
1039 IF(ZABS(X2).GT.H) THEN
1048 C EVALUATE F AT THE PREDICTED MINIMUM.
1051 40 CALL FLIN(X2,J,F,F2)
1053 IF(K.LT.NITS .AND. F2.GT.F0) THEN
1055 C NO SUCCESS, SO TRY AGAIN.
1061 IF(F0.LT.F1 .AND. X1*X2.GT.ZERO) GO TO 20
1070 C INCREMENT THE ONE-DIMENSIONAL SEARCH COUNTER.
1080 C GET A NEW ESTIMATE OF THE SECOND DERIVATIVE.
1082 IF(ZABS(X2*(X2-X1)).GT.SMALL) THEN
1083 D2=(X2*(F1-F0)-X1*(FM-F0))/(X1*X2*(X1-X2))
1085 IF(JPRINT.GE.5) WRITE(LP,50)X1,X2,F0,FM,F1,D2
1086 50 FORMAT(/' RECOMPUTE D2 IN SUBROUTINE MINX.'/
1087 * 5X,'X1 =',1PG15.7,5X,'X2 =',G15.7/
1088 * 5X,'F0 =',G15.7,5X,'FM =',G15.7,5X,'F1 =',G15.7/
1091 ELSE IF(K.GT.0) THEN
1094 IF(JPRINT.GE.5) WRITE(LP,60)
1095 60 FORMAT(/' SET D2=0 IN SUBROUTINE MINX.')
1100 IF(D2.LE.SMALL) THEN
1103 IF(JPRINT.GE.5) WRITE(LP,70)D2
1104 70 FORMAT(/' SET D2=SMALL=',1PG15.7,' IN SUBROUTINE MINX.')
1107 IF(JPRINT.GE.5) WRITE(LP,80)X1,X2,FX,FM,SF1
1108 80 FORMAT(/' SUBROUTINE MINX. X1 =',1PG15.7,5X,'X2 =',G15.7/
1109 * 5X,'FX =',G15.7,5X,'FM =',G15.7,5X,'SF1 =',G15.7)
1118 C UPDATE X FOR A LINEAR SEARCH BUT NOT FOR A PARABOLIC SEARCH.
1127 IF(JPRINT.GE.5) WRITE(LP,100)D2,X1,F1,FX
1128 100 FORMAT(/' LEAVE SUBROUTINE MINX.'/
1129 * 5X,'D2 =',1PG15.7,5X,'X1 =',G15.7,5X,'F1 =',G15.7/
1137 SUBROUTINE FLIN(XL,J,F,FLN)
1139 C FLIN IS A FUNCTION OF ONE VARIABLE XL WHICH IS MINIMIZED BY
1142 C "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
1143 C RICHARD P. BRENT, PRENTICE-HALL 1973, PAGES 159-160
1145 C SUBROUTINE FLIN IS CALLED BY SUBROUTINE MINX.
1147 INTEGER N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
1150 DOUBLE PRECISION V,X,D,Q0,Q1,DMIN,EPSMCH,FX,H,QD0,QD1,QF1,
1151 * SMALL,T,XLDT,XM2,XM4,DSEED,SCBD
1152 DOUBLE PRECISION XL,F,FLN, TT, QA,QB,QC
1156 COMMON /CPRAX/ V(128,128),X(128),D(128),Q0(128),Q1(128),
1157 * DMIN,EPSMCH,FX,H,QD0,QD1,QF1,SMALL,T,XLDT,XM2,XM4,DSEED,SCBD,
1158 * N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
1165 TT(I)=X(I)+XL*V(I,J)
1170 C SEARCH ALONG A PARABOLIC SPACE CURVE.
1172 QA=XL*(XL-QD1)/(QD0*(QD0+QD1))
1173 QB=(XL+QD0)*(QD1-XL)/(QD0*QD1)
1174 QC=XL*(XL+QD0)/(QD1*(QD0+QD1))
1177 TT(I)=QA*Q0(I)+QB*X(I)+QC*Q1(I)
1181 C INCREMENT FUNCTION EVALUATION COUNTER.
1191 SUBROUTINE MINFIT(N,EPS,TOL,AB,Q,E,NMAX,LP)
1193 C AN IMPROVED VERSION OF MINFIT, RESTRICTED TO M=N, P=0.
1194 C SEE GOLUB AND REINSCH (1970).
1196 C "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
1197 C RICHARD P. BRENT, PRENTICE-HALL 1973, PAGES 156-158
1199 C G. H. GOLUB AND C. REINSCH,
1200 C "SINGULAR VALUE DECOMPOSITION AND LEAST SQUARES SOLUTIONS',
1201 C NUMERISCHE MATHEMATIK 14 (1970) PAGES 403-420
1203 C THE SINGULAR VALUES OF THE ARRAY AB ARE RETURNED IN Q,
1204 C AND AB IS OVERWRITTEN WITH THE ORTHOGONAL MATRIX V SUCH THAT
1205 C U.DIAG(Q)=AB.V, WHERE U IS ANOTHER ORTHOGONAL MATRIX.
1207 C SUBROUTINE MINFIT IS CALLED BY SUBROUTINE PRAXIS.
1210 * I,II,J,JTHIRT,K,KK,KT,L,LL2,LPI,L2
1212 DOUBLE PRECISION EPS,TOL,AB,Q,E,
1213 * DABS,DSQRT,ZABS,ZSQRT,ARG,
1214 * C,DENOM,F,G,H,ONE,X,Y,Z,ZERO,S,TWO
1216 DIMENSION AB(NMAX,N),Q(N),E(N)
1219 ZSQRT(ARG)=DSQRT(ARG)
1227 C HOUSEHOLDER'S REDUCTION TO BIDIAGONAL FORM...
1267 AB(K,J)=AB(K,J)+F*AB(K,I)
1283 80 IF(S.LT.TOL) THEN
1296 IF(L.GT.N) GO TO 130
1309 AB(J,K)=AB(J,K)+S*E(K)
1314 130 Y=ZABS(Q(I))+ZABS(E(I))
1319 C ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS...
1327 IF(L.GT.N) GO TO 200
1340 AB(K,J)=AB(K,J)+S*AB(K,I)
1345 IF(L.GT.N) GO TO 200
1356 C DIAGONALIZATION OF THE BIDIAGONAL FORM...
1363 C LABEL TESTFSPLITTING...
1367 IF(KT.GT.JTHIRT) THEN
1370 230 FORMAT(' QR FAILED.')
1377 C IF(...) GO TO TESTFCONVERGENCE
1379 IF(ZABS(E(L)).LE.EPS) GO TO 270
1381 C IF(...) GO TO CANCELLATION
1383 IF(ZABS(Q(L-1)).LE.EPS) GO TO 250
1386 C CANCELLATION OF E(L) IF L IS GREATER THAN 1...
1387 C LABEL CANCELLATION...
1391 IF(L.GT.K) GO TO 270
1396 C IF(...) GO TO TESTFCONVERGENCE
1398 IF(ZABS(F).LE.EPS) GO TO 270
1401 IF(ZABS(F).LT.ZABS(G)) THEN
1402 H=ZABS(G)*ZSQRT(ONE+(F/G)**2)
1403 ELSE IF(F.NE.ZERO) THEN
1404 H=ZABS(F)*ZSQRT(ONE+(G/F)**2)
1416 C THE ABOVE REPLACES Q(I) AND H BY SQUARE ROOT OF (G*G+F*F)
1417 C WHICH MAY GIVE INCORRECT RESULTS IF THE SQUARES UNDERFLOW OR IF
1424 C LABEL TESTFCONVERGENCE...
1428 C IF(...) GO TO CONVERGENCE
1430 IF(L.EQ.K) GO TO 310
1432 C SHIFT FROM BOTTOM 2*2 MINOR.
1438 F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(TWO*H*Y)
1447 F=((X-Z)*(X+Z)+H*(Y/DENOM-H))/X
1449 C NEXT QR TRANSFORMATION...
1454 IF(LPI.GT.K) GO TO 300
1461 IF(ZABS(F).LT.ZABS(H)) THEN
1462 Z=ZABS(H)*ZSQRT(ONE+(F/H)**2)
1463 ELSE IF(F.NE.ZERO) THEN
1464 Z=ZABS(F)*ZSQRT(ONE+(H/F)**2)
1490 IF(ZABS(F).LT.ZABS(H)) THEN
1491 Z=ZABS(H)*ZSQRT(ONE+(F/H)**2)
1492 ELSE IF(F.NE.ZERO) THEN
1493 Z=ZABS(F)*ZSQRT(ONE+(H/F)**2)
1515 C GO TO TESTFSPLITTING
1519 C LABEL CONVERGENCE...
1521 310 IF(Z.LT.ZERO) THEN
1523 C Q(K) IS MADE NON-NEGATIVE.
1539 C THIS SUBROUTINE SORTS THE ELEMENTS OF D
1540 C AND THE CORRESPONDING COLUMNS OF V INTO DESCENDING ORDER.
1542 C "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
1543 C RICHARD P. BRENT, PRENTICE-HALL 1973, PAGES 158-159
1545 INTEGER N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
1546 INTEGER I,IPI,J,K,NMI
1548 DOUBLE PRECISION V,X,D,Q0,Q1,DMIN,EPSMCH,FX,H,QD0,QD1,QF1,
1549 * SMALL,T,XLDT,XM2,XM4,DSEED,SCBD
1552 COMMON /CPRAX/ V(128,128),X(128),D(128),Q0(128),Q1(128),
1553 * DMIN,EPSMCH,FX,H,QD0,QD1,QF1,SMALL,T,XLDT,XM2,XM4,DSEED,SCBD,
1554 * N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
1557 IF(NMI.LT.1) GO TO 50
1562 IF(IPI.GT.N) GO TO 20
1589 SUBROUTINE RANINI(RVALUE)
1591 C SUBROUTINE RANINI PERFORMS INITIALIZATION FOR SUBROUTINE RANDOM.
1593 C "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
1594 C RICHARD P. BRENT, PRENTICE-HALL 1973, PAGES 163-164
1598 DOUBLE PRECISION RVALUE,R,RAN3,DMOD,DABS,RAN1
1600 COMMON /COMRAN/ RAN3(127),RAN1,JRAN2
1602 R=DMOD(DABS(RVALUE),8190.0D0)+1
1605 10 IF(JRAN2.GT.0) THEN
1610 R=DMOD(1756.0D0*R,8191.0D0)
1611 RAN1=(RAN1+(R-DMOD(R,32.0D0))/32.0D0)/256.0D0
1623 SUBROUTINE RANDOM(RANVAL)
1625 C SUBROUTINE RANDOM RETURNS A DOUBLE PRECISION PSEUDORANDOM NUMBER
1626 C UNIFORMLY DISTRIBUTED IN (0,1) (INCLUDING 0 BUT NOT 1).
1628 C "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
1629 C RICHARD P. BRENT, PRENTICE-HALL 1973, PAGES 163-164
1631 C BEFORE THE FIRST CALL TO RANDOM, THE USER MUST
1632 C CALL RANINI(R) ONCE (ONLY) WITH R A DOUBLE PRECISION NUMBER
1633 C EQUAL TO ANY INTEGER VALUE.
1634 C BRENT (PAGE 166) USED THE EQUIVALENT OF
1635 C CALL RANINI(4.0D0) .
1637 C THE ALGORITHM USED IN SUBROUTINE RANDOM RETURNS X(N)/2**56,
1638 C WHERE X(N) = X(N-1) + X(N-127) (MOD 2**56) .
1639 C SINCE (1 + X + X**127) IS PRIMITIVE (MOD 2),
1640 C THE PERIOD IS AT LEAST (2**127 - 1), WHICH EXCEEDS 10**38.
1642 C SEE "SEMINUMERICAL ALGORITHMS", VOLUME 2 OF
1643 C "THE ART OF COMPUTER PROGRAMMING" BY DONALD E. KNUTH,
1644 C ADDISON-WESLEY 1969, PAGES 26, 34, AND 464.
1646 C X(N) IS STORED IN DOUBLE PRECISION AS RAN3 = X(N)/2**56 - 1/2,
1647 C AND ALL DOUBLE PRECISION ARITHMETIC IS EXACT.
1651 DOUBLE PRECISION RANVAL,RAN3,RAN1
1653 COMMON /COMRAN/ RAN3(127),RAN1,JRAN2
1661 RAN1=RAN1+RAN3(JRAN2+1)
1662 IF(RAN1.LT.0.0D0) THEN
1676 DOUBLE PRECISION FUNCTION DRANDM(DL)
1678 C SIMPLE PORTABLE PSEUDORANDOM NUMBER GENERATOR.
1680 C DRANDM RETURNS FUNCTION VALUES THAT ARE PSEUDORANDOM
1681 C NUMBERS UNIFORMLY DISTRIBUTED ON THE INTERVAL (0,1).
1683 C 'NUMERICAL MATHEMATICS AND COMPUTING' BY WARD CHENEY AND
1684 C DAVID KINCAID, BROOKS/COLE PUBLISHING COMPANY
1685 C (FIRST EDITION, 1980), PAGE 203
1687 C AT THE BEGINNING OF EXECUTION, OR WHENEVER A NEW SEQUENCE IS
1688 C TO BE INITIATED, SET DL EQUAL TO AN INTEGER VALUE BETWEEN
1689 C 1.0D0 AND 2147483646.0D0, INCLUSIVE. DO THIS ONLY ONCE.
1690 C THEREAFTER, DO NOT SET OR ALTER DL IN ANY WAY.
1691 C FUNCTION DRANDM WILL MODIFY DL FOR ITS OWN PURPOSES.
1693 C DRANDM USES A MULTIPLICATIVE CONGRUENTIAL METHOD.
1694 C THE NUMBERS GENERATED BY DRANDM SUFFER FROM THE PARALLEL
1695 C PLANES DEFECT DISCOVERED BY G. MARSAGLIA, AND SHOULD NOT BE
1696 C USED WHEN HIGH-QUALITY RANDOMNESS IS REQUIRED. IN THAT
1697 C CASE, USE A "SHUFFLING" METHOD.
1699 DOUBLE PRECISION DL,DMOD
1701 10 DL=DMOD(16807.0D0*DL,2147483647.0D0)
1702 DRANDM=DL/2147483647.0D0
1703 IF(DRANDM.LE.0.0D0 .OR. DRANDM.GE.1.0D0) GO TO 10