SUBROUTINE SOLVIT(A,Y,R,X,W,V,B,M,N,NR,IE,LUN)
C REVISED DEC 4,78
C
C  FOR LINEAR LEAST-SQUARES CURVE FITTING OR SIMULTANEOUS SOLUTION
C  TO LINEAR EQUATIONS
C-----------------------------------------------------------------------
C  DESCRIPTION OF PARAMETERS --
C   A--   TWO-DIMENSIONAL ARRAY, LENGTH M BY N, BUT DIMENSION IT NR BY
C          THE MAXIMUM VALUE OF N.
C         FOR CURVE FITTING, A IS THE INDEPENDENT-VARIABLE MATRIX,
C          E.G. FOR THE CURVE  LOG Y = X1 + X2/T + X3 *LOG T  PUT
C          1.0,1.0/T(1), AND LOG(T(1)) IN THE 1ST 3 COL. OF ROW 1,
C          1.0,1.0/T(2), AND LOG(T(2)) IN THE 1ST 3 COL. OF ROW 2,
C          ETC., FOR M DATA SETS THERE WILL BE M ROWS AND N IS 3.
C         FOR THE SIMULTANEOUS SOLUTION, EACH ROW CONTAINS THE
C          COEFFICIENTS FOR A PARTICULAR EQUATION.
C   Y--   ONE-DIMENSIONAL ARRAY OF LENGTH M, BUT DIMENSION IT NR.
C          FOR CURVE FITTING,STORE THE DEPENDENT VARIABLES
C          OR THEIR TRANSFORMS.
C          FOR SIMULTANEOUS SOLUTION STORE THE CONSTANT TERMS.
C   R--   ONE-DIMENSIONAL ARRAY OF LENGTH M, BUT DIMENSION IT NR.
C          SOLVIT STORES THE RESIDUALS, THE DIFFERENCE BETWEEN
C          INPUT Y VALUES AND CALCULATED VALUES.
C   X--   ONE-DIMENSIONAL ARRAY, LENGTH N.  DIMENSION IT MAXIMUM N.
C          FOR CURVE FITTING,X IS THE ARRAY OF COEFFICIENTS.
C          FOR SIMULTANEOUS SOLUTION,X IS THE SOLUTION ARRAY.
C   W,V,B  ARE ONE-DIMENSIONAL WORKING ARRAYS, DIMENSION RESPECTIVELY
C          NR,NR, AND NR TIMES THE MAXIMUM VALUE OF M.
C        TO SAVE CORE SPACE, THESE ARRAYS MAY BE EQUIVALENCED TO ARRAYS
C        USED ONLY BEFORE OR AFTER CALLING SOLVIT, E.G. YCALC.
C  THE SUM OF RESIDUALS SQUARED APPEARS IN THE FIRST WORD OF ARRAY B.
C
C CONSTANTS
C   M--   NUMBER OF ROWS BEING USED IN MATRIX A.
C   N--   NUMBER OF COLUMNS BEING USED IN MATRIX A.
C   NR--  NUMBER OF ROWS MATRIX A IS DIMENSIONED IN CALLING PROGRAM.
C   IE--  ERROR CODE-- 0-NO ERROR, 1- ERROR STATEMENT WRITTEN
C
C   TYPICAL DIMENSION AND CALL STATEMENTS
C      TO FIT TO A MAXIMUM OF 250 DATA PAIRS AND A MAXIMUM OF 10 TERMS
C          REAL A(250,10),Y(250),R(250),X(10),W(250),V(250),B(2500)
C          CALL SOLVIT(A,Y,R,X,W,V,B,N,M,250,IE)
C      TO SOLVE A MAXIMUM OF 80 SIMULTANEOUS EQUATIONS
C          REAL A(80,80),Y(80),R(80),X(80),W(80),V(80),B(6400)
C          CALL SOLVIT(A,Y,R,X,W,V,B,N,N,80,IE)
C
C  DEVELOPED BY A. R. MILLER
C-----------------------------------------------------------------------
C
	REAL A(NR,1),Y(1),R(1),X(1),W(1),V(1),B(1)
C
	IE = 0
	IF (N .GT. M) GO TO 304
	DO 10 I = 1, M
	DO 10 J = 1, N
	LL = (J - 1) * M + I
10	B(LL) = A(I,J)
	DO 70  I=1, N
	IIA = (I - 1) * M
	T = 0.
	DO 20  K = I, M
	LL = IIA + K
20	T = T + B(LL) * B(LL)
	IF (T .EQ. 0.0) GO TO 300
	LL = IIA + I
