MINICERN/mathlib/gen/e/dlsqp2.F

   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