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