C	W(I) = SIGN(ST , B(LL))
	W(I) = SQRT(T)
	IF (B(LL) .LT. 0.0) W(I)= -W(I)
	B(LL) = B(LL) + W(I)
	V(I) = - B(LL) * W(I)
	IF (I .EQ. N) GO TO 75
	IP = I + 1
	DO 70  J = IP, N
	JIA = (J - 1) * M
	T = 0.
	DO 40  K = I, M
	LL = IIA + K
	LLL = JIA + K
40	T = T + B(LL) * B(LLL)
	T = T / V(I)
	DO 70  L = I , M
	LL = JIA + L
	LLL = IIA + L
70	B(LL) = B(LL) + T * B(LLL)
75	S = 0.
	DO 100 L = 1, M
	R(L)=Y(L)
100	S = S + R(L) * R(L)
	DO 102 I = 1, N
102	X(I) = 0.
	YE = S
105	DO 130 I = 1, N
	IIA = (I - 1) * M
	T = 0.
	DO 110 K = I, M
	LL = IIA + K
110	T = T + B(LL) * R(K)
	T = T / V(I)
	DO 120 L = I, M
	LL = IIA + L
120	R(L) = R(L) + T * B(LL)
130	CONTINUE
	IN = N
	DO 190 I = 1, N
	R(IN) = R(IN) / W(IN)
	INM = IN - 1
	LLIN = INM * M
	J = INM
	DO 155 L = 1, INM
	LL = LLIN + J
	R(J) = R(J) + R(IN) * B(LL)
155	J = J - 1
190	IN = IN - 1
	DO 210 I = 1, N
210	X(I) = X(I) - R(I)
	DO 230 I = 1, M
	RD = Y(I)
	DO 220 J = 1, N
220	RD = RD - A(I,J) * X(J)
	R(I) = RD
230	CONTINUE
	RE = 0.
	DO 240 I = 1, M
240	RE = RE + R(I) * R(I)
	IF (RE .GE. S) GO TO 270
	S = RE
	GO TO 105
270	B(1) = RE
	IF (RE .LT. YE) RETURN
	WRITE (LUN,1001)
	GO TO 390
300	WRITE (LUN,1002)
	GO TO 390
304	WRITE (LUN,1003)
390	IE = 1
	RETURN
1001	FORMAT('0SOLVIT ERROR--R.GT.Y,'
     * ' LOOK FOR WRONG VALUES IN MATRIX')
1002	FORMAT(30H0SOLVIT ERROR--MATRIX SINGULAR)
1003	FORMAT(37H0SOLVIT ERROR--MORE COLUMNS THAN ROWS)
	END
	SUBROUTINE PPLOT(X,Y,YCALC,N,LI,ILINES)
C NOV 14,78  NARROW WIDTH
C-------------------------------------------------------------------
C  PURPOSE--
C   PRODUCE A PLOT ON THE LINE PRINTER OF ONE OR TWO ONE-DIMENSIONAL
C   ARRAYS, E.G., Y AND YCALC, AS A FUNCTION OF A THIRD ONE-DIMENSIONAL
C   ARRAY, E,G., X.  ARRAYS MUST BE ARRANGED IN ASCENDING OR DESCENDING
C   ORDER OF INDEPENDENT VARIABLE X.
C
C
C  DESCRIPTION OF PARAMETERS
C    X -- ARRAY OF INDEPENDENT VARIABLES.
C    Y -- ARRAY OF DEPENDENT VARIABLES, X SYMBOL,OR M WHEN MULTIPLE.
C    YCALC-- ARRAY OF CALCULATED DEPENDENT VARIABLES, * SYMBOL.
C    N -- LENGTH OF ARRAYS X,Y, AND YCALC BEING USED.
C	  MAKE N NEGATIVE IF ONLY X AND Y ARE TO BE PLOTTED.
C	  YCALC IS THEN A DUMMY VARIABLE.
C    LI - LOGICAL UNIT NUMBER FOR OUTPUT
C    ILINES - NUMBER OF LINES FOR PLOT
C
C  WRITTEN BY ALAN R. MILLER
C--------------------------------------------------------------------
C
	REAL X(1),Y(1),YCALC(1),YLABEL(6)
	INTEGER	OUT(51),PLUS,STAR,BLANK,MULT
	DATA IPLUS/1H+/,STAR/1H*/,BLANK/1H /,MULT/1HM/,LINEL/51/
	JF(P)=IFIX((P-YMIN)/YSCALE+1.0)
