C @(#)mnag.for	17.1.1.1 (ES0-DMD) 01/25/02 17:19:00
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
      SUBROUTINE E04GAF(M, N, X, F, SS, E, MODE, D, W, IW, FUNCT,LSQ,
     * MONIT, IPRINT, MAXFUN, IFAIL)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENT  subroutine  E04GAF		modified by A. Lauberts  11.10.1985
C
C     MARK 4 RE-ISSUE. NAG COPYRIGHT 1974
C     MARK 5C REVISED
C     MARK 6B REVISED  IER-118 (SEP 1978)
C---------------------------------------------------------------------------
      INTEGER PAR1, PAR2, PAR3, PAR4, PAR5, PAR6, P01AAF, N, I, M,
     *MAXFUN, IFAIL, IW, MODE, NN, NWA, NWV, NWT, NWU, NWS, NWW,IERR,
     * IR, IT, NA, NWTI, IPRINT, NB, NWUI, JJ, NWWI, NWSI,NWVI, J, II,
     * NWTJ, ND, NDJ
      DOUBLE PRECISION ACC, RHO, SIG, Q, QC, SS, Z, VW, DQ, SSP,DS,
     * Y, TR, X(N), F(M), W(IW), D(N), E(N), X02AAF
      EXTERNAL   FUNCT, LSQ, MONIT
C
      LOGICAL IFL, CHECK
C
C     E04GAF IS BASED ON HARWELL SUBROUTINE VA07A (SEE 'A MODIFIED
C     MARQU
C     SUBROUTINE FOR NON-LINEAR LEAST SQUARES' BY R. FLETCHER,
C     HARWELL R
C     AERE-R6799, H.M.S.O. (MAY,1971)). IT IS NOT QUITE A
C     STRAIGHTFORWAR
C     TRANSLATION, AS CERTAIN CHANGRS HAVE BEEN MADE IN ORDER TO
C     CONFORM
C     CONVENTIONS, AND TO USE OTHER NAG ROUTINES.
C     ALL PRINTING HAS BEEN REMOVED. IT IS ASSUMED THAT THE USER WI
C     SUPPLY A PROCEDURE MONIT WHICH WILL BE CALLED ONCE EVERY
C     IPRINT IT
C     - IF IPRINT <=0, MONIT IS NOT CALLED. INSTEAD OF PRINTING
C     ERROR MA
C     E04GAF SETS THE ERROR INDICATOR IFAIL.
C     E04GAZ IS USED FOR THE FORMATION OF INNER PRODUCTS , AND F03A
C     F04AGF ARE USED TO SOLVE THE MARQUARDT EQUATIONS.
C     SO E0JGAF ASSUMES THAT LSQ SETS THE UPPER TRIANGLE OF A, UNLIK
C     WHICH ASSUMES THAT THE LOWER TRIANGLE IS GIVEN. THUS ,
C     THROUGHOUT
C     A REFERENCE TO A(I,J), SAY, CORRESPONDS TO A(J,I) IN THE
C     ORIGINAL.
C     THERE ARE NO MACHINE DEPENDENT CONSTANTS IN THE BODY OF E04GA
C     HOWEVER, ITS BEHAVIOUR MAY BE AFFECTED BY THE ACCURACY OF THE
C     MACH
C     (AND TO A LESSER EXTENT BY THE WAY IN WHICH ARITHMETIC
C     EXPRESSIONS
C     COMPILED AND EXECUTED).
C     W IS USED TO STORE MOST OF THE INTERNAL ARRAYS OF VA07A
C--------------------------------------------------------------------
C
      IF (N.LT.1) GO TO 40
      ACC = X02AAF(1.0D0)
      CHECK = .FALSE.
      DO 20 I=1,N
         CHECK = CHECK .OR. E(I).LT.ACC
         IF (MODE.NE.3) GO TO 20
         CHECK = CHECK .OR. D(I).LT.ACC
   20 CONTINUE
      IF (M.GE.1 .AND. MAXFUN.GE.1 .AND. (IFAIL.EQ.0 .OR.IFAIL.EQ.1) .
     *AND. IW.GE.N*(N+4)+M .AND. .NOT. CHECK) GO TO 60
   40 IFAIL = P01AAF(IFAIL,1)
      RETURN
   60 IF (MODE.NE.1 .AND. MODE.NE.2 .AND. MODE.NE.3) MODE = 1
      NN = N*N
      NWA = 0
      NWV = NN
      NWT = NWV + N
      NWU = NWT + N
      NWS = NWU + N
      NWW = NWS + N
      IERR = 0
      RHO = .25D0
      SIG = .75D0
      Q = 0.D0
      QC = 1.D0
      IFL = .FALSE.
      CALL FUNCT(M, N, X, F, IFL)
      IF (.NOT.IFL) GO TO 80
      IERR = 3
      GO TO 860
   80 CONTINUE
      IR = 1
      CALL E04GAZ(F, M, 1, 2, F, M, 1, 2, 0.0D0, SS, M, 4)
      CALL LSQ(M, N, X, F, W(NWA+1), W(NWV+1))
      IT = 0
      IF (MODE.EQ.3) GO TO 120
      IF (MODE.EQ.2) GO TO 160
      NA = NWA - N
      DO 100 I=1,N
         NA = NA + N + 1
         D(I) = W(NA)
         IF (D(I).LE.0.D0) D(I) = 1.D0
  100 CONTINUE
  120 CONTINUE
      DO 140 I=1,N
         NWTI = NWT + I
         W(NWTI) = DSQRT(D(I))
  140 CONTINUE
      GO TO 200
  160 CONTINUE
      DO 180 I=1,N
         NWTI = NWT + I
         W(NWTI) = 1.D0
         D(I) = 1.D0
  180 CONTINUE
  200 CONTINUE
      IF (IPRINT.LE.0) GO TO 220
      IF (MOD(IT,IPRINT).NE.0) GO TO 220
      CALL MONIT(M, N, X, F, SS, W(NWV+1), IR)
  220 CONTINUE
      IT = IT + 1
      NB = NWA - N
      DO 240 I=1,N
         NB = NB + N + 1
         NWUI = NWU + I
         W(NWUI) = W(NB)
         IF (W(NWUI).GE.0.0D0) GO TO 240
         IERR = 4
         GO TO 800
  240 CONTINUE
  260 CONTINUE
      NA = NWA - N
      DO 280 I=1,N
         NA = NA + N + 1
         NWUI = NWU + I
         W(NA) = W(NWUI) + Q*D(I)
  280 CONTINUE
      JJ = 1
      CALL F03AEF(N, W(NWA+1), N, W(NWW+1), Z, I, JJ)
C     W NOW CONTAINS THE RECIPROCALS OF THE DIAGONAL ELEMENTS OF
C     THE CHO
C     FACTORS.
      IF (JJ.EQ.0) GO TO 300
      Q = 2.D0*Q
      IF (Q.EQ.0.D0) Q = 1.D0
      GO TO 260
  300 CONTINUE
      CALL F04AGF(N, 1, W(NWA+1), N, W(NWW+1), W(NWV+1), N,W(NWS+1), N)
C     THE DIAGONAL ELEMENTS OF A ARE NOW THE RECIPROCALS OF THE
C     THOSE TH
C     WOULD BE PRODUCED BY THE ORIGINAL VA07A
      NA = NWA - N
      DO 320 I=1,N
         NA = NA + N + 1
         NWWI = NWW + I
         NWSI = NWS + I
         W(NA) = W(NWWI)
         W(NWWI) = W(NWSI)
  320 CONTINUE
      CALL E04GAZ(W(NWV+1), N, 1, 2, W(NWW+1), M, 1, 2, 0.0D0, VW, N,4)
      IF (VW.GT.0.D0) GO TO 340
      IERR = 5
      GO TO 800
  340 DO 360 I=1,N
         NWUI = NWU + I
         NWWI = NWW + I
         Z = W(NWUI)*W(NWWI)
         PAR1 = (I-1)*N + 1
         PAR2 = PAR1 + 1
         PAR3 = I*(N+1)
         PAR4 = PAR3 + N
         PAR5 = I + 1
         PAR6 = I + 2
         IF (I.GT.1) CALL E04GAZ(W(NWA+1), NN, PAR1, PAR2,W(NWW+1), M,
     *    1, 2, Z, Z, I-1, 4)
         IF (I.LT.N) CALL E04GAZ(W(NWA+1), NN, PAR3, PAR4,W(NWW+1), M,
     *    PAR5, PAR6, Z, Z, N-I, 4)
         NWSI = NWS + I
         NWVI = NWV + I
         W(NWSI) = 2.0D0*W(NWVI) - Z
  360 CONTINUE
      DQ = 0.D0
      CALL E04GAZ(W(NWS+1), N, 1, 2, W(NWW+1), M, 1, 2, 0.0D0, DQ, N,4)
      J = 0
      DO 380 I=1,N
         NWWI = NWW + I
         NWSI = NWS + I
         IF (DABS(W(NWWI)).GT.E(I)) J = 1
         W(NWSI) = X(I) - W(NWWI)
  380 CONTINUE
      IFL = .FALSE.
      CALL FUNCT(M, N, W(NWS+1), W(NWW+1), IFL)
      IF (IFL) GO TO 440
C 390 CONTINUE
      CONTINUE
      IR = IR + 1
      CALL E04GAZ(W(NWW+1), M, 1, 2, W(NWW+1), M, 1, 2, 0.0D0, SSP,M, 4)
      DS = SS - SSP
C     ONLY THE DIAGONAL ELEMENTS OF A HAVE BEEN CHANGED SINCE
C     F04AGF WAS
C     SO , IF AN EXIT IS MADE, ONLY THE DIAGONAL ELEMENTS NEED TO
C     BE RES
      IF (J.EQ.0) GO TO 740
C     NOTICE THAT E04GAF USES IFAIL TO DISTINGUISH BETWEEN
C     DIFFERENT TYP
C     EXIT.
      IF (DQ.GT.0.0D0) GO TO 400
      IERR = 5
      GO TO 740
  400 CONTINUE
      IF (IR.LT.MAXFUN) GO TO 420
      IERR = 2
      GO TO 740
  420 CONTINUE
      IF (DS.GE.RHO*DQ) GO TO 660
      Y = .5D0
      Z = 2.D0*VW - DS
      IF (Z.GT.0.D0) Y = VW/Z
      IF (Y.GT..5D0) Y = .5D0
      IF (Y.LT..1D0) Y = .1D0
      GO TO 460
  440 Y = .1D0
      DS = 0.D0
  460 CONTINUE
      IF (Q.NE.0.D0) GO TO 640
C     THE DIAGONAL ELEMENTS OF A ARE ALREADY THE RECIPROCALS OF
C     THOSE TH
C     BE PRODUCED BY VA07A. SO, 2 LINES HAVE BEEN OMMITTED HERE IN
C     THE
C     TRANSLATION .
      Y = 2.D0*Y
      NB = NWA + 1
      DO 500 I=2,N
         II = I - 1
         NA = NWA + I - N
         NB = NB + N + 1
         DO 480 J=1,II
            PAR1 = (J-1)*N + J
            PAR2 = PAR1 + 1
            PAR3 = (J-1)*N + I
            PAR4 = PAR3 + N
            CALL E04GAZ(W(NWA+1), NN, PAR1, PAR2, W(NWA+1), NN,PAR3,
     *       PAR4, 0.0D0, Z, I-J, 3)
            Z = -Z
            NA = NA + N
            W(NA) = Z*W(NB)
  480    CONTINUE
  500 CONTINUE
      DO 540 I=1,N
         NA = NWA + (I-1)*(N+1)
         DO 520 J=I,N
            PAR1 = (I-1)*N + J
            PAR2 = PAR1 + 1
            PAR3 = (J-1)*N + J
            PAR4 = PAR3 + 1
            CALL E04GAZ(W(NWA+1), NN, PAR1, PAR2, W(NWA+1), NN,PAR3,
     *       PAR4, 0.0D0, Z, N-J+1, 4)
            NA = NA + 1
            W(NA) = DABS(Z)
  520    CONTINUE
  540 CONTINUE
      Q = 0.D0
      TR = 0.D0
      NA = NWA - N
      DO 620 I=1,N
         NA = NA + N + 1
         TR = TR + W(NA)*D(I)
         Z = 0.D0
         NB = NWA + I - N
         DO 560 J=1,I
            NB = NB + N
            NWTJ = NWT + J
            Z = Z + W(NB)*W(NWTJ)
  560    CONTINUE
         IF (I.EQ.N) GO TO 600
         II = I + 1
         ND = NWA + (I-1)*N
         DO 580 J=II,N
            NWTJ = NWT + J
            NDJ = ND + J
            Z = Z + W(NDJ)*W(NWTJ)
  580    CONTINUE
  600    CONTINUE
         NWTI = NWT + I
         Z = Z*W(NWTI)
         IF (Z.GT.Q) Q = Z
  620 CONTINUE
      IF (TR.LT.Q) Q = TR
      Q = 1.D0/Q
      QC = Q
  640 CONTINUE
      Q = Q/Y
      IF (DS) 260, 260, 680
  660 CONTINUE
      IF (DS.LE.SIG*DQ) GO TO 680
      Q = Q*.5D0
      IF (Q.LT.QC) Q = 0.D0
  680 CONTINUE
      SS = SSP
      DO 700 I=1,N
         NWSI = NWS + I
         X(I) = W(NWSI)
  700 CONTINUE
      DO 720 I=1,M
         NWWI = NWW + I
         F(I) = W(NWWI)
  720 CONTINUE
      CALL LSQ(M, N, X, F, W(NWA+1), W(NWV+1))
      GO TO 200
  740 CONTINUE
      IF (DS.LE.0.D0) GO TO 800
      SS = SSP
      DO 760 I=1,N
         NWSI = NWS + I
         X(I) = W(NWSI)
  760 CONTINUE
      DO 780 I=1,M
         NWWI = NWW + I
         F(I) = W(NWWI)
  780 CONTINUE
  800 CONTINUE
      NA = NWA - N
C     RESTORE THE DIAGONAL ELEMENTS OF A.
      DO 820 I=1,N
         NA = NA + N + 1
         NWUI = NWU + I
         W(NA) = W(NWUI)
  820 CONTINUE
      IF (IERR) 840, 840, 860
  840 IFAIL = 0
      GO TO 880
  860 IFAIL = P01AAF(IFAIL,IERR)
  880 RETURN
      END


      SUBROUTINE E04GAZ(A, M1, I, J, B, M2, K, L, X, SUM, N, IFLAG)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENT subroutine E04GAZ
C
C     MARK 2 RELEASE. NAG COPYRIGHT 1972
C---------------------------------------------------------------
      INTEGER M, N, MM, I, J, K, L, IFLAG, M1, M2
C
      DOUBLE PRECISION SX, AX, BX, XX
      DOUBLE PRECISION X, SUM, A(M1), B(M2)
      SX = 0.0D0
      DO 20 M=1,N
         MM = I + (M-1)*(J-I)
         AX = A(MM)
         MM = K + (M-1)*(L-K)
         BX = B(MM)
         SX = SX + AX*BX
   20 CONTINUE
      XX = X
      M = IFLAG + 1
      GO TO (40, 60, 80, 100, 40, 60, 80, 100), M
   40 SX = SX + XX
      GO TO 120
   60 SX = XX - SX
      GO TO 120
   80 SX = SX - XX
      GO TO 120
  100 SX = -SX - XX
  120 SUM = SX
      RETURN
      END


      SUBROUTINE F03AEF(N, A, IA, P, D1, ID, IFAIL)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENT subroutine F03AEF
C
C     MARK 2 RELEASE. NAG COPYRIGHT 1972
C     MARK 3 REVISED.
C     MARK 4.5 REVISED
C
C     CHOLDET1
C     THE UPPER TRIANGLE OF A POSITIVE DEFINITE SYMMETRIC MATRIX,
C     A, IS STORED IN THE UPPER TRIANGLE OF AN N*N ARRAY A(I,J),
C     I=1,N, J=1,N. THE CHOLESKY DECOMPOSITION A = LU, WHERE
C     U IS THE TRANSPOSE OF L, IS PERFORMED AND STORED IN THE
C     REMAINDER OF THE ARRAY A EXCEPT FOR THE RECIPROCALS OF THE
C     DIAGONAL ELEMENTS WHICH ARE STORED IN P(I), I = 1,N,
C     INSTEAD OF THE ELEMENTS THEMSELVES. A IS RETAINED SO THAT
C     THE SOLUTION OBTAINED CAN BE SUBSEQUENTLY IMPROVED. THE
C     DETERMINANT, D1 * 2.0**ID, OF A IS ALSO COMPUTED. THE
C     SUBROUTINE
C     WILL FAIL IF A, MODIFIED BY THE ROUNDING ERRORS, IS NOT
C     POSITIVE DEFINITE.
C     1ST DECEMBER 1971
C----------------------------------------------------------------------
      INTEGER ISAVE, IFAIL, IFAIL1, ID, I, N, ISTART, J, IA, P01AAF
C
      DOUBLE PRECISION D1, X, D2, P(N), A(IA,N)
      ISAVE = IFAIL
      IFAIL1 = 0
      D1 = 1.0D0
      ID = 0
      DO 120 I=1,N
         ISTART = I
         DO 100 J=ISTART,N
            CALL X03AAF(A(I,1), IA*N+1-I, A(J,1), IA*N+1-J, I-1,IA, IA,
     *       -A(I,J), 0.0D0, X, D2, .FALSE., IFAIL1)
            X = -X
            IF (J.NE.I) GO TO 80
            D1 = D1*X
            IF (X.LE.0.0D0) GO TO 140
   20       IF (DABS(D1).LT.1.0D0) GO TO 40
            D1 = D1*0.0625D0
            ID = ID + 4
            GO TO 20
   40       IF (DABS(D1).GE.0.0625D0) GO TO 60
            D1 = D1*16.0D0
            ID = ID - 4
            GO TO 40
   60       P(I) = 1.0D0/DSQRT(X)
            GO TO 100
   80       A(J,I) = X*P(I)
  100    CONTINUE
  120 CONTINUE
      IFAIL = 0
      RETURN
  140 IFAIL = P01AAF(ISAVE,1)
      RETURN
      END


      SUBROUTINE F04AGF(N, IR, A, IA, P, B, IB, X, IX)
C.IDENT subroutine F04AGF
C
C     MARK 2 RELEASE. NAG COPYRIGHT 1972
C     MARK 4.5 REVISED
C     CHOLSOL1
C     SOLVES AX=B, WHERE A IS A POSITIVE DEFINITE SYMMETRIC MATRIX
C     AND B IS AN N*IR MATRIX OF IR RIGHT-HAND SIDES. THE
C     SUBROUTINE
C     F04AGF MUST BY PRECEDED BY F03AEF IN WHICH L IS
C     PRODUCED IN A(I,J) AND P(I), FROM A. AX=B IS SOLVED IN TWO
C     STEPS, LY=B AND UX=Y, AND X IS OVERWRITTEN ON Y.
C     1ST AUGUST 1971
C-----------------------------------------------------------------------
      INTEGER IFAIL1, J, IR, I, N, I1, IA, IX, II, IB
      DOUBLE PRECISION D1, D2, A(IA,N), P(N), B(IB,IR), X(IX,IR)
      IFAIL1 = 0
      DO 60 J=1,IR
C     SOLUTION OF LY= B
         DO 20 I=1,N
            I1 = I - 1
            CALL X03AAF(A(I,1), IA*N-I+1, X(1,J), (IR-J+1)*IX, I1,IA, 1,
     *       B(I,J), 0.0D0, D1, D2, .FALSE., IFAIL1)
            X(I,J) = -P(I)*D1
   20    CONTINUE
C     SOLUTION OF UX= Y
         X(N,J) = -P(N)*X(N,J)
         IF (N.EQ.1) GO TO 60
         DO 40 II=2,N
            I = N - II + 1
            I1 = I + 1
            CALL X03AAF(A(I1,I), (N-I+1)*IA-I, X(I1,J),(IR-J+1)*IX-I, N-
     *      I, 1, 1, X(I,J), 0.0D0, D1, D2, .FALSE.,IFAIL1)
            X(I,J) = -P(I)*D1
   40    CONTINUE
   60 CONTINUE
      RETURN
      END


      INTEGER FUNCTION P01AAF(IFAIL, ERROR)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENT integer function P01AAF
C
C     MARK 1 RELEASE.  NAG COPYRIGHT 1971
C     MARK 3 REVISED
C     MARK 4A REVISED, IER-45
C     MARK 4.5 REVISED
C     MARK 7 REVISED (DEC 1978)
C     RETURNS THE VALUE OF ERROR, REVISED (AL, NOV 1982)
C------------------------------------------------------------------
      INTEGER ERROR, IFAIL, NOUT
C
C     TEST IF NO ERROR DETECTED
      IF (ERROR.EQ.0) GO TO 20
C     DETERMINE OUTPUT UNIT FOR MESSAGE
      CALL X04AAF (0,NOUT)
C     TEST FOR SOFT FAILURE
      IF (MOD(IFAIL,10).EQ.1) GO TO 10
C     HARD FAILURE
C     WRITE (NOUT,99999) ERROR
C     STOPPING MECHANISM NOT APPLIED (REVISED/AL)
C     CALL P01AAZ
      P01AAF=0                                     ! inserted RHW
      RETURN
C     SOFT FAIL
C     TEST IF ERROR MESSAGES SUPPRESSED
   10 IF (MOD(IFAIL/10,10).EQ.0) GO TO 20
C     WRITE (NOUT,99999) ERROR
   20 P01AAF = ERROR
      RETURN
C99999FORMAT (1H0, 38HERROR_DETECTED_BY_LIBRARY_ROUTINE_, A8,
C    * 11H - IFAIL = , I5//)
      END


      DOUBLE PRECISION FUNCTION X02AAF(X)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENT double precision function X02AAF
C
C     NAG COPYRIGHT 1975
C     MARK 4.5 RELEASE
C-----------------------------------------------------------------
      DOUBLE PRECISION X
C+
C     * EPS *
C     RETURNS THE VALUE EPS WHERE EPS IS THE SMALLEST
C     POSITIVE
C     NUMBER SUCH THAT 1.0 + EPS > 1.0
C     THE X PARAMETER IS NOT USED
C     FOR DEC VAX 11/780
C     X02AAF = 2.0**(-56.0)
C-
      X02AAF = 2.0D0**(-55)
      RETURN
      END


      SUBROUTINE X03AAF(A, ISIZEA, B, ISIZEB, N, ISTEPA, ISTEPB,C1, C2,
     * D1, D2, SW, IFAIL)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENT subroutine X03AAF
C
C     NAG COPYRIGHT 1975
C     MARK 4.5 RELEASE
C     MARK 6 REVISED
C
C     CALCULATES THE VALUE OF A SCALAR PRODUCT USING BASIC
C     OR ADDITIONAL PRECISION AND ADDS IT TO A BASIC OR ADDITIONAL
C     PRECISION INITIAL VALUE.
C--------------------------------------------------------------------
C.... DELETE THIS LINE WHEN ADDITIONAL PRECISION HAS BEEN DEALT WITH
      INTEGER P01AAF, ISAVE, ISIZEA, ISIZEB, ISTEPA, ISTEPB, IFAIL,IS,
     * IT, N, I
C
      DOUBLE PRECISION SUM
      DOUBLE PRECISION A(ISIZEA), B(ISIZEB), C1, C2, D1, D2, X
      LOGICAL SW
      ISAVE = IFAIL
      IFAIL = 0
      IF (ISTEPA.GT.0 .AND. ISTEPB.GT.0) GO TO 20
      IFAIL = P01AAF(ISAVE,1)
      RETURN
   20 IS = 1 - ISTEPA
      IT = 1 - ISTEPB
      IF (SW) GO TO 80
      X = 0.0D0
      IF (N.LT.1) GO TO 60
      DO 40 I=1,N
         IS = IS + ISTEPA
         IT = IT + ISTEPB
         X = X + A(IS)*B(IT)
   40 CONTINUE
   60 D1 = X + (C1+C2)
      D2 = 0.0D0
      RETURN
C
   80 SUM = 0.0D0
      IF (N.LT.1) GO TO 120
      DO 100 I=1,N
         IS = IS + ISTEPA
         IT = IT + ISTEPB
C
C        SUM = SUM + DBLE(A(IS))*B(IT)
         SUM = SUM + (A(IS))*B(IT)
  100 CONTINUE
C
C 120 SUM = SUM + (DBLE(C1)+C2)
  120 SUM = SUM + ((C1)+C2)
C
C     D1 = SUM + SUM - DBLE(SNGL(SUM))
      D1 = SUM + SUM - ((SUM))
C     THE LAST STATEMENT ASSUMES THAT THE MACHINE SIMPLY
C     TRUNCATES WHEN ASSIGNING A DOUBLE PRECISION QUANTITY
C     TO A SINGLE PRECISION VARIABLE. IF INSTEAD THE MACHINE
C     ROUNDS, REPLACE THE LAST STATEMENT BY
C     D1 = SUM
      D2 = SUM - D1
      RETURN
      END


      SUBROUTINE X04AAF(I,NERR)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENT subroutine X04AAF
C
C     MARK 7 RELEASE. NAG COPYRIGHT 1978
C     MARK 7C REVISED IER-190 (MAY 1979)
C     IF I = 0, SETS NERR TO CURRENT ERROR MESSAGE UNIT NUMBER
C     (STORED IN NERR1).
C     IF I = 1, CHANGES CURRENT ERROR MESSAGE UNIT NUMBER TO
C     VALUE SPECIFIED BY NERR.
C
C     *** NOTE ***
C     THIS ROUTINE ASSUMES THAT THE VALUE OF NERR1 IS SAVED
C     BETWEEN CALLS.  IN SOME IMPLEMENTATIONS IT MAY BE
C     NECESSARY TO STORE NERR1 IN A LABELLED COMMON
C     BLOCK /AX04AA/ TO ACHIEVE THIS.
C-------------------------------------------------------------------
C     .. SCALAR ARGUMENTS ..
      INTEGER I, NERR
C     ..
C     .. LOCAL SCALARS ..
      INTEGER NERR1
C     ..
      DATA NERR1 /6/
      IF (I.EQ.0) NERR = NERR1
      IF (I.EQ.1) NERR1 = NERR
      RETURN
      END