1C
2C Copyright 2002 Free Software Foundation, Inc.
3C
4C This file is part of GNU Radio
5C
6C GNU Radio is free software; you can redistribute it and/or modify
7C it under the terms of the GNU General Public License as published by
8C the Free Software Foundation; either version 3, or (at your option)
9C any later version.
10C
11C GNU Radio is distributed in the hope that it will be useful,
12C but WITHOUT ANY WARRANTY; without even the implied warranty of
13C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14C GNU General Public License for more details.
15C
16C You should have received a copy of the GNU General Public License
17C along with GNU Radio; see the file COPYING.  If not, write to
18C the Free Software Foundation, Inc., 51 Franklin Street,
19C Boston, MA 02110-1301, USA.
20C
21      DOUBLE PRECISION FUNCTION PRAX2(F,INITV,NDIM,OUT)
22      DOUBLE PRECISION INITV(128),OUT(128), F
23      INTEGER NDIM
24      EXTERNAL F
25C
26      DOUBLE PRECISION V,X,D,Q0,Q1,DMIN,EPSMCH,FX,H,QD0,QD1,QF1,
27     *   SMALL,T,XLDT,XM2,XM4,DSEED,SCBD
28C
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
32
33C
34      N=NDIM
35      do 10 I=1,N
36 10      X(I) = INITV(I)
37
38C
39      call praset
40
41C -1 produces no diagnostic output
42      jprint = -1
43      nfmax = 3000
44C tighter tolerance
45      T=1.0D-6
46C
47      call praxis(f)
48C
49      do 30 I=1,N
50 30      OUT(I) = X(I)
51C
52      prax2 = fx
53      return
54      end
55
56
57      SUBROUTINE PRASET
58C
59C  PRASET 1.0           JUNE 1995
60C
61C  SET INITIAL VALUES FOR SOME QUANTITIES USED IN SUBROUTINE PRAXIS.
62C  THE USER CAN RESET THESE, IF DESIRED,
63C  AFTER CALLING PRASET AND BEFORE CALLING PRAXIS.
64C
65C  J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT,
66C     OKLAHOMA STATE UNIVERSITY
67C
68C  ON MANY MACHINES, SUBROUTINE PRAXIS WILL CAUSE UNDERFLOW AND/OR
69C  DIVIDE CHECK WHEN COMPUTING EPSMCH**4 AND EPSMCH**(-4).
70C  IN THAT CASE, SET EPSMCH=1.0D-9 (OR POSSIBLY EPSMCH=1.0D-8)
71C  AFTER CALLING SUBROUTINE PRASET.
72C
73      INTEGER N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
74      INTEGER J
75C
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
79C
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
83C
84      RZERO=0.0D0
85      UNITR=1.0D0
86      RTWO=2.0D0
87C
88C  NMAX IS THE DIMENSION OF THE ARRAYS V(*,*), X(*), D(*),
89C  Q0(*), AND Q1(*).
90C
91      NMAX=128
92C
93C  NFMAX IS THE MAXIMUM NUMBER OF FUNCTION EVALUATIONS PERMITTED.
94C
95      NFMAX=100000
96C
97C  LP IS THE LOGICAL UNIT NUMBER FOR PRINTED OUTPUT.
98C
99      LP=6
100C
101C  T IS A CONVERGENCE TOLERANCE USED IN SUBROUTINE PRAXIS.
102C
103      T=1.0D-5
104C
105C  JPRINT CONTROLS PRINTED OUTPUT IN PRAXIS.
106C
107      JPRINT=4
108C
109C  H IS AN ESTIMATE OF THE DISTANCE FROM THE INITIAL POINT
110C  TO THE SOLUTION.
111C
112      H=1.0D0
113C
114C  USE BISECTION TO COMPUTE THE VALUE OF EPSMCH, "MACHINE EPSILON".
115C  EPSMCH IS THE SMALLEST FLOATING POINT (REAL OR DOUBLE PRECISION)
116C  NUMBER WHICH, WHEN ADDED TO ONE, GIVES A RESULT GREATER THAN ONE.
117C
118      A=RZERO
119      B=UNITR
120   10 XMID=A+(B-A)/RTWO
121      IF(XMID.LE.A .OR. XMID.GE.B) GO TO 20
122      XPLUS=UNITR+XMID
123      IF(XPLUS.GT.UNITR) THEN
124         B=XMID
125      ELSE
126         A=XMID
127      ENDIF
128      GO TO 10
129C
130   20 EPSMCH=B
131C
132      DO 30 J=1,NMAX
133         X(J)=RZERO
134   30 CONTINUE
135C
136C  JRANCH = 1 TO USE BRENT'S RANDOM,
137C  JRANCH = 2 TO USE FUNCTION DRANDM.
138C
139      JRANCH=1
140C
141      CALL RANINI(4.0D0)
142C
143C  DSEED IS AN INITIAL SEED FOR DRANDM,
144C  A SUBROUTINE THAT GENERATES PSEUDORANDOM NUMBERS
145C  UNIFORMLY DISTRIBUTED ON (0,1).
146C
147      DSEED=1234567.0D0
148C
149C  SCBD IS AN UPPER BOUND ON THE SCALE FACTORS IN PRAXIS.
150C  IF THE AXES MAY BE BADLY SCALED (WHICH IS TO BE AVOIDED IF
151C  POSSIBLE) THEN SET SCBD = 10, OTHERWISE 1.
152C
153      SCBD=1.0D0
154C
155C  ILLCIN IS THE INITIAL VALUE OF ILLC,
156C  THE FLAG THAT SIGNALS AN ILL-CONDITIONED PROBLEM.
157C  IF THE PROBLEM IS KNOWN TO BE ILL-CONDITIONED SET ILLCIN=1,
158C  OTHERWISE 0.
159C
160      ILLCIN=0
161C
162C  KTM IS A CONVERGENCE SWITCH USED IN PRAXIS.
163C  KTM+1 IS THE NUMBER OF ITERATIONS WITHOUT IMPROVEMENT
164C  BEFORE THE ALGORITHM TERMINATES.
165C  KTM=4 IS VERY CAUTIOUS.
166C  USUALLY KTM=1 IS SATISFACTORY.
167C
168      KTM=1
169C
170      RETURN
171C
172C  END PRASET
173C
174      END
175      SUBROUTINE PRAXIS(F)
176C
177C  PRAXIS 2.0        JUNE 1995
178C
179C  THE PRAXIS PACKAGE MINIMIZES THE FUNCTION F(X,N) OF N
180C  VARIABLES X(1),...,X(N), USING THE PRINCIPAL AXIS METHOD.
181C  F MUST BE A SMOOTH (CONTINUOUSLY DIFFERENTIABLE) FUNCTION.
182C
183C  "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
184C     RICHARD P. BRENT, PRENTICE-HALL 1973 (ISBN 0-13-022335-2),
185C     PAGES 156-167
186C
187C  TRANSLATED FROM ALGOL W TO A.N.S.I. 1966 STANDARD BASIC FORTRAN
188C     BY ROSALEE TAYLOR AND SUE PINSKI, COMPUTER SCIENCE DEPARTMENT,
189C     OKLAHOMA STATE UNIVERSITY (DECEMBER 1973).
190C
191C  UPDATED TO A.N.S.I. STANDARD FORTRAN 77 BY J. P. CHANDLER
192C     COMPUTER SCIENCE DEPARTMENT, OKLAHOMA STATE UNIVERSITY
193C
194C
195C  SUBROUTINE PRAXIS CALLS SUBPROGRAMS
196C     F, MINX, RANDOM (OR DRANDM), QUAD, MINFIT, SORT.
197C
198C  SUBROUTINE QUAD CALLS MINX.
199C
200C  SUBROUTINE MINX CALLS FLIN.
201C
202C  SUBROUTINE FLIN CALLS F.
203C
204C
205C  INPUT QUANTITIES (SET IN THE CALLING PROGRAM)...
206C
207C     F        FUNCTION F(X,N) TO BE MINIMIZED
208C
209C     X(*)     INITIAL GUESS OF MINIMUM
210C
211C     N        DIMENSION OF X  (NOTE...  N MUST BE .GE. 2)
212C
213C     H        MAXIMUM STEP SIZE
214C
215C     T        TOLERANCE
216C
217C     EPSMCH   MACHINE PRECISION
218C
219C     JPRINT   PRINT SWITCH
220C
221C
222C  OUTPUT QUANTITIES...
223C
224C     X(*)     ESTIMATED POINT OF MINIMUM
225C
226C     FX       VALUE OF F AT X
227C
228C     NL       NUMBER OF LINEAR SEARCHES
229C
230C     NF       NUMBER OF FUNCTION EVALUATIONS
231C
232C     V(*,*)   EIGENVECTORS OF A
233C                 NEW DIRECTIONS
234C
235C     D(*)     EIGENVALUES OF A
236C                 NEW D
237C
238C     Z(*)     SCALE FACTORS
239C
240C
241C  ON ENTRY X(*) HOLDS A GUESS.  ON RETURN IT HOLDS THE ESTIMATED
242C  POINT OF MINIMUM, WITH (HOPEFULLY)
243C  ABS(ERROR) LESS THAN SQRT(EPSMCH)*ABS(X) + T, WHERE
244C  EPSMCH IS THE MACHINE PRECISION, THE SMALLEST NUMBER SUCH THAT
245C  (1 + EPSMCH) IS GREATER THAN 1.
246C
247C  T IS A TOLERANCE.
248C
249C  H IS THE MAXIMUM STEP SIZE, SET TO ABOUT THE MAXIMUM EXPECTED
250C  DISTANCE FROM THE GUESS TO THE MINIMUM.  IF H IS SET TOO
251C  SMALL OR TOO LARGE THEN THE INITIAL RATE OF CONVERGENCE WILL
252C  BE SLOW.
253C
254C  THE USER SHOULD OBSERVE THE COMMENT ON HEURISTIC NUMBERS
255C  AT THE BEGINNING OF THE SUBROUTINE.
256C
257C  JPRINT CONTROLS THE PRINTING OF INTERMEDIATE RESULTS.
258C  IT USES SUBROUTINES FLIN, MINX, QUAD, SORT, AND MINFIT.
259C  IF JPRINT = 1, F IS PRINTED AFTER EVERY N+1 OR N+2 LINEAR
260C  MINIMIZATIONS, AND FINAL X IS PRINTED, BUT INTERMEDIATE
261C  X ONLY IF N IS LESS THAN OR EQUAL TO 4.
262C  IF JPRINT = 2, EIGENVALUES OF A AND SCALE FACTORS ARE ALSO PRINTED.
263C  IF JPRINT = 3, F AND X ARE PRINTED AFTER EVERY FEW LINEAR
264C  MINIMIZATIONS.
265C  IF JPRINT = 4, EIGENVECTORS ARE ALSO PRINTED.
266C  IF JPRINT = 5, ADDITIONAL DEBUGGING INFORMATION IS ALSO PRINTED.
267C
268C  RANDOM RETURNS A RANDOM NUMBER UNIFORMLY DISTRIBUTED IN (0, 1).
269C
270C  THIS SUBROUTINE IS MACHINE-INDEPENDENT, APART FROM THE
271C  SPECIFICATION OF EPSMCH.  WE ASSUME THAT EPSMCH**(-4) DOES NOT
272C  OVERFLOW (IF IT DOES THEN EPSMCH MUST BE INCREASED), AND THAT ON
273C  FLOATING-POINT UNDERFLOW THE RESULT IS SET TO ZERO.
274C
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
277C
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
284C
285      EXTERNAL F
286C
287      DIMENSION Y(128),Z(128),E(128)
288C
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
292C
293      ZABS(ARG)=DABS(ARG)
294      ZSQRT(ARG)=DSQRT(ARG)
295C
296C  INITIALIZATION...
297C
298      RHALF=0.5D0
299      ONE=1.0D0
300      TENTH=0.1D0
301      HUNDTH=0.01D0
302      HUNDRD=100.0D0
303      ZERO=0.0D0
304      PT9=0.9D0
305      TEN=10.0D0
306      TWO=2.0D0
307C
308C  MACHINE DEPENDENT NUMBERS...
309C
310C  ON MANY COMPUTERS, VSMALL WILL UNDERFLOW,
311C  AND COMPUTING XLARGE MAY CAUSE A DIVISION BY ZERO.
312C  IN THAT CASE, EPSMCH SHOULD BE SET EQUAL TO 1.0D-9
313C  (OR POSSIBLY 1.0D-8) BEFORE CALLING PRAXIS.
314C
315      SMALL=EPSMCH*EPSMCH
316      VSMALL=SMALL*SMALL
317      XLARGE=ONE/SMALL
318      VLARGE=ONE/VSMALL
319      XM2=ZSQRT(EPSMCH)
320      XM4=ZSQRT(XM2)
321C
322C  HEURISTIC NUMBERS...
323C
324C  IF THE AXES MAY BE BADLY SCALED (WHICH IS TO BE AVOIDED IF
325C  POSSIBLE) THEN SET SCBD = 10, OTHERWISE 1.
326C
327C  IF THE PROBLEM IS KNOWN TO BE ILL-CONDITIONED SET ILLC = 1,
328C  OTHERWISE 0.
329C
330C  KTM+1 IS THE NUMBER OF ITERATIONS WITHOUT IMPROVEMENT
331C  BEFORE THE ALGORITHM TERMINATES.
332C  KTM=4 IS VERY CAUTIOUS.
333C  USUALLY KTM=1 IS SATISFACTORY.
334C
335C  BRENT RECOMMENDED THE FOLLOWING VALUES FOR MOST PROBLEMS...
336C
337C                    SCBD=1.0
338C                    ILLC=0
339C                    KTM=1
340C
341C  SCBD, ILLCIN, AND KTM ARE NOW IN COMMON.
342C  THEY ARE INITIALIZED IN SUBROUTINE PRASET,
343C  AND CAN BE RESET BY THE USER AFTER CALLING PRASET.
344C
345      ILLC=ILLCIN
346C
347      IF(ILLC.EQ.1) THEN
348         DLDFAC=TENTH
349      ELSE
350         DLDFAC=HUNDTH
351      ENDIF
352C
353      KT=0
354      NL=0
355      NF=1
356      FX=F(X,N)
357      QF1=FX
358      T=SMALL+ZABS(T)
359      T2=T
360      DMIN=SMALL
361      IF(H.LT.HUNDRD*T) H=HUNDRD*T
362      XLDT=H
363C
364      DO 20 I=1,N
365         DO 10 J=1,N
366            V(I,J)=ZERO
367   10    CONTINUE
368         V(I,I)=ONE
369   20 CONTINUE
370C
371      QD0=ZERO
372      D(1)=ZERO
373C
374C  Q0(*) AND Q1(*) ARE PREVIOUS X(*) POINTS,
375C  INITIALIZED IN PRAXIS, USED IN FLIN,
376C  AND CHANGED IN QUAD.
377C
378      DO 30 I=1,N
379         Q1(I)=X(I)
380C
381C  Q0(*) WAS NOT INITIALIZED IN BRENT'S ALGOL PROCEDURE.
382C
383         Q0(I)=X(I)
384   30 CONTINUE
385C
386      IF(JPRINT.GT.0) THEN
387         WRITE(LP,40)NL,NF,FX
388   40    FORMAT(/' NL =',I10,5X,'NF =',I10/5X,'FX =',1PG15.7)
389C
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))
393         ENDIF
394      ENDIF
395C
396C  MAIN LOOP...
397C  LABEL L0...
398C
399   60 SF=D(1)
400      S=ZERO
401      D(1)=ZERO
402C
403C  MINIMIZE ALONG THE FIRST DIRECTION.
404C
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)
408C
409      FXVALU=FX
410      CALL MINX(1,2,D(1),S,FXVALU,0,F)
411C
412      IF(S.LE.ZERO) THEN
413         DO 80 I=1,N
414            V(I,1)=-V(I,1)
415   80    CONTINUE
416      ENDIF
417C
418      IF(SF.LE.PT9*D(1) .OR. PT9*SF.GE.D(1)) THEN
419C
420         IF(N.GE.2) THEN
421            DO 90 I=2,N
422               D(I)=ZERO
423   90       CONTINUE
424         ENDIF
425C
426      ENDIF
427C
428      IF(N.LT.2) GO TO 320
429      DO 310 K=2,N
430C
431         DO 100 I=1,N
432            Y(I)=X(I)
433  100    CONTINUE
434C
435         SF=FX
436         IF(KT.GT.0) ILLC=1
437C
438C  LABEL L1...
439C
440  110    KL=K
441         DF=ZERO
442C
443         IF(ILLC.EQ.1) THEN
444C
445C  TAKE A RANDOM STEP TO GET OUT OF A RESOLUTION VALLEY.
446C
447C  PRAXIS ASSUMES THAT RANDOM (OR DRANDM) RETURNS
448C  A PSEUDORANDOM NUMBER UNIFORMLY DISTRIBUTED IN (0,1),
449C  AND THAT ANY INITIALIZATION OF THE RANDOM NUMBER GENERATOR
450C  HAS ALREADY BEEN DONE.
451C
452            DO 130 I=1,N
453C
454               IF(JRANCH.EQ.1) THEN
455                  CALL RANDOM(RANVAL)
456               ELSE
457                  RANVAL=DRANDM(DSEED)
458               ENDIF
459C
460               S=(TENTH*XLDT+T2*TEN**KT)*(RANVAL-RHALF)
461               Z(I)=S
462C
463               DO 120 J=1,N
464                  X(J)=X(J)+S*V(J,I)
465  120          CONTINUE
466  130       CONTINUE
467C
468            FX=F(X,N)
469            NF=NF+1
470C
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)
474         ENDIF
475C
476         IF(K.GT.N) GO TO 170
477         DO 160 K2=K,N
478            SL=FX
479            S=ZERO
480C
481C  MINIMIZE ALONG NON-CONJUGATE DIRECTIONS.
482C
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)
487C
488            FXVALU=FX
489            CALL MINX(K2,2,D(K2),S,FXVALU,0,F)
490C
491            IF(ILLC.EQ.1) THEN
492               S=D(K2)*(S+Z(K2))**2
493            ELSE
494               S=SL-FX
495            ENDIF
496C
497            IF(DF.LT.S) THEN
498               DF=S
499               KL=K2
500            ENDIF
501  160    CONTINUE
502C
503  170    IF(ILLC.EQ.0 .AND. DF.LT.ZABS(HUNDRD*EPSMCH*FX)) THEN
504C
505C  NO SUCCESS WITH ILLC=0, SO TRY ONCE WITH ILLC=1 .
506C
507            ILLC=1
508C
509C  GO TO L1.
510C
511            GO TO 110
512         ENDIF
513C
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))
517         ENDIF
518C
519         KM1=K-1
520         IF(KM1.LT.1) GO TO 210
521         DO 200 K2=1,KM1
522C
523C  MINIMIZE ALONG CONJUGATE DIRECTIONS.
524C
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)
529C
530            S=ZERO
531            FXVALU=FX
532            CALL MINX(K2,2,D(K2),S,FXVALU,0,F)
533  200    CONTINUE
534C
535  210    F1=FX
536         FX=SF
537C
538         XLDS=ZERO
539         DO 220 I=1,N
540            SL=X(I)
541            X(I)=Y(I)
542            SL=SL-Y(I)
543            Y(I)=SL
544            XLDS=XLDS+SL*SL
545  220    CONTINUE
546C
547         XLDS=ZSQRT(XLDS)
548         IF(XLDS.GT.SMALL) THEN
549C
550C  THROW AWAY THE DIRECTION KL AND MINIMIZE ALONG
551C  THE NEW CONJUGATE DIRECTION.
552C
553            IK=KL-1
554            IF(K.GT.IK) GO TO 250
555            DO 240 IM=K,IK
556               I=IK-IM+K
557C
558               DO 230 J=1,N
559                  V(J,I+1)=V(J,I)
560  230          CONTINUE
561C
562               D(I+1)=D(I)
563  240       CONTINUE
564C
565  250       D(K)=ZERO
566C
567            DO 260 I=1,N
568               V(I,K)=Y(I)/XLDS
569  260       CONTINUE
570C
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)
575C
576            F1VALU=F1
577            CALL MINX(K,4,D(K),XLDS,F1VALU,1,F)
578C
579            IF(XLDS.LE.ZERO) THEN
580               XLDS=-XLDS
581C
582               DO 280 I=1,N
583                  V(I,K)=-V(I,K)
584  280          CONTINUE
585            ENDIF
586         ENDIF
587C
588         XLDT=DLDFAC*XLDT
589         IF(XLDT.LT.XLDS) XLDT=XLDS
590C
591         IF(JPRINT.GT.0) THEN
592            WRITE(LP,40)NL,NF,FX
593            IF(N.LE.4 .OR. JPRINT.GT.2) THEN
594               WRITE(LP,50)(X(I),I=1,N)
595            ENDIF
596         ENDIF
597C
598         T2=ZERO
599         DO 290 I=1,N
600            T2=T2+X(I)**2
601  290    CONTINUE
602         T2=XM2*ZSQRT(T2)+T
603C
604C  SEE IF THE STEP LENGTH EXCEEDS HALF THE TOLERANCE.
605C
606         IF(XLDT.GT.RHALF*T2) THEN
607            KT=0
608         ELSE
609            KT=KT+1
610         ENDIF
611C
612C  IF(...) GO TO L2
613C
614         IF(KT.GT.KTM) GO TO 550
615C
616         IF(NF.GE.NFMAX) THEN
617            WRITE(LP,300)NFMAX
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.')
622            GO TO 550
623         ENDIF
624C
625  310 CONTINUE
626C
627C  TRY QUADRATIC EXTRAPOLATION IN CASE WE ARE STUCK IN A CURVED VALLEY.
628C
629  320 CALL QUAD(F)
630C
631      DN=ZERO
632      DO 330 I=1,N
633         D(I)=ONE/ZSQRT(D(I))
634         IF(DN.LT.D(I)) DN=D(I)
635  330 CONTINUE
636C
637      IF(JPRINT.GT.3) THEN
638C
639         WRITE(LP,340)
640  340    FORMAT(/' NEW DIRECTIONS')
641C
642         DO 360 I=1,N
643            WRITE(LP,350)I,(V(I,J),J=1,N)
644  350       FORMAT(1X,I5,4X,1PG15.7,4G15.7/(10X,5G15.7))
645  360    CONTINUE
646      ENDIF
647C
648      DO 380 J=1,N
649C
650         S=D(J)/DN
651         DO 370 I=1,N
652            V(I,J)=S*V(I,J)
653  370    CONTINUE
654  380 CONTINUE
655C
656      IF(SCBD.GT.ONE) THEN
657C
658C  SCALE THE AXES TO TRY TO REDUCE THE CONDITION NUMBER.
659C
660         S=VLARGE
661         DO 400 I=1,N
662C
663            SL=ZERO
664            DO 390 J=1,N
665               SL=SL+V(I,J)**2
666  390       CONTINUE
667C
668            Z(I)=ZSQRT(SL)
669            IF(Z(I).LT.XM4) Z(I)=XM4
670            IF(S.GT.Z(I)) S=Z(I)
671  400    CONTINUE
672C
673         DO 410 I=1,N
674            SL=S/Z(I)
675            Z(I)=ONE/SL
676C
677            IF(Z(I).GT.SCBD) THEN
678               SL=ONE/SCBD
679               Z(I)=SCBD
680            ENDIF
681C
682C  IT APPEARS THAT THERE ARE TWO MISSING END; STATEMENTS
683C  AT THIS POINT IN BRENT'S LISTING.
684C
685  410    CONTINUE
686      ENDIF
687C
688C  TRANSPOSE V FOR MINFIT.
689C
690      IF(N.LT.2) GO TO 440
691      DO 430 I=2,N
692C
693         IMU=I-1
694         DO 420 J=1,IMU
695            S=V(I,J)
696            V(I,J)=V(J,I)
697            V(J,I)=S
698  420    CONTINUE
699  430 CONTINUE
700C
701C  FIND THE SINGULAR VALUE DECOMPOSITION OF V.
702C  THIS GIVES THE EIGENVALUES AND PRINCIPAL AXES
703C  OF THE APPROXIMATING QUADRATIC FORM
704C  WITHOUT SQUARING THE CONDITION NUMBER.
705C
706  440 CALL MINFIT(N,EPSMCH,VSMALL,V,D,E,NMAX,LP)
707C
708      IF(SCBD.GT.ONE) THEN
709C
710C  UNSCALING...
711C
712         DO 460 I=1,N
713C
714            S=Z(I)
715            DO 450 J=1,N
716               V(I,J)=S*V(I,J)
717  450       CONTINUE
718  460    CONTINUE
719C
720         DO 490 I=1,N
721C
722            S=ZERO
723            DO 470 J=1,N
724               S=S+V(J,I)**2
725  470       CONTINUE
726            S=ZSQRT(S)
727C
728            D(I)=S*D(I)
729C
730            S=ONE/S
731            DO 480 J=1,N
732               V(J,I)=S*V(J,I)
733  480       CONTINUE
734  490    CONTINUE
735      ENDIF
736C
737      DO 500 I=1,N
738C
739         IF(DN*D(I).GT.XLARGE) THEN
740            D(I)=VSMALL
741         ELSE IF(DN*D(I).LT.SMALL) THEN
742            D(I)=VLARGE
743         ELSE
744            D(I)=ONE/(DN*D(I))**2
745         ENDIF
746  500 CONTINUE
747C
748C  SORT THE NEW EIGENVALUES AND EIGENVECTORS.
749C
750      CALL SORT
751C
752      DMIN=D(N)
753      IF(DMIN.LT.SMALL) DMIN=SMALL
754C
755      IF(XM2*D(1).GT.DMIN) THEN
756         ILLC=1
757      ELSE
758         ILLC=0
759      ENDIF
760C
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))
764      ENDIF
765C
766      IF(JPRINT.GT.1) THEN
767         WRITE(LP,520)(D(I),I=1,N)
768  520    FORMAT(/' EIGENVALUES OF A'/(1X,1PG15.7,4G15.7))
769      ENDIF
770C
771      IF(JPRINT.GT.3) THEN
772C
773         WRITE(LP,530)
774  530    FORMAT(/' EIGENVECTORS OF A')
775C
776         DO 540 I=1,N
777            WRITE(LP,350)I,(V(I,J),J=1,N)
778  540    CONTINUE
779      ENDIF
780C
781C  GO BACK TO THE MAIN LOOP.
782C  GO TO L0
783C
784C  HANDLE THE CASE N .EQ. 1 IN AN AD HOC WAY.
785C  (BRENT DID NOT PROVIDE FOR THIS CASE.)
786C
787      IF(N.GE.2) GO TO 60
788C
789C  LABEL L2...
790C
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))
794      ENDIF
795C
796      FX=F(X,N)
797C
798      IF(JPRINT.GE.0) WRITE(LP,570)FX,NL,NF
799  570 FORMAT(/' EXIT PRAXIS.   FX =',1PG25.17,5X,'NL =',I8,
800     *   5X,'NF =',I9)
801C
802      RETURN
803C
804C  END PRAXIS
805C
806      END
807      SUBROUTINE QUAD(F)
808C
809C  THIS SUBROUTINE LOOKS FOR THE MINIMUM ALONG
810C  A CURVE DEFINED BY Q0, Q1, AND X.
811C
812C  "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
813C  RICHARD P. BRENT, PRENTICE-HALL 1973, PAGE 161
814C
815C  SUBROUTINE QUAD IS CALLED BY SUBROUTINE PRAXIS.
816C
817      INTEGER N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
818      INTEGER I
819C
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
824C
825      EXTERNAL F
826C
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
830C
831      ZSQRT(ARG)=DSQRT(ARG)
832C
833      ZERO=0.0D0
834      ONE=1.0D0
835C
836      S=FX
837      FX=QF1
838      QF1=S
839      QD1=ZERO
840C
841      DO 10 I=1,N
842         S=X(I)
843         XL=Q1(I)
844         X(I)=XL
845         Q1(I)=S
846         QD1=QD1+(S-XL)**2
847   10 CONTINUE
848C
849      QD1=ZSQRT(QD1)
850      XL=QD1
851      S=ZERO
852C
853      IF(QD0.GT.ZERO .AND. QD1.GT.ZERO .AND. NL.GE.3*N*N) THEN
854C
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)
859C
860         QF1VAL=QF1
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))
865      ELSE
866         FX=QF1
867         QA=ZERO
868         QB=ZERO
869         QC=ONE
870      ENDIF
871C
872      QD0=QD1
873C
874      DO 30 I=1,N
875         S=Q0(I)
876         Q0(I)=X(I)
877         X(I)=QA*S+QB*X(I)+QC*Q1(I)
878   30 CONTINUE
879C
880      RETURN
881C
882C  END QUAD
883C
884      END
885      SUBROUTINE MINX(J,NITS,D2,X1,F1,IFK,F)
886C
887C  SUBROUTINE MINX MINIMIZES F FROM X IN THE DIRECTION V(*,J)
888C  UNLESS J IS LESS THAN 1, WHEN A QUADRATIC SEARCH IS DONE IN
889C  THE PLANE DEFINED BY Q0, Q1, AND X.
890C
891C  "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
892C  RICHARD P. BRENT, PRENTICE-HALL 1973, PAGES 159-160
893C
894C  SUBROUTINE MINX IS CALLED BY SUBROUTINES PRAXIS AND QUAD.
895C
896C  D2 AND X1 RETURN RESULTS.
897C  J, NITS, F1 AND IFK ARE VALUE PARAMETERS THAT RETURN NOTHING.
898C  DO NOT SEND A COMMON VARIABLE TO MINX FOR PARAMETER F1.
899C
900C
901C  D2 IS AN APPROXIMATION TO HALF OF
902C  THE SECOND DERIVATIVE OF F (OR ZERO).
903C
904C  X1 IS AN ESTIMATE OF DISTANCE TO MINIMUM,
905C  RETURNED AS THE DISTANCE FOUND.
906C
907C  IF IFK = 1 THEN F1 IS FLIN(X1), OTHERWISE X1 AND F1 ARE
908C  IGNORED ON ENTRY UNLESS FINAL FX IS GREATER THAN F1.
909C
910C  NITS CONTROLS THE NUMBER OF TIMES AN ATTEMPT IS MADE TO
911C  HALVE THE INTERVAL.
912C
913      EXTERNAL F
914C
915      INTEGER N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
916      INTEGER IFK,J,NITS,     I,IDZ,K
917C
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
924C
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
928C
929      ZSQRT(ARG)=DSQRT(ARG)
930      ZABS(ARG)=DABS(ARG)
931C
932      HUNDTH=0.01D0
933      ZERO=0.0D0
934      TWO=2.0D0
935      RHALF=0.5D0
936C
937      SF1=F1
938      SX1=X1
939      K=0
940      XM=ZERO
941      FM=FX
942      F0=FX
943C
944      IF(D2.LT.EPSMCH) THEN
945         IDZ=1
946      ELSE
947         IDZ=0
948      ENDIF
949C
950C  FIND THE STEP SIZE.
951C
952      S=ZERO
953      DO 10 I=1,N
954         S=S+X(I)**2
955   10 CONTINUE
956      S=ZSQRT(S)
957C
958      IF(IDZ.EQ.1) THEN
959         DENOM=DMIN
960      ELSE
961         DENOM=D2
962      ENDIF
963C
964      T2=XM4*ZSQRT(ZABS(FX)/DENOM+S*XLDT)+XM2*XLDT
965      S=XM4*S+T
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
969C
970      IF(IFK.EQ.1 .AND. F1.LE.FM) THEN
971         XM=X1
972         FM=F1
973      ENDIF
974C
975      IF(IFK.EQ.0 .OR. ZABS(X1).LT.T2) THEN
976C
977         IF(X1.GE.ZERO) THEN
978            X1=T2
979         ELSE
980            X1=-T2
981         ENDIF
982C
983         CALL FLIN(X1,J,F,F1)
984      ENDIF
985C
986      IF(F1.LT.FM) THEN
987         XM=X1
988         FM=F1
989      ENDIF
990C
991C  LABEL L0...
992C
993   20 IF(IDZ.EQ.1) THEN
994C
995C  EVALUATE FLIN AT ANOTHER POINT,
996C  AND ESTIMATE THE SECOND DERIVATIVE.
997C
998         IF(F0.LT.F1) THEN
999            X2=-X1
1000         ELSE
1001            X2=TWO*X1
1002         ENDIF
1003C
1004         CALL FLIN(X2,J,F,F2)
1005C
1006         IF(F2.LE.FM) THEN
1007            XM=X2
1008            FM=F2
1009         ENDIF
1010C
1011         D2=(X2*(F1-F0)-X1*(F2-F0))/(X1*X2*(X1-X2))
1012C
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/
1017     *      5X,'D2 =',G15.7)
1018      ENDIF
1019C
1020C  ESTIMATE THE FIRST DERIVATIVE AT 0.
1021C
1022      D1=(F1-F0)/X1-X1*D2
1023      IDZ=1
1024C
1025C  PREDICT THE MINIMUM.
1026C
1027      IF(D2.LE.SMALL) THEN
1028C
1029         IF(D1.LT.ZERO) THEN
1030            X2=H
1031         ELSE
1032            X2=-H
1033         ENDIF
1034C
1035      ELSE
1036         X2=-RHALF*D1/D2
1037      ENDIF
1038C
1039      IF(ZABS(X2).GT.H) THEN
1040C
1041         IF(X2.GT.ZERO) THEN
1042            X2=H
1043         ELSE
1044            X2=-H
1045         ENDIF
1046      ENDIF
1047C
1048C  EVALUATE F AT THE PREDICTED MINIMUM.
1049C  LABEL L1...
1050C
1051   40 CALL FLIN(X2,J,F,F2)
1052C
1053      IF(K.LT.NITS .AND. F2.GT.F0) THEN
1054C
1055C  NO SUCCESS, SO TRY AGAIN.
1056C
1057         K=K+1
1058C
1059C  IF(...) GO TO L0
1060C
1061         IF(F0.LT.F1 .AND. X1*X2.GT.ZERO) GO TO 20
1062         X2=X2/TWO
1063C
1064C  GO TO L1
1065C
1066         GO TO 40
1067C
1068      ENDIF
1069C
1070C  INCREMENT THE ONE-DIMENSIONAL SEARCH COUNTER.
1071C
1072      NL=NL+1
1073C
1074      IF(F2.GT.FM) THEN
1075         X2=XM
1076      ELSE
1077         FM=F2
1078      ENDIF
1079C
1080C  GET A NEW ESTIMATE OF THE SECOND DERIVATIVE.
1081C
1082      IF(ZABS(X2*(X2-X1)).GT.SMALL) THEN
1083         D2=(X2*(F1-F0)-X1*(FM-F0))/(X1*X2*(X1-X2))
1084C
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/
1089     *      5X,'D2 =',G15.7)
1090C
1091      ELSE IF(K.GT.0) THEN
1092         D2=ZERO
1093C
1094         IF(JPRINT.GE.5) WRITE(LP,60)
1095   60    FORMAT(/' SET D2=0 IN SUBROUTINE MINX.')
1096      ELSE
1097         D2=D2
1098      ENDIF
1099C
1100      IF(D2.LE.SMALL) THEN
1101         D2=SMALL
1102C
1103         IF(JPRINT.GE.5) WRITE(LP,70)D2
1104   70    FORMAT(/' SET D2=SMALL=',1PG15.7,' IN SUBROUTINE MINX.')
1105      ENDIF
1106C
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)
1110C
1111      X1=X2
1112      FX=FM
1113      IF(SF1.LT.FX) THEN
1114         FX=SF1
1115         X1=SX1
1116      ENDIF
1117C
1118C  UPDATE X FOR A LINEAR SEARCH BUT NOT FOR A PARABOLIC SEARCH.
1119C
1120      IF(J.GT.0) THEN
1121C
1122         DO 90 I=1,N
1123            X(I)=X(I)+X1*V(I,J)
1124   90    CONTINUE
1125      ENDIF
1126C
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/
1130     *   5X,'FX =',G15.7)
1131C
1132      RETURN
1133C
1134C  END MINX
1135C
1136      END
1137      SUBROUTINE FLIN(XL,J,F,FLN)
1138C
1139C  FLIN IS A FUNCTION OF ONE VARIABLE XL WHICH IS MINIMIZED BY
1140C  SUBROUTINE MINX.
1141C
1142C  "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
1143C  RICHARD P. BRENT, PRENTICE-HALL 1973, PAGES 159-160
1144C
1145C  SUBROUTINE FLIN IS CALLED BY SUBROUTINE MINX.
1146C
1147      INTEGER N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
1148      INTEGER J,   I
1149C
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
1153C
1154      DIMENSION TT(128)
1155C
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
1159C
1160      IF(J.GT.0) THEN
1161C
1162C  LINEAR SEARCH...
1163C
1164         DO 10 I=1,N
1165            TT(I)=X(I)+XL*V(I,J)
1166   10    CONTINUE
1167C
1168      ELSE
1169C
1170C  SEARCH ALONG A PARABOLIC SPACE CURVE.
1171C
1172         QA=XL*(XL-QD1)/(QD0*(QD0+QD1))
1173         QB=(XL+QD0)*(QD1-XL)/(QD0*QD1)
1174         QC=XL*(XL+QD0)/(QD1*(QD0+QD1))
1175C
1176         DO 20 I=1,N
1177            TT(I)=QA*Q0(I)+QB*X(I)+QC*Q1(I)
1178   20    CONTINUE
1179      ENDIF
1180C
1181C  INCREMENT FUNCTION EVALUATION COUNTER.
1182C
1183      NF=NF+1
1184      FLN=F(TT,N)
1185C
1186      RETURN
1187C
1188C  END FLIN
1189C
1190      END
1191      SUBROUTINE MINFIT(N,EPS,TOL,AB,Q,E,NMAX,LP)
1192C
1193C  AN IMPROVED VERSION OF MINFIT, RESTRICTED TO M=N, P=0.
1194C  SEE GOLUB AND REINSCH (1970).
1195C
1196C  "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
1197C  RICHARD P. BRENT, PRENTICE-HALL 1973, PAGES 156-158
1198C
1199C  G. H. GOLUB AND C. REINSCH,
1200C  "SINGULAR VALUE DECOMPOSITION AND LEAST SQUARES SOLUTIONS',
1201C  NUMERISCHE MATHEMATIK 14 (1970) PAGES 403-420
1202C
1203C  THE SINGULAR VALUES OF THE ARRAY AB ARE RETURNED IN Q,
1204C  AND AB IS OVERWRITTEN WITH THE ORTHOGONAL MATRIX V SUCH THAT
1205C  U.DIAG(Q)=AB.V, WHERE U IS ANOTHER ORTHOGONAL MATRIX.
1206C
1207C  SUBROUTINE MINFIT IS CALLED BY SUBROUTINE PRAXIS.
1208C
1209      INTEGER N,NMAX,LP,
1210     *   I,II,J,JTHIRT,K,KK,KT,L,LL2,LPI,L2
1211C
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
1215C
1216      DIMENSION AB(NMAX,N),Q(N),E(N)
1217C
1218      ZABS(ARG)=DABS(ARG)
1219      ZSQRT(ARG)=DSQRT(ARG)
1220C
1221      JTHIRT=30
1222C
1223      ZERO=0.0D0
1224      ONE=1.0D0
1225      TWO=2.0D0
1226C
1227C  HOUSEHOLDER'S REDUCTION TO BIDIAGONAL FORM...
1228C
1229      X=ZERO
1230      G=ZERO
1231C
1232      DO 140 I=1,N
1233         E(I)=G
1234         S=ZERO
1235         L=I+1
1236C
1237         DO 10 J=I,N
1238            S=S+AB(J,I)**2
1239   10    CONTINUE
1240C
1241         IF(S.LT.TOL) THEN
1242            G=ZERO
1243         ELSE
1244            F=AB(I,I)
1245C
1246            IF(F.LT.ZERO) THEN
1247               G=ZSQRT(S)
1248            ELSE
1249               G=-ZSQRT(S)
1250            ENDIF
1251C
1252            H=F*G-S
1253            AB(I,I)=F-G
1254C
1255            IF(L.GT.N) GO TO 60
1256            DO 50 J=L,N
1257C
1258               F=ZERO
1259               IF(I.GT.N) GO TO 30
1260               DO 20 K=I,N
1261                  F=F+AB(K,I)*AB(K,J)
1262   20          CONTINUE
1263   30          F=F/H
1264C
1265               IF(I.GT.N) GO TO 50
1266               DO 40 K=I,N
1267                  AB(K,J)=AB(K,J)+F*AB(K,I)
1268   40          CONTINUE
1269   50       CONTINUE
1270         ENDIF
1271C
1272   60    Q(I)=G
1273         S=ZERO
1274C
1275         IF(I.LE.N) THEN
1276C
1277            IF(L.GT.N) GO TO 80
1278            DO 70 J=L,N
1279               S=S+AB(I,J)**2
1280   70       CONTINUE
1281         ENDIF
1282C
1283   80    IF(S.LT.TOL) THEN
1284            G=ZERO
1285         ELSE
1286            F=AB(I,I+1)
1287C
1288            IF(F.LT.ZERO) THEN
1289               G=ZSQRT(S)
1290            ELSE
1291               G=-ZSQRT(S)
1292            ENDIF
1293C
1294            H=F*G-S
1295            AB(I,I+1)=F-G
1296            IF(L.GT.N) GO TO 130
1297            DO 90 J=L,N
1298               E(J)=AB(I,J)/H
1299   90       CONTINUE
1300C
1301            DO 120 J=L,N
1302C
1303               S=ZERO
1304               DO 100 K=L,N
1305                  S=S+AB(J,K)*AB(I,K)
1306  100          CONTINUE
1307C
1308               DO 110 K=L,N
1309                  AB(J,K)=AB(J,K)+S*E(K)
1310  110          CONTINUE
1311  120       CONTINUE
1312         ENDIF
1313C
1314  130    Y=ZABS(Q(I))+ZABS(E(I))
1315C
1316         IF(Y.GT.X) X=Y
1317  140 CONTINUE
1318C
1319C  ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS...
1320C
1321      DO 210 II=1,N
1322         I=N-II+1
1323C
1324         IF(G.NE.ZERO) THEN
1325            H=AB(I,I+1)*G
1326C
1327            IF(L.GT.N) GO TO 200
1328            DO 150 J=L,N
1329               AB(J,I)=AB(I,J)/H
1330  150       CONTINUE
1331C
1332            DO 180 J=L,N
1333C
1334               S=ZERO
1335               DO 160 K=L,N
1336                  S=S+AB(I,K)*AB(K,J)
1337  160          CONTINUE
1338C
1339               DO 170 K=L,N
1340                  AB(K,J)=AB(K,J)+S*AB(K,I)
1341  170          CONTINUE
1342  180       CONTINUE
1343         ENDIF
1344C
1345         IF(L.GT.N) GO TO 200
1346         DO 190 J=L,N
1347            AB(J,I)=ZERO
1348            AB(I,J)=ZERO
1349  190    CONTINUE
1350C
1351  200    AB(I,I)=ONE
1352         G=E(I)
1353         L=I
1354  210 CONTINUE
1355C
1356C  DIAGONALIZATION OF THE BIDIAGONAL FORM...
1357C
1358      EPS=EPS*X
1359      DO 330 KK=1,N
1360         K=N-KK+1
1361         KT=0
1362C
1363C  LABEL TESTFSPLITTING...
1364C
1365  220    KT=KT+1
1366C
1367         IF(KT.GT.JTHIRT) THEN
1368            E(K)=ZERO
1369            WRITE(LP,230)
1370  230       FORMAT(' QR FAILED.')
1371         ENDIF
1372C
1373         DO 240 LL2=1,K
1374            L2=K-LL2+1
1375            L=L2
1376C
1377C  IF(...) GO TO TESTFCONVERGENCE
1378C
1379            IF(ZABS(E(L)).LE.EPS) GO TO 270
1380C
1381C  IF(...) GO TO CANCELLATION
1382C
1383            IF(ZABS(Q(L-1)).LE.EPS) GO TO 250
1384  240    CONTINUE
1385C
1386C  CANCELLATION OF E(L) IF L IS GREATER THAN 1...
1387C  LABEL CANCELLATION...
1388C
1389  250    C=ZERO
1390         S=ONE
1391         IF(L.GT.K) GO TO 270
1392         DO 260 I=L,K
1393            F=S*E(I)
1394            E(I)=C*E(I)
1395C
1396C  IF(...) GO TO TESTFCONVERGENCE
1397C
1398            IF(ZABS(F).LE.EPS) GO TO 270
1399            G=Q(I)
1400C
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)
1405            ELSE
1406               H=ZERO
1407            ENDIF
1408C
1409            Q(I)=H
1410C
1411            IF(H.EQ.ZERO) THEN
1412               H=ONE
1413               G=ONE
1414            ENDIF
1415C
1416C  THE ABOVE REPLACES Q(I) AND H BY SQUARE ROOT OF (G*G+F*F)
1417C  WHICH MAY GIVE INCORRECT RESULTS IF THE SQUARES UNDERFLOW OR IF
1418C  F = G = 0 .
1419C
1420            C=G/H
1421            S=-F/H
1422  260    CONTINUE
1423C
1424C  LABEL TESTFCONVERGENCE...
1425C
1426  270    Z=Q(K)
1427C
1428C  IF(...) GO TO CONVERGENCE
1429C
1430         IF(L.EQ.K) GO TO 310
1431C
1432C  SHIFT FROM BOTTOM 2*2 MINOR.
1433C
1434         X=Q(L)
1435         Y=Q(K-1)
1436         G=E(K-1)
1437         H=E(K)
1438         F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(TWO*H*Y)
1439         G=ZSQRT(F*F+ONE)
1440C
1441         IF(F.LT.ZERO) THEN
1442            DENOM=F-G
1443         ELSE
1444            DENOM=F+G
1445         ENDIF
1446C
1447         F=((X-Z)*(X+Z)+H*(Y/DENOM-H))/X
1448C
1449C  NEXT QR TRANSFORMATION...
1450C
1451         S=ONE
1452         C=ONE
1453         LPI=L+1
1454         IF(LPI.GT.K) GO TO 300
1455         DO 290 I=LPI,K
1456            G=E(I)
1457            Y=Q(I)
1458            H=S*G
1459            G=G*C
1460C
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)
1465            ELSE
1466               Z=ZERO
1467            ENDIF
1468C
1469            E(I-1)=Z
1470C
1471            IF(Z.EQ.ZERO) THEN
1472               F=ONE
1473               Z=ONE
1474            ENDIF
1475C
1476            C=F/Z
1477            S=H/Z
1478            F=X*C+G*S
1479            G=-X*S+G*C
1480            H=Y*S
1481            Y=Y*C
1482C
1483            DO 280 J=1,N
1484               X=AB(J,I-1)
1485               Z=AB(J,I)
1486               AB(J,I-1)=X*C+Z*S
1487               AB(J,I)=-X*S+Z*C
1488  280       CONTINUE
1489C
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)
1494            ELSE
1495               Z=ZERO
1496            ENDIF
1497C
1498            Q(I-1)=Z
1499C
1500            IF(Z.EQ.ZERO) THEN
1501               F=ONE
1502               Z=ONE
1503            ENDIF
1504C
1505            C=F/Z
1506            S=H/Z
1507            F=C*G+S*Y
1508            X=-S*G+C*Y
1509  290    CONTINUE
1510C
1511  300    E(L)=ZERO
1512         E(K)=F
1513         Q(K)=X
1514C
1515C  GO TO TESTFSPLITTING
1516C
1517         GO TO 220
1518C
1519C  LABEL CONVERGENCE...
1520C
1521  310    IF(Z.LT.ZERO) THEN
1522C
1523C  Q(K) IS MADE NON-NEGATIVE.
1524C
1525            Q(K)=-Z
1526            DO 320 J=1,N
1527               AB(J,K)=-AB(J,K)
1528  320       CONTINUE
1529         ENDIF
1530  330 CONTINUE
1531C
1532      RETURN
1533C
1534C  END MINFIT
1535C
1536      END
1537      SUBROUTINE SORT
1538C
1539C  THIS SUBROUTINE SORTS THE ELEMENTS OF D
1540C  AND THE CORRESPONDING COLUMNS OF V INTO DESCENDING ORDER.
1541C
1542C  "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
1543C  RICHARD P. BRENT, PRENTICE-HALL 1973, PAGES 158-159
1544C
1545      INTEGER N,NL,NF,LP,JPRINT,NMAX,ILLCIN,KTM,NFMAX,JRANCH
1546      INTEGER I,IPI,J,K,NMI
1547C
1548      DOUBLE PRECISION V,X,D,Q0,Q1,DMIN,EPSMCH,FX,H,QD0,QD1,QF1,
1549     *   SMALL,T,XLDT,XM2,XM4,DSEED,SCBD
1550      DOUBLE PRECISION S
1551C
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
1555C
1556      NMI=N-1
1557      IF(NMI.LT.1) GO TO 50
1558      DO 40 I=1,NMI
1559         K=I
1560         S=D(I)
1561         IPI=I+1
1562         IF(IPI.GT.N) GO TO 20
1563C
1564         DO 10 J=IPI,N
1565C
1566            IF(D(J).GT.S) THEN
1567               K=J
1568               S=D(J)
1569            ENDIF
1570   10    CONTINUE
1571C
1572   20    IF(K.GT.I) THEN
1573            D(K)=D(I)
1574            D(I)=S
1575C
1576            DO 30 J=1,N
1577               S=V(J,I)
1578               V(J,I)=V(J,K)
1579               V(J,K)=S
1580   30       CONTINUE
1581         ENDIF
1582   40 CONTINUE
1583C
1584   50 RETURN
1585C
1586C  END SORT
1587C
1588      END
1589      SUBROUTINE RANINI(RVALUE)
1590C
1591C  SUBROUTINE RANINI PERFORMS INITIALIZATION FOR SUBROUTINE RANDOM.
1592C
1593C  "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
1594C  RICHARD P. BRENT, PRENTICE-HALL 1973, PAGES 163-164
1595C
1596      INTEGER JRAN2,I
1597C
1598      DOUBLE PRECISION RVALUE,R,RAN3,DMOD,DABS,RAN1
1599C
1600      COMMON /COMRAN/ RAN3(127),RAN1,JRAN2
1601C
1602      R=DMOD(DABS(RVALUE),8190.0D0)+1
1603      JRAN2=127
1604C
1605   10 IF(JRAN2.GT.0) THEN
1606         JRAN2=JRAN2-1
1607         RAN1=-2.0D0**55
1608C
1609         DO 20 I=1,7
1610            R=DMOD(1756.0D0*R,8191.0D0)
1611            RAN1=(RAN1+(R-DMOD(R,32.0D0))/32.0D0)/256.0D0
1612   20    CONTINUE
1613C
1614         RAN3(JRAN2+1)=RAN1
1615         GO TO 10
1616      ENDIF
1617C
1618      RETURN
1619C
1620C  END RANINI
1621C
1622      END
1623      SUBROUTINE RANDOM(RANVAL)
1624C
1625C  SUBROUTINE RANDOM RETURNS A DOUBLE PRECISION PSEUDORANDOM NUMBER
1626C  UNIFORMLY DISTRIBUTED IN (0,1) (INCLUDING 0 BUT NOT 1).
1627C
1628C  "ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES",
1629C  RICHARD P. BRENT, PRENTICE-HALL 1973, PAGES 163-164
1630C
1631C  BEFORE THE FIRST CALL TO RANDOM, THE USER MUST
1632C  CALL RANINI(R) ONCE (ONLY) WITH R A DOUBLE PRECISION NUMBER
1633C  EQUAL TO ANY INTEGER VALUE.
1634C  BRENT (PAGE 166) USED THE EQUIVALENT OF
1635C     CALL RANINI(4.0D0) .
1636C
1637C  THE ALGORITHM USED IN SUBROUTINE RANDOM RETURNS X(N)/2**56,
1638C  WHERE   X(N) = X(N-1) + X(N-127)  (MOD 2**56) .
1639C  SINCE (1 + X + X**127) IS PRIMITIVE (MOD 2),
1640C  THE PERIOD IS AT LEAST (2**127 - 1), WHICH EXCEEDS 10**38.
1641C
1642C  SEE "SEMINUMERICAL ALGORITHMS", VOLUME 2 OF
1643C  "THE ART OF COMPUTER PROGRAMMING" BY DONALD E. KNUTH,
1644C  ADDISON-WESLEY 1969, PAGES 26, 34, AND 464.
1645C
1646C  X(N) IS STORED IN DOUBLE PRECISION AS  RAN3 = X(N)/2**56 - 1/2,
1647C  AND ALL DOUBLE PRECISION ARITHMETIC IS EXACT.
1648C
1649      INTEGER JRAN2
1650C
1651      DOUBLE PRECISION RANVAL,RAN3,RAN1
1652C
1653      COMMON /COMRAN/ RAN3(127),RAN1,JRAN2
1654C
1655      IF(JRAN2.EQ.0) THEN
1656         JRAN2=126
1657      ELSE
1658         JRAN2=JRAN2-1
1659      ENDIF
1660C
1661      RAN1=RAN1+RAN3(JRAN2+1)
1662      IF(RAN1.LT.0.0D0) THEN
1663         RAN1=RAN1+0.5D0
1664      ELSE
1665         RAN1=RAN1-0.5D0
1666      ENDIF
1667C
1668      RAN3(JRAN2+1)=RAN1
1669      RANVAL=RAN1+0.5D0
1670C
1671      RETURN
1672C
1673C  END RANDOM
1674C
1675      END
1676      DOUBLE PRECISION FUNCTION DRANDM(DL)
1677C
1678C  SIMPLE PORTABLE PSEUDORANDOM NUMBER GENERATOR.
1679C
1680C  DRANDM RETURNS FUNCTION VALUES THAT ARE PSEUDORANDOM
1681C  NUMBERS UNIFORMLY DISTRIBUTED ON THE INTERVAL (0,1).
1682C
1683C  'NUMERICAL MATHEMATICS AND COMPUTING' BY WARD CHENEY AND
1684C  DAVID KINCAID, BROOKS/COLE PUBLISHING COMPANY
1685C  (FIRST EDITION, 1980), PAGE 203
1686C
1687C  AT THE BEGINNING OF EXECUTION, OR WHENEVER A NEW SEQUENCE IS
1688C  TO BE INITIATED, SET DL EQUAL TO AN INTEGER VALUE BETWEEN
1689C  1.0D0 AND 2147483646.0D0, INCLUSIVE.  DO THIS ONLY ONCE.
1690C  THEREAFTER, DO NOT SET OR ALTER DL IN ANY WAY.
1691C  FUNCTION DRANDM WILL MODIFY DL FOR ITS OWN PURPOSES.
1692C
1693C  DRANDM USES A MULTIPLICATIVE CONGRUENTIAL METHOD.
1694C  THE NUMBERS GENERATED BY DRANDM SUFFER FROM THE PARALLEL
1695C  PLANES DEFECT DISCOVERED BY G. MARSAGLIA, AND SHOULD NOT BE
1696C  USED WHEN HIGH-QUALITY RANDOMNESS IS REQUIRED.  IN THAT
1697C  CASE, USE A "SHUFFLING" METHOD.
1698C
1699      DOUBLE PRECISION DL,DMOD
1700C
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
1704      RETURN
1705      END
1706