cmptr3.f 3.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113
  1. *DECK CMPTR3
  2. SUBROUTINE CMPTR3 (M, A, B, C, K, Y1, Y2, Y3, TCOS, D, W1, W2, W3)
  3. C***BEGIN PROLOGUE CMPTR3
  4. C***SUBSIDIARY
  5. C***PURPOSE Subsidiary to CMGNBN
  6. C***LIBRARY SLATEC
  7. C***TYPE COMPLEX (TRI3-S, CMPTR3-C)
  8. C***AUTHOR (UNKNOWN)
  9. C***DESCRIPTION
  10. C
  11. C Subroutine to solve tridiagonal systems.
  12. C
  13. C***SEE ALSO CMGNBN
  14. C***ROUTINES CALLED (NONE)
  15. C***REVISION HISTORY (YYMMDD)
  16. C 801001 DATE WRITTEN
  17. C 890206 REVISION DATE from Version 3.2
  18. C 891214 Prologue converted to Version 4.0 format. (BAB)
  19. C 900402 Added TYPE section. (WRB)
  20. C***END PROLOGUE CMPTR3
  21. COMPLEX A ,B ,C ,Y1 ,
  22. 1 Y2 ,Y3 ,TCOS ,D ,
  23. 2 W1 ,W2 ,W3 ,X ,
  24. 3 XX ,Z
  25. DIMENSION A(*) ,B(*) ,C(*) ,K(4) ,
  26. 1 TCOS(*) ,Y1(*) ,Y2(*) ,Y3(*) ,
  27. 2 D(*) ,W1(*) ,W2(*) ,W3(*)
  28. INTEGER K1P1, K2P1, K3P1, K4P1
  29. C
  30. C***FIRST EXECUTABLE STATEMENT CMPTR3
  31. MM1 = M-1
  32. K1 = K(1)
  33. K2 = K(2)
  34. K3 = K(3)
  35. K4 = K(4)
  36. K1P1 = K1+1
  37. K2P1 = K2+1
  38. K3P1 = K3+1
  39. K4P1 = K4+1
  40. K2K3K4 = K2+K3+K4
  41. IF (K2K3K4 .EQ. 0) GO TO 101
  42. L1 = K1P1/K2P1
  43. L2 = K1P1/K3P1
  44. L3 = K1P1/K4P1
  45. LINT1 = 1
  46. LINT2 = 1
  47. LINT3 = 1
  48. KINT1 = K1
  49. KINT2 = KINT1+K2
  50. KINT3 = KINT2+K3
  51. 101 CONTINUE
  52. DO 115 N=1,K1
  53. X = TCOS(N)
  54. IF (K2K3K4 .EQ. 0) GO TO 107
  55. IF (N .NE. L1) GO TO 103
  56. DO 102 I=1,M
  57. W1(I) = Y1(I)
  58. 102 CONTINUE
  59. 103 IF (N .NE. L2) GO TO 105
  60. DO 104 I=1,M
  61. W2(I) = Y2(I)
  62. 104 CONTINUE
  63. 105 IF (N .NE. L3) GO TO 107
  64. DO 106 I=1,M
  65. W3(I) = Y3(I)
  66. 106 CONTINUE
  67. 107 CONTINUE
  68. Z = 1./(B(1)-X)
  69. D(1) = C(1)*Z
  70. Y1(1) = Y1(1)*Z
  71. Y2(1) = Y2(1)*Z
  72. Y3(1) = Y3(1)*Z
  73. DO 108 I=2,M
  74. Z = 1./(B(I)-X-A(I)*D(I-1))
  75. D(I) = C(I)*Z
  76. Y1(I) = (Y1(I)-A(I)*Y1(I-1))*Z
  77. Y2(I) = (Y2(I)-A(I)*Y2(I-1))*Z
  78. Y3(I) = (Y3(I)-A(I)*Y3(I-1))*Z
  79. 108 CONTINUE
  80. DO 109 IP=1,MM1
  81. I = M-IP
  82. Y1(I) = Y1(I)-D(I)*Y1(I+1)
  83. Y2(I) = Y2(I)-D(I)*Y2(I+1)
  84. Y3(I) = Y3(I)-D(I)*Y3(I+1)
  85. 109 CONTINUE
  86. IF (K2K3K4 .EQ. 0) GO TO 115
  87. IF (N .NE. L1) GO TO 111
  88. I = LINT1+KINT1
  89. XX = X-TCOS(I)
  90. DO 110 I=1,M
  91. Y1(I) = XX*Y1(I)+W1(I)
  92. 110 CONTINUE
  93. LINT1 = LINT1+1
  94. L1 = (LINT1*K1P1)/K2P1
  95. 111 IF (N .NE. L2) GO TO 113
  96. I = LINT2+KINT2
  97. XX = X-TCOS(I)
  98. DO 112 I=1,M
  99. Y2(I) = XX*Y2(I)+W2(I)
  100. 112 CONTINUE
  101. LINT2 = LINT2+1
  102. L2 = (LINT2*K1P1)/K3P1
  103. 113 IF (N .NE. L3) GO TO 115
  104. I = LINT3+KINT3
  105. XX = X-TCOS(I)
  106. DO 114 I=1,M
  107. Y3(I) = XX*Y3(I)+W3(I)
  108. 114 CONTINUE
  109. LINT3 = LINT3+1
  110. L3 = (LINT3*K1P1)/K4P1
  111. 115 CONTINUE
  112. RETURN
  113. END