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