C
	LINES=MAX0(ILINES,MIN0(N,200))
C	GOTO 1
C	ENTRY PPLOTL(X,Y,YCALC,N,JLINE)
C	LINES=JLINE
1	LONE=0
	PLUS=IPLUS
	IF (N.GT.0) GOTO 2
C  IF N IS NEGATIVE PLOT ONLY X AND Y
	N=-N
	LONE=1
	PLUS=STAR
	WRITE(LI,100) N,PLUS,MULT
2	XLOW=X(1)
	IF (LONE.EQ.0) WRITE(LI,101) N,PLUS,STAR,MULT
	XHI=X(N)
	YMAX=Y(1)
	YMIN=YMAX
	XSCALE=(XHI-XLOW)/FLOAT(LINES-1)
	IF (LONE.EQ.1) GOTO 5
	DO 3 I=1,N
	YMIN=AMIN1(YMIN,Y(I),YCALC(I))
3	YMAX=AMAX1(YMAX,Y(I),YCALC(I))
	GOTO 7
5	DO 6 I=1,N
	YMIN=AMIN1(YMIN,Y(I))
6	YMAX=AMAX1(YMAX,Y(I))
7	YSCALE=(YMAX-YMIN)/50
	YS10=YSCALE*10.0
	YLABEL(1)=YMIN
	DO 9  KN=1,4
9	YLABEL(KN+1)=YLABEL(KN)+YS10
	YLABEL(6)=YMAX
	IPRINT=0
	IF (ABS(XHI).GE.1.E5.OR. ABS(XHI).LT.1.E-2)  IPRINT=1
	DO 10 I=1,LINEL
10	OUT(I)=BLANK
C  FIRST LINE.
	JP=JF(Y(1))
	OUT(JP)=PLUS
	IF (LONE.EQ.1) GOTO 12
	JP=JF(YCALC(1))
	OUT(JP)=STAR
12	L=1
	XLABEL=XLOW
	ISKIP=0
	DO 80 I=2,LINES
	XNEXT =XLOW+XSCALE*FLOAT(I-1)
	IF (ISKIP.EQ.1) GOTO 25
13	L=L+1
	CHANGE=XNEXT-0.5*XSCALE
	IF ((X(L)-CHANGE)*SIGN(1.0,XSCALE).GT.0.0) GOTO 30
C  MULTIPLE POINT
	JP=JF(Y(L))
	IF (OUT(JP).EQ.PLUS.OR.OUT(JP).EQ.MULT) GOTO 15
	OUT(JP)=PLUS
	GOTO 20
15	OUT(JP)=MULT
20	IF (LONE.EQ.1) GOTO 13
	JP=JF(YCALC(L))
	OUT(JP)=STAR
	GOTO 13
C  SKIP LINE
25	WRITE(LI,103)
	GOTO 40
C  PRINT LINE.
30	DO 31 IX=1,LINEL
	JX = LINEL-IX+1
	IF (OUT(JX).NE.BLANK) GOTO 32
31	CONTINUE
	JX = 1
32	IF (IPRINT.EQ.0) WRITE (LI,104) XLABEL,(OUT(IX),IX=1,JX)
	IF (IPRINT.EQ.1) WRITE (LI,102) XLABEL,(OUT(IX),IX=1,JX)
	DO 35 IX=1,LINEL
35	OUT(IX)=BLANK
40	CHANGE=XNEXT + 0.5*XSCALE
	IF ((X(L)-CHANGE)*SIGN(1.0,XSCALE).LT.0.0) GOTO 50
	ISKIP=1
	GOTO 80
50	ISKIP=0
	XLABEL=XNEXT
	JP=JF(Y(L))
	OUT(JP)=PLUS
	IF (LONE.EQ.1) GOTO 80
	JP=JF(YCALC(L))
	OUT(JP)=STAR
80	CONTINUE
81	IF (L.GE.N) GOTO 90
C  MULTIPLE POINT FOR LAST LINE
	L=L+1
	JP=JF(Y(L))
	IF (OUT(JP).EQ.PLUS.OR.OUT(JP).EQ.MULT) GOTO 83
	OUT(JP)=PLUS
	GOTO 81
83	OUT(JP)=MULT
	GOTO 81
90	DO 91 IX=1,LINEL
	JX = LINEL-IX+1
	IF (OUT(JX).NE.BLANK) GOTO 92
91	CONTINUE
	JX = 1
