!     Last change:  J     8 Oct 2000   11:40 am
!                                       *** copyright 2000 ***
! *** filename sum3.f95                 *** John F. Monahan **
!                                       **********************
      PROGRAM PSUM3
!  TEST OF ADDITION ALGORITHM THAT ADDS FROM 1 TO N BEST WAY
!                                                                       
      IMPLICIT NONE
      INTEGER, PARAMETER :: N = 32768   ! HOW MANY TO ADD
      REAL, DIMENSION(N) :: X,Y,Z
      REAL XX,YY,ZZ
      REAL, PARAMETER :: G = .5772157   ! EULER'S CONSTANT
      INTEGER I
      INTEGER IXX,IYY,IZZ,IX1,IY1,IZ1
!                   some compilers only allow integers to be
!                    written in z (internal) format
      EQUIVALENCE (XX,IXX),(YY,IYY),(ZZ,IZZ),(IX1,X(1)),    &
     &          (IY1,Y(1)),(IZ1,Z(1))
!
  21  FORMAT(' N =',I12,4X,'ANALYTIC SUM',E20.10,Z12/       &
     &       'SUM3',E20.10,Z12)
  25  FORMAT(/' FIRST EXAMPLE - ADD NUMBERS FROM 1 TO N')
  26  FORMAT(/' SECOND EXAMPLE - ADD SQUARES FROM 1 TO N')
  27  FORMAT(/' THIRD EXAMPLE - HARMONIC SUM')
!
      OPEN (UNIT=6, FILE='sum3.out')
!                                                                       
!   LOAD UP VECTOR OF VALUES FOR THREE EXAMPLES
!
      DO I = 1,N                                                      
      X(I) = I                                                          
      Y(I) = X(I)*X(I)                                                  
      Z(I) = 1.0/REAL(N+1-I)
      END DO  ! LOOP ON I
!
!   ADD THEM UP USING SUM3 AND COMPARE TO ANALYTICAL SUMMATION
! 
      WRITE(6,25)
      CALL SUM3(X,N)          ! FIRST PROBLEM
      XX = REAL(N)*REAL(N+1)/2.
      WRITE(6,21) N,XX,IXX,X(1),IX1                                     
      WRITE(6,26)
      CALL SUM3(Y,N)          ! SECOND PROBLEM
      YY = REAL(N)*REAL(N+1)*REAL(2*N+1)/6.
      WRITE(6,21) N,YY,IYY,Y(1),IY1                                     
      WRITE(6,27)
      CALL SUM3(Z,N)          ! THIRD PROBLEM
      ZZ = G + LOG(REAL(N))
      WRITE(6,21) N,ZZ,IZZ,Z(1),IZ1                                     
      STOP                                                              
      END PROGRAM PSUM3
      SUBROUTINE SUM3(X,N)                                              
!  ALGORITHM TO SUM AN ORDERED (NONDECREASING) LIST OF NUMBERS
! IN THE BEST POSSIBLE WAY - ALWAYS ADDING SMALLEST TOGETHER
      INTEGER , INTENT(IN) :: N
      REAL, DIMENSION( N ), INTENT(IN OUT) :: X
      INTEGER XS,XL,SS,SL,XID,SID
!
      IF( N .EQ. 1 ) THEN
                          RETURN
                     ENDIF
!                         START   XS = X START    XL = X LENGTH
!                                 SS = SUM START  SL = SUM LENGTH
      XS = 0                                                            
      SS = 0                                                            
      XL = N                                                            
      SL = 0
      CALL SUMX1X2
!                   MAIN LOOP IS DO WITH EXIT
      DO
!                   CASE BASED ON XID AND SID -- LENGTH OF LISTS
      SELECT CASE (XID+XID+SID)
          CASE (1)
!                   ONLY ONE ELEMENT LEFT -- WE'RE DONE
      X(1) = X(SS+1)                                                    
      EXIT
          CASE (2)
!                   XID = 0, SID=2  XLIST EMPTY
!                   RESET AFTER FINISHING LIST - ADD SUMLIST
      XS = SS                                                           
      XL = SL                                                           
      SS = 0                                                            
      SL = 0                                                            
      CALL SUMX1X2
          CASE (3)
!                   XID = 1, SID = 1
      CALL SUMX1S1
          CASE (4)
!                   XID = 1, SID = 2
!                   SPECIAL - 1 IN X LIST, MORE IN SUMLIST
      IF( X(XS+1) .GT. X(SS+2) ) THEN
                                      CALL SUMS1S2
                                 ELSE
                                      CALL SUMX1S1
                                 ENDIF
          CASE (5)
!                   XID = 2, SID = 1
!                   SPECIAL - 1 IN SUMLIST, MORE IN X LIST
      IF( X(SS+1) .GT. X(XS+2) ) THEN
                                      CALL SUMX1X2
                                 ELSE
                                      CALL SUMX1S1
                                 ENDIF
          CASE (6)
!                   XID = 2, SID = 2
!                   REGULAR SITUATION - BOTH LISTS LONG
      IF( X(XS+2) .LE. X(SS+1) ) THEN
                                      CALL SUMX1X2
           ELSEIF( X(XS+1) .LE. X(SS+2) ) THEN
                                               CALL SUMX1S1
                                          ELSE
                                               CALL SUMS1S2
                                          ENDIF
              CASE DEFAULT
              EXIT
      END SELECT
      END DO        ! END OF MAIN DO LOOP -- EXIT TO RETURN
      RETURN
      CONTAINS
      SUBROUTINE SUMS1S2
      X(SS+SL+1) = X(SS+1) + X(SS+2)
      SL = SL - 1                                                       
      SS = SS + 2                                                       
      SID = MIN(SL,2)                                                  
      END SUBROUTINE SUMS1S2

      SUBROUTINE SUMX1X2
      SL = SL +1                                                        
      X(SS+SL) = X(XS+1) + X(XS+2)                                      
      XS = XS + 2                                                       
      XL = XL - 2                                                       
      XID = MIN(XL,2)                                                  
      SID = MIN(SL,2)
      END SUBROUTINE SUMX1X2

      SUBROUTINE SUMX1S1
      SS = SS + 1
      XS = XS + 1                                                       
      X(SS+SL) = X(SS) + X(XS)                                          
      XL = XL - 1                                                       
      XID = MIN(XL,2)                                                  
      END SUBROUTINE SUMX1S1
      END SUBROUTINE SUM3
! *** end of filename sum3.f95          *********************