]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:02:24 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | SUBROUTINE DLSQP2(N,X,Y,A0,A1,A2,SD,IFAIL) | |
11 | #if !defined(CERNLIB_DOUBLE) | |
12 | CHARACTER*6 NAME | |
13 | NAME = 'DLSQP2' | |
14 | CALL MTLPRT(NAME,'E201', | |
15 | +'not available on this machine - see documentation') | |
16 | RETURN | |
17 | END | |
18 | #endif | |
19 | #if defined(CERNLIB_DOUBLE) | |
20 | IMPLICIT DOUBLE PRECISION (A-H,O-Z) | |
21 | ||
22 | DIMENSION X(*),Y(*) | |
23 | ||
24 | PARAMETER (R0 = 0) | |
25 | ||
26 | A0=0 | |
27 | A1=0 | |
28 | A2=0 | |
29 | SD=0 | |
30 | IF(N .LE. 2) THEN | |
31 | IFAIL=1 | |
32 | ELSE | |
33 | FN=N | |
34 | XM=0 | |
35 | DO 1 K = 1,N | |
36 | XM=XM+X(K) | |
37 | 1 CONTINUE | |
38 | XM=XM/FN | |
39 | SX=0 | |
40 | SXX=0 | |
41 | SXXX=0 | |
42 | SXXXX=0 | |
43 | SY=0 | |
44 | SYY=0 | |
45 | SXY=0 | |
46 | SXXY=0 | |
47 | DO 2 K = 1,N | |
48 | XK=X(K)-XM | |
49 | YK=Y(K) | |
50 | XK2=XK**2 | |
51 | SX=SX+XK | |
52 | SXX=SXX+XK2 | |
53 | SXXX=SXXX+XK2*XK | |
54 | SXXXX=SXXXX+XK2**2 | |
55 | SY=SY+YK | |
56 | SYY=SYY+YK**2 | |
57 | SXY=SXY+XK*YK | |
58 | SXXY=SXXY+XK2*YK | |
59 | 2 CONTINUE | |
60 | DET=(FN*SXXXX-SXX**2)*SXX-FN*SXXX**2 | |
61 | IF(DET .GT. 0) THEN | |
62 | A2=(SXX*(FN*SXXY-SXX*SY)-FN*SXXX*SXY)/DET | |
63 | A1=(SXY-SXXX*A2)/SXX | |
64 | A0=(SY-SXX*A2)/FN | |
65 | IFAIL=0 | |
66 | ELSE | |
67 | IFAIL=-1 | |
68 | ENDIF | |
69 | ENDIF | |
70 | IF(IFAIL .EQ. 0 .AND. N .GT. 3) | |
71 | 1 SD=SQRT(MAX(R0,SYY-A0*SY-A1*SXY-A2*SXXY)/(N-3)) | |
72 | A0=A0+XM*(XM*A2-A1) | |
73 | A1=A1-2*XM*A2 | |
74 | RETURN | |
75 | END | |
76 | #endif |