92	IF (IPRINT.EQ.0) WRITE (LI,104) XHI,(OUT(IX),IX=1,JX)
	IF (IPRINT.EQ.1) WRITE (LI,102) XHI,(OUT(IX),IX=1,JX)
	WRITE (LI,107)
	IF (ABS(YMAX).LT.1.E4.AND.ABS(YMAX).GE.1.E-2) GOTO 96
	WRITE (LI,106) YLABEL
	GOTO 97
96	WRITE (LI,108) YLABEL
97	WRITE (LI,109)
	RETURN
100	FORMAT(1H1,I3,10H DATA SETS,6X,
     $ A1,7H-DATA, ,A1,15H-MULTIPLE POINT/)
101	FORMAT(1H1,I3,10H DATA SETS,6X,
     $ A1,7H-DATA, ,A1,15H-FITTED CURVE, ,A1,15H-MULTIPLE POINT/)
102	FORMAT(1X ,1PE11.4,5X,101A1)
103	FORMAT(1X)
104	FORMAT(1X ,F11.4,5X,101A1)
106	FORMAT(1H0,8X,1P11E10.2)
107	FORMAT(17X,5(1H.9X),1H.)
108	FORMAT(1H0,9X,11F10.4)
109	FORMAT(/30X,18HDEPENDENT VARIABLE)
	END
	SUBROUTINE SORT2(X,Y,N)
C
C SHELL-METZNER SORT
C
	REAL X(1),Y(1)
C
	IREV = 0
	IF (N .GT. 0) GOTO 5
	N = -N
	IREV = 1
5	J4 = N
10	J4 = J4 / 2
	IF (J4.EQ.0) RETURN
	J2 = N - J4
	J = 1
20	I = J
30	J3 = I + J4
	IF (IREV .EQ. 0 .AND. X(I) .LE. X(J3)) GOTO 40
	IF (IREV .EQ. 1 .AND. X(I) .GE. X(J3)) GOTO 40
	HOLD = X(I)
	X(I) = X(J3)
	X(J3) = HOLD
	HOLD = Y(I)
	Y(I) = Y(J3)
	Y(J3) = HOLD
	I = I - J4
	IF (I .GE. 1) GOTO 30
40	J = J + 1
	IF (J .GT. J2) GOTO 10
	GOTO 20
	END
      SUBROUTINE SQUARE(A,B,M,N,NA,MB,E,F,G)
C  PROGRAM TO CONVERT CURVE FITTING DATA ARRAY TO SQUARE ARRAY
C  OBTAIN B = TRANSPOSE A * A
C
C  A - M-BY-N CURVE-FITTING ARRAY
C  B - N-BY-N SQUARE ARRAY, A TRANSPOSE * A
C  M - ROWS OF A
C  N - COLUMNS OF A
C  NA- NUMBER OF ROWS ARRAY A IS DIMENSIONED
C  NB- NUMBER OF ROWS ARRAY B IS DIMENSIONED
C  E - DUMMY ARRAY, LENGTH M*M
C  F,G DUMMY ARRAYS, LENGTH M
C
C SUBROUTINE MATINV IS NEEDED
C
	REAL A(NA,1),B(MB,1),E(1),F(1),G(1)
C
	DO 50 K=1,N
	DO 50 L=1,K
	B(K,L)=0.0
	DO 30 I=1,M
30	B(K,L)=B(K,L)+A(I,L)*A(I,K)
	IF (K.NE.L) B(L,K)=B(K,L)
50	CONTINUE
C  INVERT MATRIX B
	CALL MATINV(B,N,MB,E,F,G)
	RETURN
	END
	SUBROUTINE MATINV(B,N,NR,A,L,M)
C
C PURPOSE - TO INVERT A MATRIX
C
C  B - ORIGINAL N-BY-N MATRIX, REPLACED BY INVERSE
C  N - ORDER OF MATRIX
C  NR- NUMBER OF ROWS B IS DIMENSIONED
C  D - RESULTANT DETERMINANT
C  A - WORK ARRAY, LENGTH N*N
C  L,M - WORK ARRAYS, LENGTH N
C
	DIMENSION B(NR,1),A(1),L(1),M(1)
C
C CONVERT SQUARE ARRAY B TO LINEAR ARRAY A
C
	K=0
	DO 5 J=1,N
	DO 5 I=1,N
	K=K+1
5	A(K)=B(I,J)
C
C SEARCH FOR LARGEST ELEMENT
C
	D=1.0
	NK=-N
	DO 80 K=1,N
	NK=NK+N
	L(K)=K
	M(K)=K
	KK=NK+K
	BIGA=A(KK)
	DO 20 J=K,N
	IZ=N*(J-1)
	DO 20 I=K,N
	IJ=IZ+I
