]>
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 DLSQPM(N,X,Y,M,A,SD,IFAIL) | |
11 | #if !defined(CERNLIB_DOUBLE) | |
12 | CHARACTER*6 NAME | |
13 | NAME = 'DLSQPM' | |
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 | PARAMETER (IDIM = 21, R0 = 0) | |
23 | ||
24 | DIMENSION X(*),Y(*),A(0:*),B(IDIM,IDIM),XY(0:IDIM) | |
25 | ||
26 | M1=M+1 | |
27 | IFAIL=0 | |
28 | IF(N .LE. 1 .OR. M .LT. 0 .OR. M1 .GT. IDIM .OR. M1 .GT. N) THEN | |
29 | IFAIL=1 | |
30 | ELSEIF(M .EQ. 0) THEN | |
31 | XY(0)=DVSUM(N,Y(1),Y(2)) | |
32 | A(0)=XY(0)/N | |
33 | SYY=DVMPY(N,Y(1),Y(2),Y(1),Y(2)) | |
34 | ELSE | |
35 | DO 11 J = 1,M1 | |
36 | A(J-1)=0 | |
37 | B(J,1)=0 | |
38 | B(M1,J)=0 | |
39 | 11 CONTINUE | |
40 | B(1,1)=N | |
41 | SYY=0 | |
42 | DO 4 K = 1,N | |
43 | XK=X(K) | |
44 | YK=Y(K) | |
45 | SYY=SYY+Y(K)**2 | |
46 | P=1 | |
47 | A(0)=A(0)+YK | |
48 | DO 2 J = 2,M1 | |
49 | P=P*XK | |
50 | B(J,1)=B(J,1)+P | |
51 | A(J-1)=A(J-1)+P*YK | |
52 | 2 CONTINUE | |
53 | DO 3 J = 2,M1 | |
54 | P=P*XK | |
55 | B(M1,J)=B(M1,J)+P | |
56 | 3 CONTINUE | |
57 | 4 CONTINUE | |
58 | DO 5 I = 2,M | |
59 | DO 5 K = I,M1 | |
60 | B(K-1,I)=B(K,I-1) | |
61 | 5 CONTINUE | |
62 | DO 6 J = 0,M | |
63 | XY(J)=A(J) | |
64 | 6 CONTINUE | |
65 | CALL DSEQN(M1,B,IDIM,IFAIL,1,A) | |
66 | ENDIF | |
67 | SD=0 | |
68 | IF(IFAIL .EQ. 0) THEN | |
69 | IF(N .GT. M1) THEN | |
70 | SD=SYY | |
71 | DO 7 J = 0,M | |
72 | SD=SD-A(J)*XY(J) | |
73 | 7 CONTINUE | |
74 | SD=SQRT(MAX(R0,SD)/(N-M1)) | |
75 | ENDIF | |
76 | ELSE | |
77 | CALL DVSET(M1,R0,A(0),A(1)) | |
78 | M=0 | |
79 | ENDIF | |
80 | RETURN | |
81 | END | |
82 | #endif |