home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / FORTH / 4THPROG.ZIP / SSPLIB.FOR < prev    next >
Encoding:
Text File  |  1982-05-31  |  12.0 KB  |  527 lines
  1.     SUBROUTINE SOLVIT(A,Y,R,X,W,V,B,M,N,NR,IE,LUN)
  2. C REVISED DEC 4,78
  3. C
  4. C  FOR LINEAR LEAST-SQUARES CURVE FITTING OR SIMULTANEOUS SOLUTION
  5. C  TO LINEAR EQUATIONS
  6. C-----------------------------------------------------------------------
  7. C  DESCRIPTION OF PARAMETERS --
  8. C   A--   TWO-DIMENSIONAL ARRAY, LENGTH M BY N, BUT DIMENSION IT NR BY
  9. C          THE MAXIMUM VALUE OF N.
  10. C         FOR CURVE FITTING, A IS THE INDEPENDENT-VARIABLE MATRIX,
  11. C          E.G. FOR THE CURVE  LOG Y = X1 + X2/T + X3 *LOG T  PUT
  12. C          1.0,1.0/T(1), AND LOG(T(1)) IN THE 1ST 3 COL. OF ROW 1,
  13. C          1.0,1.0/T(2), AND LOG(T(2)) IN THE 1ST 3 COL. OF ROW 2,
  14. C          ETC., FOR M DATA SETS THERE WILL BE M ROWS AND N IS 3.
  15. C         FOR THE SIMULTANEOUS SOLUTION, EACH ROW CONTAINS THE
  16. C          COEFFICIENTS FOR A PARTICULAR EQUATION.
  17. C   Y--   ONE-DIMENSIONAL ARRAY OF LENGTH M, BUT DIMENSION IT NR.
  18. C          FOR CURVE FITTING,STORE THE DEPENDENT VARIABLES
  19. C          OR THEIR TRANSFORMS.
  20. C          FOR SIMULTANEOUS SOLUTION STORE THE CONSTANT TERMS.
  21. C   R--   ONE-DIMENSIONAL ARRAY OF LENGTH M, BUT DIMENSION IT NR.
  22. C          SOLVIT STORES THE RESIDUALS, THE DIFFERENCE BETWEEN
  23. C          INPUT Y VALUES AND CALCULATED VALUES.
  24. C   X--   ONE-DIMENSIONAL ARRAY, LENGTH N.  DIMENSION IT MAXIMUM N.
  25. C          FOR CURVE FITTING,X IS THE ARRAY OF COEFFICIENTS.
  26. C          FOR SIMULTANEOUS SOLUTION,X IS THE SOLUTION ARRAY.
  27. C   W,V,B  ARE ONE-DIMENSIONAL WORKING ARRAYS, DIMENSION RESPECTIVELY
  28. C          NR,NR, AND NR TIMES THE MAXIMUM VALUE OF M.
  29. C        TO SAVE CORE SPACE, THESE ARRAYS MAY BE EQUIVALENCED TO ARRAYS
  30. C        USED ONLY BEFORE OR AFTER CALLING SOLVIT, E.G. YCALC.
  31. C  THE SUM OF RESIDUALS SQUARED APPEARS IN THE FIRST WORD OF ARRAY B.
  32. C
  33. C CONSTANTS
  34. C   M--   NUMBER OF ROWS BEING USED IN MATRIX A.
  35. C   N--   NUMBER OF COLUMNS BEING USED IN MATRIX A.
  36. C   NR--  NUMBER OF ROWS MATRIX A IS DIMENSIONED IN CALLING PROGRAM.
  37. C   IE--  ERROR CODE-- 0-NO ERROR, 1- ERROR STATEMENT WRITTEN
  38. C
  39. C   TYPICAL DIMENSION AND CALL STATEMENTS
  40. C      TO FIT TO A MAXIMUM OF 250 DATA PAIRS AND A MAXIMUM OF 10 TERMS
  41. C          REAL A(250,10),Y(250),R(250),X(10),W(250),V(250),B(2500)
  42. C          CALL SOLVIT(A,Y,R,X,W,V,B,N,M,250,IE)
  43. C      TO SOLVE A MAXIMUM OF 80 SIMULTANEOUS EQUATIONS
  44. C          REAL A(80,80),Y(80),R(80),X(80),W(80),V(80),B(6400)
  45. C          CALL SOLVIT(A,Y,R,X,W,V,B,N,N,80,IE)
  46. C
  47. C  DEVELOPED BY A. R. MILLER
  48. C-----------------------------------------------------------------------
  49. C
  50.     REAL A(NR,1),Y(1),R(1),X(1),W(1),V(1),B(1)
  51. C
  52.     IE = 0
  53.     IF (N .GT. M) GO TO 304
  54.     DO 10 I = 1, M
  55.     DO 10 J = 1, N
  56.     LL = (J - 1) * M + I
  57. 10    B(LL) = A(I,J)
  58.     DO 70  I=1, N
  59.     IIA = (I - 1) * M
  60.     T = 0.
  61.     DO 20  K = I, M
  62.     LL = IIA + K
  63. 20    T = T + B(LL) * B(LL)
  64.     IF (T .EQ. 0.0) GO TO 300
  65.     LL = IIA + I
  66. C    W(I) = SIGN(ST , B(LL))
  67.     W(I) = SQRT(T)
  68.     IF (B(LL) .LT. 0.0) W(I)= -W(I)
  69.     B(LL) = B(LL) + W(I)
  70.     V(I) = - B(LL) * W(I)
  71.     IF (I .EQ. N) GO TO 75
  72.     IP = I + 1
  73.     DO 70  J = IP, N
  74.     JIA = (J - 1) * M
  75.     T = 0.
  76.     DO 40  K = I, M
  77.     LL = IIA + K
  78.     LLL = JIA + K
  79. 40    T = T + B(LL) * B(LLL)
  80.     T = T / V(I)
  81.     DO 70  L = I , M
  82.     LL = JIA + L
  83.     LLL = IIA + L
  84. 70    B(LL) = B(LL) + T * B(LLL)
  85. 75    S = 0.
  86.     DO 100 L = 1, M
  87.     R(L)=Y(L)
  88. 100    S = S + R(L) * R(L)
  89.     DO 102 I = 1, N
  90. 102    X(I) = 0.
  91.     YE = S
  92. 105    DO 130 I = 1, N
  93.     IIA = (I - 1) * M
  94.     T = 0.
  95.     DO 110 K = I, M
  96.     LL = IIA + K
  97. 110    T = T + B(LL) * R(K)
  98.     T = T / V(I)
  99.     DO 120 L = I, M
  100.     LL = IIA + L
  101. 120    R(L) = R(L) + T * B(LL)
  102. 130    CONTINUE
  103.     IN = N
  104.     DO 190 I = 1, N
  105.     R(IN) = R(IN) / W(IN)
  106.     INM = IN - 1
  107.     LLIN = INM * M
  108.     J = INM
  109.     DO 155 L = 1, INM
  110.     LL = LLIN + J
  111.     R(J) = R(J) + R(IN) * B(LL)
  112. 155    J = J - 1
  113. 190    IN = IN - 1
  114.     DO 210 I = 1, N
  115. 210    X(I) = X(I) - R(I)
  116.     DO 230 I = 1, M
  117.     RD = Y(I)
  118.     DO 220 J = 1, N
  119. 220    RD = RD - A(I,J) * X(J)
  120.     R(I) = RD
  121. 230    CONTINUE
  122.     RE = 0.
  123.     DO 240 I = 1, M
  124. 240    RE = RE + R(I) * R(I)
  125.     IF (RE .GE. S) GO TO 270
  126.     S = RE
  127.     GO TO 105
  128. 270    B(1) = RE
  129.     IF (RE .LT. YE) RETURN
  130.     WRITE (LUN,1001)
  131.     GO TO 390
  132. 300    WRITE (LUN,1002)
  133.     GO TO 390
  134. 304    WRITE (LUN,1003)
  135. 390    IE = 1
  136.     RETURN
  137. 1001    FORMAT('0SOLVIT ERROR--R.GT.Y,'
  138.      * ' LOOK FOR WRONG VALUES IN MATRIX')
  139. 1002    FORMAT(30H0SOLVIT ERROR--MATRIX SINGULAR)
  140. 1003    FORMAT(37H0SOLVIT ERROR--MORE COLUMNS THAN ROWS)
  141.     END
  142.     SUBROUTINE PPLOT(X,Y,YCALC,N,LI,ILINES)
  143. C NOV 14,78  NARROW WIDTH
  144. C-------------------------------------------------------------------
  145. C  PURPOSE--
  146. C   PRODUCE A PLOT ON THE LINE PRINTER OF ONE OR TWO ONE-DIMENSIONAL
  147. C   ARRAYS, E.G., Y AND YCALC, AS A FUNCTION OF A THIRD ONE-DIMENSIONAL
  148. C   ARRAY, E,G., X.  ARRAYS MUST BE ARRANGED IN ASCENDING OR DESCENDING
  149. C   ORDER OF INDEPENDENT VARIABLE X.
  150. C
  151. C
  152. C  DESCRIPTION OF PARAMETERS
  153. C    X -- ARRAY OF INDEPENDENT VARIABLES.
  154. C    Y -- ARRAY OF DEPENDENT VARIABLES, X SYMBOL,OR M WHEN MULTIPLE.
  155. C    YCALC-- ARRAY OF CALCULATED DEPENDENT VARIABLES, * SYMBOL.
  156. C    N -- LENGTH OF ARRAYS X,Y, AND YCALC BEING USED.
  157. C      MAKE N NEGATIVE IF ONLY X AND Y ARE TO BE PLOTTED.
  158. C      YCALC IS THEN A DUMMY VARIABLE.
  159. C    LI - LOGICAL UNIT NUMBER FOR OUTPUT
  160. C    ILINES - NUMBER OF LINES FOR PLOT
  161. C
  162. C  WRITTEN BY ALAN R. MILLER
  163. C--------------------------------------------------------------------
  164. C
  165.     REAL X(1),Y(1),YCALC(1),YLABEL(6)
  166.     INTEGER    OUT(51),PLUS,STAR,BLANK,MULT
  167.     DATA IPLUS/1H+/,STAR/1H*/,BLANK/1H /,MULT/1HM/,LINEL/51/
  168.     JF(P)=IFIX((P-YMIN)/YSCALE+1.0)
  169. C
  170.     LINES=MAX0(ILINES,MIN0(N,200))
  171. C    GOTO 1
  172. C    ENTRY PPLOTL(X,Y,YCALC,N,JLINE)
  173. C    LINES=JLINE
  174. 1    LONE=0
  175.     PLUS=IPLUS
  176.     IF (N.GT.0) GOTO 2
  177. C  IF N IS NEGATIVE PLOT ONLY X AND Y
  178.     N=-N
  179.     LONE=1
  180.     PLUS=STAR
  181.     WRITE(LI,100) N,PLUS,MULT
  182. 2    XLOW=X(1)
  183.     IF (LONE.EQ.0) WRITE(LI,101) N,PLUS,STAR,MULT
  184.     XHI=X(N)
  185.     YMAX=Y(1)
  186.     YMIN=YMAX
  187.     XSCALE=(XHI-XLOW)/FLOAT(LINES-1)
  188.     IF (LONE.EQ.1) GOTO 5
  189.     DO 3 I=1,N
  190.     YMIN=AMIN1(YMIN,Y(I),YCALC(I))
  191. 3    YMAX=AMAX1(YMAX,Y(I),YCALC(I))
  192.     GOTO 7
  193. 5    DO 6 I=1,N
  194.     YMIN=AMIN1(YMIN,Y(I))
  195. 6    YMAX=AMAX1(YMAX,Y(I))
  196. 7    YSCALE=(YMAX-YMIN)/50
  197.     YS10=YSCALE*10.0
  198.     YLABEL(1)=YMIN
  199.     DO 9  KN=1,4
  200. 9    YLABEL(KN+1)=YLABEL(KN)+YS10
  201.     YLABEL(6)=YMAX
  202.     IPRINT=0
  203.     IF (ABS(XHI).GE.1.E5.OR. ABS(XHI).LT.1.E-2)  IPRINT=1
  204.     DO 10 I=1,LINEL
  205. 10    OUT(I)=BLANK
  206. C  FIRST LINE.
  207.     JP=JF(Y(1))
  208.     OUT(JP)=PLUS
  209.     IF (LONE.EQ.1) GOTO 12
  210.     JP=JF(YCALC(1))
  211.     OUT(JP)=STAR
  212. 12    L=1
  213.     XLABEL=XLOW
  214.     ISKIP=0
  215.     DO 80 I=2,LINES
  216.     XNEXT =XLOW+XSCALE*FLOAT(I-1)
  217.     IF (ISKIP.EQ.1) GOTO 25
  218. 13    L=L+1
  219.     CHANGE=XNEXT-0.5*XSCALE
  220.     IF ((X(L)-CHANGE)*SIGN(1.0,XSCALE).GT.0.0) GOTO 30
  221. C  MULTIPLE POINT
  222.     JP=JF(Y(L))
  223.     IF (OUT(JP).EQ.PLUS.OR.OUT(JP).EQ.MULT) GOTO 15
  224.     OUT(JP)=PLUS
  225.     GOTO 20
  226. 15    OUT(JP)=MULT
  227. 20    IF (LONE.EQ.1) GOTO 13
  228.     JP=JF(YCALC(L))
  229.     OUT(JP)=STAR
  230.     GOTO 13
  231. C  SKIP LINE
  232. 25    WRITE(LI,103)
  233.     GOTO 40
  234. C  PRINT LINE.
  235. 30    DO 31 IX=1,LINEL
  236.     JX = LINEL-IX+1
  237.     IF (OUT(JX).NE.BLANK) GOTO 32
  238. 31    CONTINUE
  239.     JX = 1
  240. 32    IF (IPRINT.EQ.0) WRITE (LI,104) XLABEL,(OUT(IX),IX=1,JX)
  241.     IF (IPRINT.EQ.1) WRITE (LI,102) XLABEL,(OUT(IX),IX=1,JX)
  242.     DO 35 IX=1,LINEL
  243. 35    OUT(IX)=BLANK
  244. 40    CHANGE=XNEXT + 0.5*XSCALE
  245.     IF ((X(L)-CHANGE)*SIGN(1.0,XSCALE).LT.0.0) GOTO 50
  246.     ISKIP=1
  247.     GOTO 80
  248. 50    ISKIP=0
  249.     XLABEL=XNEXT
  250.     JP=JF(Y(L))
  251.     OUT(JP)=PLUS
  252.     IF (LONE.EQ.1) GOTO 80
  253.     JP=JF(YCALC(L))
  254.     OUT(JP)=STAR
  255. 80    CONTINUE
  256. 81    IF (L.GE.N) GOTO 90
  257. C  MULTIPLE POINT FOR LAST LINE
  258.     L=L+1
  259.     JP=JF(Y(L))
  260.     IF (OUT(JP).EQ.PLUS.OR.OUT(JP).EQ.MULT) GOTO 83
  261.     OUT(JP)=PLUS
  262.     GOTO 81
  263. 83    OUT(JP)=MULT
  264.     GOTO 81
  265. 90    DO 91 IX=1,LINEL
  266.     JX = LINEL-IX+1
  267.     IF (OUT(JX).NE.BLANK) GOTO 92
  268. 91    CONTINUE
  269.     JX = 1
  270. 92    IF (IPRINT.EQ.0) WRITE (LI,104) XHI,(OUT(IX),IX=1,JX)
  271.     IF (IPRINT.EQ.1) WRITE (LI,102) XHI,(OUT(IX),IX=1,JX)
  272.     WRITE (LI,107)
  273.     IF (ABS(YMAX).LT.1.E4.AND.ABS(YMAX).GE.1.E-2) GOTO 96
  274.     WRITE (LI,106) YLABEL
  275.     GOTO 97
  276. 96    WRITE (LI,108) YLABEL
  277. 97    WRITE (LI,109)
  278.     RETURN
  279. 100    FORMAT(1H1,I3,10H DATA SETS,6X,
  280.      $ A1,7H-DATA, ,A1,15H-MULTIPLE POINT/)
  281. 101    FORMAT(1H1,I3,10H DATA SETS,6X,
  282.      $ A1,7H-DATA, ,A1,15H-FITTED CURVE, ,A1,15H-MULTIPLE POINT/)
  283. 102    FORMAT(1X ,1PE11.4,5X,101A1)
  284. 103    FORMAT(1X)
  285. 104    FORMAT(1X ,F11.4,5X,101A1)
  286. 106    FORMAT(1H0,8X,1P11E10.2)
  287. 107    FORMAT(17X,5(1H.9X),1H.)
  288. 108    FORMAT(1H0,9X,11F10.4)
  289. 109    FORMAT(/30X,18HDEPENDENT VARIABLE)
  290.     END
  291.     SUBROUTINE SORT2(X,Y,N)
  292. C
  293. C SHELL-METZNER SORT
  294. C
  295.     REAL X(1),Y(1)
  296. C
  297.     IREV = 0
  298.     IF (N .GT. 0) GOTO 5
  299.     N = -N
  300.     IREV = 1
  301. 5    J4 = N
  302. 10    J4 = J4 / 2
  303.     IF (J4.EQ.0) RETURN
  304.     J2 = N - J4
  305.     J = 1
  306. 20    I = J
  307. 30    J3 = I + J4
  308.     IF (IREV .EQ. 0 .AND. X(I) .LE. X(J3)) GOTO 40
  309.     IF (IREV .EQ. 1 .AND. X(I) .GE. X(J3)) GOTO 40
  310.     HOLD = X(I)
  311.     X(I) = X(J3)
  312.     X(J3) = HOLD
  313.     HOLD = Y(I)
  314.     Y(I) = Y(J3)
  315.     Y(J3) = HOLD
  316.     I = I - J4
  317.     IF (I .GE. 1) GOTO 30
  318. 40    J = J + 1
  319.     IF (J .GT. J2) GOTO 10
  320.     GOTO 20
  321.     END
  322.       SUBROUTINE SQUARE(A,B,M,N,NA,MB,E,F,G)
  323. C  PROGRAM TO CONVERT CURVE FITTING DATA ARRAY TO SQUARE ARRAY
  324. C  OBTAIN B = TRANSPOSE A * A
  325. C
  326. C  A - M-BY-N CURVE-FITTING ARRAY
  327. C  B - N-BY-N SQUARE ARRAY, A TRANSPOSE * A
  328. C  M - ROWS OF A
  329. C  N - COLUMNS OF A
  330. C  NA- NUMBER OF ROWS ARRAY A IS DIMENSIONED
  331. C  NB- NUMBER OF ROWS ARRAY B IS DIMENSIONED
  332. C  E - DUMMY ARRAY, LENGTH M*M
  333. C  F,G DUMMY ARRAYS, LENGTH M
  334. C
  335. C SUBROUTINE MATINV IS NEEDED
  336. C
  337.     REAL A(NA,1),B(MB,1),E(1),F(1),G(1)
  338. C
  339.     DO 50 K=1,N
  340.     DO 50 L=1,K
  341.     B(K,L)=0.0
  342.     DO 30 I=1,M
  343. 30    B(K,L)=B(K,L)+A(I,L)*A(I,K)
  344.     IF (K.NE.L) B(L,K)=B(K,L)
  345. 50    CONTINUE
  346. C  INVERT MATRIX B
  347.     CALL MATINV(B,N,MB,E,F,G)
  348.     RETURN
  349.     END
  350.     SUBROUTINE MATINV(B,N,NR,A,L,M)
  351. C
  352. C PURPOSE - TO INVERT A MATRIX
  353. C
  354. C  B - ORIGINAL N-BY-N MATRIX, REPLACED BY INVERSE
  355. C  N - ORDER OF MATRIX
  356. C  NR- NUMBER OF ROWS B IS DIMENSIONED
  357. C  D - RESULTANT DETERMINANT
  358. C  A - WORK ARRAY, LENGTH N*N
  359. C  L,M - WORK ARRAYS, LENGTH N
  360. C
  361.     DIMENSION B(NR,1),A(1),L(1),M(1)
  362. C
  363. C CONVERT SQUARE ARRAY B TO LINEAR ARRAY A
  364. C
  365.     K=0
  366.     DO 5 J=1,N
  367.     DO 5 I=1,N
  368.     K=K+1
  369. 5    A(K)=B(I,J)
  370. C
  371. C SEARCH FOR LARGEST ELEMENT
  372. C
  373.     D=1.0
  374.     NK=-N
  375.     DO 80 K=1,N
  376.     NK=NK+N
  377.     L(K)=K
  378.     M(K)=K
  379.     KK=NK+K
  380.     BIGA=A(KK)
  381.     DO 20 J=K,N
  382.     IZ=N*(J-1)
  383.     DO 20 I=K,N
  384.     IJ=IZ+I
  385. 10    IF(ABS(BIGA)- ABS(A(IJ))) 15,20,20
  386. 15    BIGA=A(IJ)
  387.     L(K)=I
  388.     M(K)=J
  389. 20    CONTINUE
  390. C
  391. C INTERCHANGE ROWS
  392. C
  393.     J=L(K)
  394.     IF(J-K) 35,35,25
  395. 25    KI=K-N
  396.     DO 30 I=1,N
  397.     KI=KI+N
  398.     HOLD=-A(KI)
  399.     JI=KI-K+J
  400.     A(KI)=A(JI)
  401. 30    A(JI) =HOLD
  402. C
  403. C INTERCHANGE COLUMNS
  404. C
  405. 35    I=M(K)
  406.     IF(I-K) 45,45,38
  407. 38    JP=N*(I-1)
  408.     DO 40 J=1,N
  409.     JK=NK+J
  410.     JI=JP+J
  411.     HOLD=-A(JK)
  412.     A(JK)=A(JI)
  413. 40    A(JI) =HOLD
  414. C
  415. C DIVIDE COLUMN BY MINUS PIVOT
  416. C PIVOT IS CONTAINED IN BIGA
  417. C
  418. 45    IF(BIGA) 48,46,48
  419. 46    D=0.0
  420.     RETURN
  421. 48    DO 55 I=1,N
  422.     IF(I-K) 50,55,50
  423. 50    IK=NK+I
  424.     A(IK)=A(IK)/(-BIGA)
  425. 55    CONTINUE
  426. C
  427. C REDUCE MATRIX
  428. C
  429.     DO 65 I=1,N
  430.     IK=NK+I
  431.     IJ=I-N
  432.     DO 65 J=1,N
  433.     IJ=IJ+N
  434.     IF(I-K) 60,65,60
  435. 60    IF(J-K) 62,65,62
  436. 62    KJ=IJ-I+K
  437.     A(IJ)=A(IK)*A(KJ)+A(IJ)
  438. 65    CONTINUE
  439. C
  440. C DIVIDE ROW BY PIVOT
  441. C
  442.     KJ=K-N
  443.     DO 75 J=1,N
  444.     KJ=KJ+N
  445.     IF(J-K) 70,75,70
  446. 70    A(KJ)=A(KJ)/BIGA
  447. 75    CONTINUE
  448. C
  449. C PRODUCT OF PIVOTS
  450. C
  451. C *** DONT USE D
  452. C    D=D*BIGA
  453. C
  454. C REPLACE PIVOT BY RECIPROCAL
  455. C
  456.     A(KK)=1.0/BIGA
  457. 80    CONTINUE
  458. C
  459. C FINAL ROW AND COLUMN INTERCHANGE
  460. C
  461.     K=N
  462. 100    K=K-1
  463.     IF(K) 150,150,105
  464. 105    I=L(K)
  465.     IF(I-K) 120,120,108
  466. 108    JQ=N*(K-1)
  467.     JR=N*(I-1)
  468.     DO 110 J=1,N
  469.     JK=JQ+J
  470.     HOLD=A(JK)
  471.     JI=JR+J
  472.     A(JK)=-A(JI)
  473. 110    A(JI) =HOLD
  474. 120    J=M(K)
  475.     IF(J-K) 100,100,125
  476. 125    KI=K-N
  477.     DO 130 I=1,N
  478.     KI=KI+N
  479.     HOLD=A(KI)
  480.     JI=KI-K+J
  481.     A(KI)=-A(JI)
  482. 130    A(JI) =HOLD
  483.     GOTO 100
  484. C
  485. C CONVERT INVERTED MATRIX A TO N-BY-N MATRIX B
  486. C
  487. 150    K=0
  488.     DO 160 J=1,N
  489.     DO 160 I=1,N
  490.     K=K+1
  491. 160    B(I,J)=A(K)
  492.     RETURN
  493.     END
  494.     SUBROUTINE LINFIT(X,Y,YC,R,N,A,B,COR,SA,SB)
  495. C LINEAR FIT Y=A + B*X, CORR. COEFF
  496. C STD ERROR ON A AND B
  497. C
  498.     REAL X(1),Y(1),YC(1),R(1)
  499. C
  500.     SUMX  = 0
  501.     SUMY  = 0
  502.     SUMXY = 0
  503.     SUMY2 = 0
  504.     SUMX2 = 0
  505.     DO 10 I=1,N
  506.     XI = X(I)
  507.     YI = Y(I)
  508.     SUMX = SUMX+XI
  509.     SUMY = SUMY+YI
  510.     SUMXY = SUMXY + XI*YI
  511.     SUMX2 = SUMX2 + XI*XI
  512. 10    SUMY2 = SUMY2 + YI*YI
  513.     DENOM = SUMX2 - SUMX*SUMX/N
  514.     C = SUMXY - SUMX*SUMY/N
  515.     B = C/DENOM
  516.     A = (SUMY - B*SUMX)/N
  517.     COR = C / SQRT(DENOM * (SUMY2 - SUMY*SUMY / N))
  518.     SEE = SQRT((SUMY2 - A*SUMY - B*SUMXY) / (N-2))
  519.     SB = SEE / SQRT(DENOM)
  520.     SA = SB * SQRT(SUMX2/N)
  521.     DO 20 I=1,N
  522.     YC(I) = A + B * X(I)
  523.     R(I) = Y(I) - YC(I)
  524. 20    CONTINUE
  525.     RETURN
  526.     END
  527.