10	IF(ABS(BIGA)- ABS(A(IJ))) 15,20,20
15	BIGA=A(IJ)
	L(K)=I
	M(K)=J
20	CONTINUE
C
C INTERCHANGE ROWS
C
	J=L(K)
	IF(J-K) 35,35,25
25	KI=K-N
	DO 30 I=1,N
	KI=KI+N
	HOLD=-A(KI)
	JI=KI-K+J
	A(KI)=A(JI)
30	A(JI) =HOLD
C
C INTERCHANGE COLUMNS
C
35	I=M(K)
	IF(I-K) 45,45,38
38	JP=N*(I-1)
	DO 40 J=1,N
	JK=NK+J
	JI=JP+J
	HOLD=-A(JK)
	A(JK)=A(JI)
40	A(JI) =HOLD
C
C DIVIDE COLUMN BY MINUS PIVOT
C PIVOT IS CONTAINED IN BIGA
C
45	IF(BIGA) 48,46,48
46	D=0.0
	RETURN
48	DO 55 I=1,N
	IF(I-K) 50,55,50
50	IK=NK+I
	A(IK)=A(IK)/(-BIGA)
55	CONTINUE
C
C REDUCE MATRIX
C
	DO 65 I=1,N
	IK=NK+I
	IJ=I-N
	DO 65 J=1,N
	IJ=IJ+N
	IF(I-K) 60,65,60
60	IF(J-K) 62,65,62
62	KJ=IJ-I+K
	A(IJ)=A(IK)*A(KJ)+A(IJ)
65	CONTINUE
C
C DIVIDE ROW BY PIVOT
C
	KJ=K-N
	DO 75 J=1,N
	KJ=KJ+N
	IF(J-K) 70,75,70
70	A(KJ)=A(KJ)/BIGA
75	CONTINUE
C
C PRODUCT OF PIVOTS
C
C *** DONT USE D
C	D=D*BIGA
C
C REPLACE PIVOT BY RECIPROCAL
C
	A(KK)=1.0/BIGA
80	CONTINUE
C
C FINAL ROW AND COLUMN INTERCHANGE
C
	K=N
100	K=K-1
	IF(K) 150,150,105
105	I=L(K)
	IF(I-K) 120,120,108
108	JQ=N*(K-1)
	JR=N*(I-1)
	DO 110 J=1,N
	JK=JQ+J
	HOLD=A(JK)
	JI=JR+J
	A(JK)=-A(JI)
110	A(JI) =HOLD
120	J=M(K)
	IF(J-K) 100,100,125
125	KI=K-N
	DO 130 I=1,N
	KI=KI+N
	HOLD=A(KI)
	JI=KI-K+J
	A(KI)=-A(JI)
130	A(JI) =HOLD
	GOTO 100
C
C CONVERT INVERTED MATRIX A TO N-BY-N MATRIX B
C
150	K=0
	DO 160 J=1,N
	DO 160 I=1,N
	K=K+1
160	B(I,J)=A(K)
	RETURN
	END
	SUBROUTINE LINFIT(X,Y,YC,R,N,A,B,COR,SA,SB)
C LINEAR FIT Y=A + B*X, CORR. COEFF
C STD ERROR ON A AND B
C
	REAL X(1),Y(1),YC(1),R(1)
C
	SUMX  = 0
	SUMY  = 0
	SUMXY = 0
	SUMY2 = 0
	SUMX2 = 0
	DO 10 I=1,N
	XI = X(I)
	YI = Y(I)
	SUMX = SUMX+XI
	SUMY = SUMY+YI
	SUMXY = SUMXY + XI*YI
	SUMX2 = SUMX2 + XI*XI
10	SUMY2 = SUMY2 + YI*YI
	DENOM = SUMX2 - SUMX*SUMX/N
	C = SUMXY - SUMX*SUMY/N
	B = C/DENOM
	A = (SUMY - B*SUMX)/N
	COR = C / SQRT(DENOM * (SUMY2 - SUMY*SUMY / N))
	SEE = SQRT((SUMY2 - A*SUMY - B*SUMXY) / (N-2))
	SB = SEE / SQRT(DENOM)
	SA = SB * SQRT(SUMX2/N)
	DO 20 I=1,N
	YC(I) = A + B * X(I)
	R(I) = Y(I) - YC(I)
20	CONTINUE
	RETURN
	END