home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / villad / dfopr.f next >
Text File  |  1996-09-28  |  4KB  |  177 lines
  1.       SUBROUTINE DFOPR
  2.      +  (
  3.      +  ND, N, N0, N1, I, ID, DIF1, DIF2, DIF3, ROOT, VECT
  4.      +  )
  5.       INTEGER           ND, N, N0, N1, I, ID
  6.       DOUBLE PRECISION  DIF1(ND), DIF2(ND), DIF3(ND), ROOT(ND), VECT(ND)
  7. C
  8. C***********************************************************************
  9. C
  10. C     VILLADSEN AND MICHELSEN, PAGES 133-134, 419
  11. C
  12. C     INPUT PARAMETERS:
  13. C
  14. C       ND     : THE DIMENSION OF THE VECTORS DIF1, DIF2, DIF3, AND ROOT
  15. C
  16. C       N      : THE DEGREE OF THE JACOBI POLYNOMIAL, (i.e. THE NUMBER
  17. C                OF INTERIOR INTERPOLATION POINTS)
  18. C
  19. C       N0     : DETERMINES WHETHER X = 0 IS INCLUDED AS AN
  20. C                INTERPOLATION POINT
  21. C
  22. C                  N0 = 0  ==>  X = 0 IS NOT INCLUDED
  23. C                  N0 = 1  ==>  X = 0 IS INCLUDED
  24. C
  25. C       N1     : DETERMINES WHETHER X = 1 IS INCLUDED AS AN
  26. C                INTERPOLATION POINT
  27. C
  28. C                  N1 = 0  ==>  X = 1 IS NOT INCLUDED
  29. C                  N1 = 1  ==>  X = 1 IS INCLUDED
  30. C
  31. C       I      : THE INDEX OF THE NODE FOR WHICH THE WEIGHTS ARE TO BE
  32. C                CALCULATED
  33. C
  34. C       ID     : INDICATOR
  35. C
  36. C                  ID = 1  ==>  FIRST DERIVATIVE WEIGHTS ARE COMPUTED
  37. C                  ID = 2  ==>  SECOND DERIVATIVE WEIGHTS ARE COMPUTED
  38. C                  ID = 3  ==>  GAUSSIAN WEIGHTS ARE COMPUTED (IN THIS
  39. C                               CASE, THE VALUE OF I IS IRRELEVANT)
  40. C
  41. C     OUTPUT PARAMETERS:
  42. C
  43. C       DIF1   : ONE DIMENSIONAL VECTOR CONTAINING THE FIRST DERIVATIVE
  44. C                OF THE NODE POLYNOMIAL AT THE ZEROS
  45. C
  46. C       DIF2   : ONE DIMENSIONAL VECTOR CONTAINING THE SECOND DERIVATIVE
  47. C                OF THE NODE POLYNOMIAL AT THE ZEROS
  48. C
  49. C       DIF3   : ONE DIMENSIONAL VECTOR CONTAINING THE THIRD DERIVATIVE
  50. C                OF THE NODE POLYNOMIAL AT THE ZEROS
  51. C
  52. C       VECT   : ONE DIMENSIONAL VECTOR OF COMPUTED WEIGHTS
  53. C
  54. C     COMMON BLOCKS:      NONE
  55. C
  56. C     REQUIRED ROUTINES:  VILERR
  57. C
  58. C***********************************************************************
  59. C
  60.       INTEGER           J,NT,IER
  61.       DOUBLE PRECISION  AX,X,Y
  62.       DOUBLE PRECISION  ZERO,ONE,TWO,THREE
  63.       LOGICAL          LSTOP
  64. C
  65.       PARAMETER ( ZERO = 0.0D+00, ONE    = 1.0D+00,
  66.      +            TWO  = 2.0D+00, THREE  = 3.0D+00 )
  67. C
  68. C -- ERROR CHECKING
  69. C
  70.       IF ((N0 .NE. 0) .AND. (N0 .NE. 1)) THEN
  71.         IER   = 1
  72.         LSTOP = .TRUE.
  73.         CALL VILERR(IER,LSTOP)
  74.       ELSE
  75.       END IF
  76. C
  77.       IF ((N1 .NE. 0) .AND. (N1 .NE. 1)) THEN
  78.         IER   = 2
  79.         LSTOP = .TRUE.
  80.         CALL VILERR(IER,LSTOP)
  81.       ELSE
  82.       END IF
  83. C
  84.       IF (ND .LT. (N + N0 + N1)) THEN
  85.         IER   = 3
  86.         LSTOP = .TRUE.
  87.         CALL VILERR(IER,LSTOP)
  88.       ELSE
  89.       END IF
  90. C
  91.       IF ((ID .NE. 1) .AND. (ID.NE. 2) .AND. (ID .NE. 3)) THEN
  92.         IER   = 6
  93.         LSTOP = .TRUE.
  94.         CALL VILERR(IER,LSTOP)
  95.       ELSE
  96.       END IF
  97. C
  98.       IF (ID .NE. 3) THEN
  99.         IF (I .LT. 1) THEN
  100.           IER   = 4
  101.           LSTOP = .TRUE.
  102.           CALL VILERR(IER,LSTOP)
  103.         ELSE
  104.         END IF
  105. C
  106.         IF (I .GT. (N + N0 + N1)) THEN
  107.           IER   = 5
  108.           LSTOP = .TRUE.
  109.           CALL VILERR(IER,LSTOP)
  110.         ELSE
  111.         END IF
  112.       ELSE
  113.       END IF
  114. C
  115.       IF ((N + N0 + N1) .LT. 1) THEN
  116.         IER   = 7
  117.         LSTOP = .TRUE.
  118.         CALL VILERR(IER,LSTOP)
  119.       ELSE
  120.       END IF
  121. C
  122. C -- EVALUATE DISCRETIZATION MATRICES AND GAUSSIAN QUADRATURE
  123. C -- WEIGHTS.  QUADRATURE WEIGHTS ARE NORMALIZED TO SUM TO ONE.
  124. C
  125.       NT = N + N0 + N1
  126. C
  127.       IF (ID .NE. 3) THEN
  128.         DO 20 J = 1,NT
  129. C
  130.           IF (J .EQ. I) THEN
  131.             IF (ID .EQ. 1) THEN
  132.               VECT(I) = DIF2(I)/DIF1(I)/TWO
  133.             ELSE
  134.               VECT(I) = DIF3(I)/DIF1(I)/THREE
  135.             END IF
  136.           ELSE
  137.             Y       = ROOT(I) - ROOT(J)
  138.             VECT(J) = DIF1(I)/DIF1(J)/Y
  139.             IF (ID .EQ. 2) THEN
  140.               VECT(J) = VECT(J)*(DIF2(I)/DIF1(I) - TWO/Y)
  141.             ELSE
  142.             END IF
  143.           END IF
  144. C
  145.    20   CONTINUE
  146.       ELSE
  147.         Y = ZERO
  148. C
  149.         DO 25 J = 1,NT
  150. C
  151.           X  = ROOT(J)
  152.           AX = X*(ONE - X)
  153. C
  154.           IF(N0 .EQ. 0) THEN
  155.             AX = AX/X/X
  156.           ELSE
  157.           END IF
  158. C
  159.           IF(N1 .EQ. 0) THEN
  160.             AX = AX/(ONE - X)/(ONE - X)
  161.           ELSE
  162.           END IF
  163. C
  164.           VECT(J) = AX/DIF1(J)**2
  165.           Y       = Y + VECT(J)
  166. C
  167.    25   CONTINUE
  168. C
  169.         DO 60 J = 1,NT
  170.           VECT(J) = VECT(J)/Y
  171.    60   CONTINUE
  172. C
  173.       END IF
  174. C
  175.       RETURN
  176.       END
  177.