home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
ftp.barnyard.co.uk
/
2015.02.ftp.barnyard.co.uk.tar
/
ftp.barnyard.co.uk
/
cpm
/
walnut-creek-CDROM
/
SIMTEL
/
CPMUG
/
CPMUG041.ARK
/
SSPLIB.FOR
< prev
next >
Wrap
Text File
|
1985-02-10
|
12KB
|
527 lines
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