Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/nlineq/nlin2.for
There are 2 other files named nlin2.for in the archive. Click here to see a list.
C WESTERN MICHIGAN UNIVERSITY
C NLIN2.FOR (FILE NAME ON LIBRARY DECTAPE)
C NLIN2.FOR IS CALLED (BY RUNUUO) FROM NLINEQ.FOR
C FORWMU PROGS. USED: RUNUUO, MINVSQ
C INTERNAL SUBR. USED: NS01A
C EXTERNAL SUBR. USED: (GENERATED BY NLINEQ.FOR)
C CALFUN.F4
C ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C
C
C***********************************************************************
C***********************************************************************
C
DIMENSION X(50),F(50),AJINV(50,50),W(5250),NORM1(50),NORM2(50)
ITEMP=20
C
C***********************************************************************
C READ IN NECESSARY INFORMATION FROM THE TEMPORARY FILE
C GENERATED BY NLINEQ.F4
C***********************************************************************
C
READ(ITEMP) IDLG,ICC,IOUT,IDSK,N
C
C***********************************************************************
C GATHER ALL THE PARAMETERS TO BE USED IN NS01A
C***********************************************************************
C
300 WRITE(IDLG,30)
30 FORMAT('-ESTIMATED SOLUTION'/' 10 NUMBERS PER LINE,
1 SEPARATED BY COMMAS'/)
READ(ICC,31)(X(I),I=1,N)
31 FORMAT(10F)
WRITE(IDLG,32)
32 FORMAT(/' STEPLENGTH TO BE USED TO APPROXIMATE FIRST DERIVATIVE
1 : ',$)
READ(ICC,31) DSTEP
WRITE(IDLG,33)
33 FORMAT(/' APPROXIMATE DISTANCE OF SOLUTION FROM ESTIMATED
1 SOLUTION: ',
1 $)
READ(ICC,31) DMAX
WRITE(IDLG,34)
34 FORMAT(/' ACCURACY REQUIRED: ',$)
READ(ICC,31) ACC
WRITE(IDLG,35)
35 FORMAT(/' NUMBER OF CALCULATIONS TO BE PERFORMED: ',$)
READ(ICC,36) MAXFUN
36 FORMAT(I)
WRITE(IDLG,37)
37 FORMAT(/' TYPE OUT OPTION: 1 FOR ALL CALCULATIONS'/
1 18X, '0 FOR FINAL CALCULATION ONLY'/)
READ(ICC,36) IPRINT
CALL NS01A(N,X,F,AJINV,DSTEP,DMAX,ACC,MAXFUN,IPRINT,W,NORM1,
1 NORM2,IDLG,ICC,IOUT)
C
C***********************************************************************
C UPON RETURN FROM THE SUBROUTINE, THE USER INTENDS TO
C USE THE SAME FUNCTIONS BUT CHANGE THE PARAMETERS
C***********************************************************************
C
GO TO 300
END
C
C---------------N, X, DSTEP, DMAX, ACC, MAXFUN, IPRINT, IDLG,
C--------------- ICC, IOUT ARE INPUT. X IS MODIFIED AND RETURNED.
C--------------- AJINV, W, ARE OUTPUT. SPACE IS PROVIDED FOR
C--------------- NORM1, NORM2 TO BE USED IN SUBR. MINVSQ.
C SUBROUTINE NS01A(N,X,F,AJINV,DSTEP,DMAX,ACC,MAXFUN,IPRINT,W,
C NORM1,NORM2,IDLG,ICC,IOUT)
C
C
C THIS SUBROUTINE IS ESSENTIALLY THE SAME AS THE ORIGINAL, THE
C ONLY CHANGES MADE ARE:
C (1) REPORT FORMAT
C (2) SUBROUTINE MINVSQ WAS SUBSTITUTED SINCE THE ORIGINAL
C ONE WAS NOT AVAILABLE
C
C***********************************************************************
C
SUBROUTINE NS01A (N,X,F,AJINV,DSTEP,DMAX,ACC,MAXFUN,IPRINT,W,
1 NORM1,NORM2,IDLG,ICC,IOUT)
DIMENSION X(1),F(1),AJINV(N,N),W(1),NORM1(1),NORM2(1)
DATA STARS/'*****'/
WRITE(IOUT,100)(STARS,J=1,13)
100 FORMAT('1'/' CALCU*'/' LATION* I',8X,'X(I)',15X,'F(I)',
1 10X,'SUM OF SQUARES'/' ****',13A5)
C SET VARIOUS PARAMETERS
MAXC=0
C 'MAXC' COUNTS THE NUMBER OF CALLS OF CALFUN
NT=N+4
NTEST=NT
C 'NT' AND 'NTEST' CAUSE AN ERROR RETURN IF F(X) DOES NOT DECREASE
DTEST=FLOAT(N+N)-0.5
C 'DTEST' IS USED TO MAINTAIN LINEAR INDEPENDENCE
NX=N*N
NF=NX+N
NW=NF+N
MW=NW+N
NDC=MW+N
ND=NDC+N
C THESE PARAMETERS SEPARATE THE WORKING SPACE ARRAY W
FMIN=0.
C USUALLY 'FMIN' IS THE LEAST CALCULATED VALUE OF F(X),
C AND THE BEST X IS IN W(NX+1) TO W(NX+N)
DD=0.
C USUALLY DD IS THE SQUARE OF THE CURRENT STEP LENGTH
DSS=DSTEP*DSTEP
DM=DMAX*DMAX
DMM=4.*DM
C---------------SEE ST. 4.
IS=5
C 'IS' CONTROLS A 'GO TO' STATEMENT FOLLOWING A CALL OF CALFUN
TINC=1.
C 'TINC' IS USED IN THE CRITERION TO INCREASE THE STEP LENGTH
C CALL THE SUBROUTINE CALFUN
1 MAXC=MAXC+1
C---------------N, X ARE INPUT. F IS RETURNED.. THIS IS AN EXTERNAL
C--------------- SUBR. GENERATED BY NLINEQ.FOR.
CALL CALFUN (N,X,F)
C TEST FOR CONVERGENCE
FSQ=0.
DO 2 I=1,N
FSQ=FSQ+F(I)*F(I)
2 CONTINUE
IF (FSQ.GT.ACC) GO TO 4
C PROVIDE PRINTING OF FINAL SOLUTION IF REQUESTED
3 WRITE(IOUT,300) MAXC,(I,X(I),F(I),I=1,2)
300 FORMAT(7X,'*'/1X,I5,' *',I4,2(2X,E17.8)/1X,'FINAL *',I4,
1 2(2X,E17.8))
IF (N.GT.2) WRITE(IOUT,261)(I,X(I),F(I),I=3,N)
WRITE(IOUT,262) FSQ
C
C***********************************************************************
C END OF ONE SOLUTION, DETERMINE USER'S INTENTION
C
C IF LAST = 0 TERMINATE THE JOB
C = 1 BRANCH TO BEGINNING OF THE PROGRAM
C = 2 RETAIN THE SAME FUNCTIONS BUT CHANGE PARAMETERS
C***********************************************************************
C
5 WRITE(IDLG,500)
500 FORMAT('-END OF JOB, TYPE 0 TO TERMINATE'/18X,'1 TO BRANCH TO
1 BEGINNING OF PROGRAM'/18X,'2 TO CHANGE PARAMETERS'/)
READ(ICC,501) LAST
501 FORMAT(I)
IF (LAST-1) 502,504,503
502 CALL EXIT
503 RETURN
C---------------NLINEQ.EXE IS IN SYS:
504 CALL RUNUUO('R NLINEQ')
C TEST FOR ERROR RETURN BECAUSE F(X) DOES NOT DECREASE
4 GO TO (10,11,11,10,11),IS
10 IF (FSQ-FMIN) 15,20,20
20 IF (DD-DSS) 12,12,11
12 NTEST=NTEST-1
IF (NTEST) 13,14,11
14 WRITE(IOUT,16) NT
16 FORMAT('-'/5X,'ERROR: ',I5,' CALCULATIONS FAILED TO IMPROVE
1 THE RESIDUALS')
17 DO 18 I=1,N
X(I)=W(NX+I)
F(I)=W(NF+I)
18 CONTINUE
FSQ=FMIN
GO TO 3
C ERROR RETURN BECAUSE A NEW JACOBIAN IS UNSUCCESSFUL
13 WRITE(IOUT,19)
19 FORMAT('-'/5X,'ERROR: F(X) FAILED TO DECREASE USING A NEW
1 JACOBIAN')
GO TO 17
15 NTEST=NT
C TEST WHETHER THERE HAVE BEEN MAXFUN CALLS OF CALFUN
11 IF (MAXFUN-MAXC) 21,21,22
21 WRITE(IOUT,23) MAXC
23 FORMAT('-'/5X,'ERROR: THERE HAVE BEEN ',I5, ' CALCULATIONS
1 PERFORMED')
IF (FSQ.GE.FMIN) GO TO 17
GO TO 3
C PROVIDE PRINTING IF REQUESTED
22 IF (IPRINT) 24,24,25
25 WRITE(IOUT,26) MAXC,X(1),F(1)
26 FORMAT(7X,'*'/1X,I5,' * 1',2(2X,E17.8))
WRITE(IOUT,261) (I,X(I),F(I),I=2,N)
261 FORMAT(7X,'*',I4,2X,E17.8,2X,E17.8)
WRITE(IOUT,262) FSQ
262 FORMAT(7X,'*',44X,E17.8)
24 GO TO (27,28,29,87,30),IS
C STORE THE RESULT OF THE INITIAL CALL OF CALFUN
30 FMIN=FSQ
DO 31 I=1,N
W(NX+I)=X(I)
W(NF+I)=F(I)
31 CONTINUE
C CALCULATE A NEW JACOBIAN APPROXIMATION
32 IC=0
IS=3
33 IC=IC+1
X(IC)=X(IC)+DSTEP
GO TO 1
29 K=IC
DO 34 I=1,N
W(K)=(F(I)-W(NF+I))/DSTEP
K=K+N
34 CONTINUE
X(IC)=W(NX+IC)
IF (IC-N) 33,35,35
C CALCULATE THE INVERSE OF THE JACOBIAN AND SET THE DIRECTION MATRIX
35 K=0
DO 36 I=1,N
DO 37 J=1,N
K=K+1
AJINV(I,J)=W(K)
W(ND+K)=0.
37 CONTINUE
W(NDC+K+I)=1.
W(NDC+I)=1.+FLOAT(N-I)
36 CONTINUE
C
C***********************************************************************
C WMU COMPUTER CENTER'S OWN MATRIX INVERSE SUBROUTINE
C*********************************************************************
C
C---------------AJINV, N, MESS ARE INPUT. DET, IEXP ARE OUTPUT.
C--------------- SPACE IS PROVIDED FOR NORM1, NORM2 BY ST. 100-3,ARG.
C--------------- LIST OF SUBR. NS01A.
CALL MINVSQ(AJINV,N,1.,NORM1,NORM2,N,MESS,1,DET,IEXP)
C START ITERATION BY PREDICTING THE DESCENT AND NEWTON MINIMA
38 DS=0.
DN=0.
SP=0.
DO 39 I=1,N
X(I)=0.
F(I)=0.
K=I
DO 40 J=1,N
X(I)=X(I)-W(K)*W(NF+J)
F(I)=F(I)-AJINV(I,J)*W(NF+J)
K=K+N
40 CONTINUE
DS=DS+X(I)*X(I)
DN=DN+F(I)*F(I)
SP=SP+X(I)*F(I)
39 CONTINUE
C TEST WHETHER A NEARBY STATIONARY POINT IS PREDICTED
IF (FMIN*FMIN-DMM*DS) 41,41,42
C IF SO THEN RETURN OR REVISE JACOBIAN
42 GO TO (43,43,44),IS
44 WRITE(IOUT,45)
45 FORMAT('-'/5X,'ERROR: A NEARBY STATIONARY POINT OF F(X) IS
1 PREDICTED')
GO TO 17
43 NTEST=0
DO 46 I=1,N
X(I)=W(NX+I)
46 CONTINUE
GO TO 32
C TEST WHETHER TO APPLY THE FULL NEWTON CORRECTION
41 IS=2
IF (DN-DD) 47,47,48
47 DD=AMAX1(DN,DSS)
DS=0.25*DN
TINC=1.
IF (DN-DSS) 49,58,58
49 IS=4
GO TO 80
C CALCULATE THE LENGTH OF THE STEEPEST DESCENT STEP
48 K=0
DMULT=0.
DO 51 I=1,N
DW=0.
DO 52 J=1,N
K=K+1
DW=DW+W(K)*X(J)
52 CONTINUE
DMULT=DMULT+DW*DW
51 CONTINUE
DMULT=DS/DMULT
DS=DS*DMULT*DMULT
C TEST WHETHER TO USE THE STEEPEST DESCENT DIRECTION
IF (DS-DD) 53,54,54
C TEST WHETHER THE INITIAL VALUE OF DD HAS BEEN SET
54 IF (DD) 55,55,56
55 DD=AMAX1(DSS,AMIN1(DM,DS))
DS=DS/(DMULT*DMULT)
GO TO 41
C SET THE MULTIPLIER OF THE STEEPEST DESCENT DIRECTION
56 ANMULT=0.
DMULT=DMULT*SQRT(DD/DS)
GO TO 98
C INTERPOLATE BETWEEN THE STEEPEST DESCENT AND THE NEWTON DIRECTIONS
53 SP=SP*DMULT
ANMULT=(DD-DS)/((SP-DS)+SQRT((SP-DD)**2+(DN-DD)*(DD-DS)))
DMULT=DMULT*(1.-ANMULT)
C CALCULATE THE CHANGE IN X AND ITS ANGLE WITH THE FIRST DIRECTION
98 DN=0.
SP=0.
DO 57 I=1,N
F(I)=DMULT*X(I)+ANMULT*F(I)
DN=DN+F(I)*F(I)
SP=SP+F(I)*W(ND+I)
57 CONTINUE
DS=0.25*DN
C TEST WHETHER AN EXTRA STEP IS NEEDED FOR INDEPENDENCE
IF (W(NDC+1)-DTEST) 58,58,59
59 IF (SP*SP-DS) 60,58,58
C TAKE THE EXTRA STEP AND UPDATE THE DIRECTION MATRIX
50 IS=2
60 DO 61 I=1,N
X(I)=W(NX+I)+DSTEP*W(ND+I)
W(NDC+I)=W(NDC+I+1)+1.
61 CONTINUE
W(ND)=1.
DO 62 I=1,N
K=ND+I
SP=W(K)
DO 63 J=2,N
W(K)=W(K+N)
K=K+N
63 CONTINUE
W(K)=SP
62 CONTINUE
GO TO 1
C EXPRESS THE NEW DIRECTION IN TERMS OF THOSE OF THE DIRECTION
C MATRIX, AND UPDATE THE COUNTS IN W(NDC+1) ETC.
58 SP=0.
K=ND
DO 64 I=1,N
X(I)=DW
DW=0.
DO 65 J=1,N
K=K+1
DW=DW+F(J)*W(K)
65 CONTINUE
GO TO (68,66),IS
66 W(NDC+I)=W(NDC+I)+1.
SP=SP+DW*DW
IF (SP-DS) 64,64,67
67 IS=1
KK=I
X(1)=DW
GO TO 69
68 X(I)=DW
69 W(NDC+I)=W(NDC+I+1)+1.
64 CONTINUE
W(ND)=1.
C REORDER THE DIRECTIONS SO THAT KK IS FIRST
IF (KK-1) 70,70,71
71 KS=NDC+KK*N
DO 72 I=1,N
K=KS+I
SP=W(K)
DO 73 J=2,KK
W(K)=W(K-N)
K=K-N
73 CONTINUE
W(K)=SP
72 CONTINUE
C GENERATE THE NEW ORTHOGONAL DIRECTION MATRIX
70 DO 74 I=1,N
W(NW+I)=0.
74 CONTINUE
SP=X(1)*X(1)
K=ND
DO 75 I=2,N
DS=SQRT(SP*(SP+X(I)*X(I)))
DW=SP/DS
DS=X(I)/DS
SP=SP+X(I)*X(I)
DO 76 J=1,N
K=K+1
W(NW+J)=W(NW+J)+X(I-1)*W(K)
W(K)=DW*W(K+N)-DS*W(NW+J)
76 CONTINUE
75 CONTINUE
SP=1./SQRT(DN)
DO 77 I=1,N
K=K+1
W(K)=SP*F(I)
77 CONTINUE
C CALCULATE THE NEXT VECTOR X, AND PREDICT THE RIGHT HAND SIDES
80 FNP=0.
K=0
DO 78 I=1,N
X(I)=W(NX+I)+F(I)
W(NW+I)=W(NF+I)
DO 79 J=1,N
K=K+1
W(NW+I)=W(NW+I)+W(K)*F(J)
79 CONTINUE
FNP=FNP+W(NW+I)**2
78 CONTINUE
C CALL CALFUN USING THE NEW VECTOR OF VARIABLES
GO TO 1
C UPDATE THE STEP SIZE
27 DMULT=0.9*FMIN+0.1*FNP-FSQ
IF (DMULT) 82,81,81
82 DD=AMAX1(DSS,0.25*DD)
TINC=1.
IF (FSQ-FMIN) 83,28,28
C TRY THE TEST TO DECIDE WHETHER TO INCREASE THE STEP LENGTH
81 SP=0.
SS=0.
DO 84 I=1,N
SP=SP+ABS(F(I)*(F(I)-W(NW+I)))
SS=SS+(F(I)-W(NW+I))**2
84 CONTINUE
PJ=1.+DMULT/(SP+SQRT(SP*SP+DMULT*SS))
SP=AMIN1(4.,TINC,PJ)
TINC=PJ/SP
DD=AMIN1(DM,SP*DD)
GO TO 83
C IF F(X) IMPROVES STORE THE NEXT VALUE OF X
87 IF (FSQ-FMIN) 83,50,50
83 FMIN=FSQ
DO 88 I=1,N
SP=X(I)
X(I)=W(NX+I)
W(NX+I)=SP
SP=F(I)
F(I)=W(NF+I)
W(NF+I)=SP
W(NW+I)=-W(NW+I)
88 CONTINUE
IF (IS-1) 28,28,50
C CALCULATE THE CHANGES IN F AND IN X
28 DO 89 I=1,N
X(I)=X(I)-W(NX+I)
F(I)=F(I)-W(NF+I)
89 CONTINUE
C UPDATE THE APPROXIMATIONS TO J AND TO AJINV
K=0
DO 90 I=1,N
W(MW+I)=X(I)
W(NW+I)=F(I)
DO 91 J=1,N
W(MW+I)=W(MW+I)-AJINV(I,J)*F(J)
K=K+1
W(NW+I)=W(NW+I)-W(K)*X(J)
91 CONTINUE
90 CONTINUE
SP=0.
SS=0.
DO 92 I=1,N
DS=0.
DO 93 J=1,N
DS=DS+AJINV(J,I)*X(J)
93 CONTINUE
SP=SP+DS*F(I)
SS=SS+X(I)*X(I)
F(I)=DS
92 CONTINUE
DMULT=1.
IF (ABS(SP)-0.1*SS) 94,95,95
94 DMULT=0.8
95 PJ=DMULT/SS
PA=DMULT/(DMULT*SP+(1.-DMULT)*SS)
K=0
DO 96 I=1,N
SP=PJ*W(NW+I)
SS=PA*W(MW+I)
DO 97 J=1,N
K=K+1
W(K)=W(K)+SP*X(J)
AJINV(I,J)=AJINV(I,J)+SS*F(J)
97 CONTINUE
96 CONTINUE
GO TO 38
END