outlets/simple.f

      PROGRAM OUTLETS
      IMPLICIT NONE
      INTEGER I, SUCC, PRED
      INTEGER N
      PARAMETER (N = 9)
      DOUBLE PRECISION L, XP(9), C, P(9), S(9), PRF(9), PRFTOT(9)
      INTEGER F(9)

      C = 1.0D0 
      L = 10.D0
      DO 3, I = 1, 9
  3   XP(I) = L / 9.D0
      F(1) = 1
      F(2) = 2
      F(3) = 3
      F(4) = 3
      F(5) = 2
      F(6) = 1
      F(7) = 2
      F(8) = 1
      F(9) = 3
     
      CALL OUTLET(N, F, XP, C, P, S, PRF, PRFTOT)

      END

      SUBROUTINE OUTLET(N, F, XP, C, P, S, PRF, PRFTOT)
      IMPLICIT NONE
      INTEGER N, F(N)  
      INTEGER IPIV(N), WORK(N), LWORK, INFO, PRED
      DOUBLE PRECISION II(N)
      DOUBLE PRECISION XP(N), C, P(N), S(N), PRF(N), PRFTOT(N)
      DOUBLE PRECISION A(N,N), B(N)
      DOUBLE PRECISION ZERO, ONE, TWO
      INTEGER I, J

      ONE = 1.0D0
      TWO = 2.0D0

      LWORK = N

      II = ONE
      DO 1, I = 1, N - 1
   1  IF(F(I).EQ.F(I+1)) II(I) = TWO
      
      A = ZERO 
      DO 5, I = 1, N
      A(I,I) = TWO*( XP(PRED(I,N)) + XP(I) ) / ( XP(PRED(I,N)) * XP(I) )
   5  B(I) = C * ( XP(PRED(I,N)) + XP(I) )
      DO 7, I = 1, N - 1
      A(I+1,I) = - II(I) / XP(I)
   7  A(I,I+1) = - II(I) / XP(I)
      A(1,N)   = - II(N) / XP(1)
      A(N,1)   = - II(N) / XP(1)

      DO 8, I = 1, N
   8  PRINT '(I5,10F7.1)', I, (A(I,J), J=1,N), B(I)

C     CALL DSYSV('U', N, 1, A, N, IPIV, B, N, WORK, LWORK, INFO)
      CALL DGESV(N,1,A,N,IPIV,B,N,INFO)

      PRINT*
      PRINT*, INFO

      DO 10, I = 1, N
      S(I) = (XP(PRED(I,N)) + XP(I))/TWO +
     &   (B(PRED(I,N))/XP(PRED(I,N))+B(PRED(I,N))/XP(PRED(I,N)))/(TWO*C)
  10  PRF(I) = (B(I)-C) * S(I)


      DO 30, I = 1, N
  30  PRINT '(2I5,3F9.3)', I, F(I), B(I), S(I), PRF(I)

      RETURN
      END




      FUNCTION SUCC(I,N)
      INTEGER SUCC, I, N

      IF (I.LT.N) SUCC = I + 1
      IF (I.EQ.N) SUCC = 1 

      RETURN
      END


      FUNCTION PRED(I,N)
      INTEGER PRED, I, N

      IF (I.EQ.1) PRED = N
      IF (I.GT.1) PRED = I - 1 

      RETURN
      END