1      SUBROUTINE BLKTR1 (N, AN, BN, CN, M, AM, BM, CM, IDIMY, Y, B, W1,
2     +   W2, W3, WD, WW, WU, PRDCT, CPRDCT)
3C***BEGIN PROLOGUE  BLKTR1
4C***SUBSIDIARY
5C***PURPOSE  Subsidiary to BLKTRI
6C***LIBRARY   SLATEC
7C***TYPE      SINGLE PRECISION (BLKTR1-S, CBLKT1-C)
8C***AUTHOR  (UNKNOWN)
9C***DESCRIPTION
10C
11C BLKTR1 solves the linear system set up by BLKTRI.
12C
13C B  contains the roots of all the B polynomials.
14C W1,W2,W3,WD,WW,WU  are all working arrays.
15C PRDCT  is either PRODP or PROD depending on whether the boundary
16C conditions in the M direction are periodic or not.
17C CPRDCT is either CPRODP or CPROD which are the complex versions
18C of PRODP and PROD. These are called in the event that some
19C of the roots of the B sub P polynomial are complex.
20C
21C***SEE ALSO  BLKTRI
22C***ROUTINES CALLED  INDXA, INDXB, INDXC
23C***COMMON BLOCKS    CBLKT
24C***REVISION HISTORY  (YYMMDD)
25C   801001  DATE WRITTEN
26C   891214  Prologue converted to Version 4.0 format.  (BAB)
27C   900402  Added TYPE section.  (WRB)
28C***END PROLOGUE  BLKTR1
29C
30      DIMENSION       AN(*)      ,BN(*)      ,CN(*)      ,AM(*)      ,
31     1                BM(*)      ,CM(*)      ,B(*)       ,W1(*)      ,
32     2                W2(*)      ,W3(*)      ,WD(*)      ,WW(*)      ,
33     3                WU(*)      ,Y(IDIMY,*)
34      COMMON /CBLKT/  NPP        ,K          ,EPS        ,CNV        ,
35     1                NM         ,NCMPLX     ,IK
36C***FIRST EXECUTABLE STATEMENT  BLKTR1
37      KDO = K-1
38      DO 109 L=1,KDO
39         IR = L-1
40         I2 = 2**IR
41         I1 = I2/2
42         I3 = I2+I1
43         I4 = I2+I2
44         IRM1 = IR-1
45         CALL INDXB (I2,IR,IM2,NM2)
46         CALL INDXB (I1,IRM1,IM3,NM3)
47         CALL INDXB (I3,IRM1,IM1,NM1)
48         CALL PRDCT (NM2,B(IM2),NM3,B(IM3),NM1,B(IM1),0,DUM,Y(1,I2),W3,
49     1               M,AM,BM,CM,WD,WW,WU)
50         IF = 2**K
51         DO 108 I=I4,IF,I4
52            IF (I-NM) 101,101,108
53  101       IPI1 = I+I1
54            IPI2 = I+I2
55            IPI3 = I+I3
56            CALL INDXC (I,IR,IDXC,NC)
57            IF (I-IF) 102,108,108
58  102       CALL INDXA (I,IR,IDXA,NA)
59            CALL INDXB (I-I1,IRM1,IM1,NM1)
60            CALL INDXB (IPI2,IR,IP2,NP2)
61            CALL INDXB (IPI1,IRM1,IP1,NP1)
62            CALL INDXB (IPI3,IRM1,IP3,NP3)
63            CALL PRDCT (NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),W3,W1,M,AM,
64     1                  BM,CM,WD,WW,WU)
65            IF (IPI2-NM) 105,105,103
66  103       DO 104 J=1,M
67               W3(J) = 0.
68               W2(J) = 0.
69  104       CONTINUE
70            GO TO 106
71  105       CALL PRDCT (NP2,B(IP2),NP1,B(IP1),NP3,B(IP3),0,DUM,
72     1                  Y(1,IPI2),W3,M,AM,BM,CM,WD,WW,WU)
73            CALL PRDCT (NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),W3,W2,M,AM,
74     1                  BM,CM,WD,WW,WU)
75  106       DO 107 J=1,M
76               Y(J,I) = W1(J)+W2(J)+Y(J,I)
77  107       CONTINUE
78  108    CONTINUE
79  109 CONTINUE
80      IF (NPP) 132,110,132
81C
82C     THE PERIODIC CASE IS TREATED USING THE CAPACITANCE MATRIX METHOD
83C
84  110 IF = 2**K
85      I = IF/2
86      I1 = I/2
87      CALL INDXB (I-I1,K-2,IM1,NM1)
88      CALL INDXB (I+I1,K-2,IP1,NP1)
89      CALL INDXB (I,K-1,IZ,NZ)
90      CALL PRDCT (NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,Y(1,I),W1,M,AM,
91     1            BM,CM,WD,WW,WU)
92      IZR = I
93      DO 111 J=1,M
94         W2(J) = W1(J)
95  111 CONTINUE
96      DO 113 LL=2,K
97         L = K-LL+1
98         IR = L-1
99         I2 = 2**IR
100         I1 = I2/2
101         I = I2
102         CALL INDXC (I,IR,IDXC,NC)
103         CALL INDXB (I,IR,IZ,NZ)
104         CALL INDXB (I-I1,IR-1,IM1,NM1)
105         CALL INDXB (I+I1,IR-1,IP1,NP1)
106         CALL PRDCT (NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),W1,W1,M,AM,BM,
107     1               CM,WD,WW,WU)
108         DO 112 J=1,M
109            W1(J) = Y(J,I)+W1(J)
110  112    CONTINUE
111         CALL PRDCT (NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,W1,W1,M,AM,
112     1               BM,CM,WD,WW,WU)
113  113 CONTINUE
114      DO 118 LL=2,K
115         L = K-LL+1
116         IR = L-1
117         I2 = 2**IR
118         I1 = I2/2
119         I4 = I2+I2
120         IFD = IF-I2
121         DO 117 I=I2,IFD,I4
122            IF (I-I2-IZR) 117,114,117
123  114       IF (I-NM) 115,115,118
124  115       CALL INDXA (I,IR,IDXA,NA)
125            CALL INDXB (I,IR,IZ,NZ)
126            CALL INDXB (I-I1,IR-1,IM1,NM1)
127            CALL INDXB (I+I1,IR-1,IP1,NP1)
128            CALL PRDCT (NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),W2,W2,M,AM,
129     1                  BM,CM,WD,WW,WU)
130            DO 116 J=1,M
131               W2(J) = Y(J,I)+W2(J)
132  116       CONTINUE
133            CALL PRDCT (NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,W2,W2,M,
134     1                  AM,BM,CM,WD,WW,WU)
135            IZR = I
136            IF (I-NM) 117,119,117
137  117    CONTINUE
138  118 CONTINUE
139  119 DO 120 J=1,M
140         Y(J,NM+1) = Y(J,NM+1)-CN(NM+1)*W1(J)-AN(NM+1)*W2(J)
141  120 CONTINUE
142      CALL INDXB (IF/2,K-1,IM1,NM1)
143      CALL INDXB (IF,K-1,IP,NP)
144      IF (NCMPLX) 121,122,121
145  121 CALL CPRDCT (NM+1,B(IP),NM1,B(IM1),0,DUM,0,DUM,Y(1,NM+1),
146     1             Y(1,NM+1),M,AM,BM,CM,W1,W3,WW)
147      GO TO 123
148  122 CALL PRDCT (NM+1,B(IP),NM1,B(IM1),0,DUM,0,DUM,Y(1,NM+1),
149     1            Y(1,NM+1),M,AM,BM,CM,WD,WW,WU)
150  123 DO 124 J=1,M
151         W1(J) = AN(1)*Y(J,NM+1)
152         W2(J) = CN(NM)*Y(J,NM+1)
153         Y(J,1) = Y(J,1)-W1(J)
154         Y(J,NM) = Y(J,NM)-W2(J)
155  124 CONTINUE
156      DO 126 L=1,KDO
157         IR = L-1
158         I2 = 2**IR
159         I4 = I2+I2
160         I1 = I2/2
161         I = I4
162         CALL INDXA (I,IR,IDXA,NA)
163         CALL INDXB (I-I2,IR,IM2,NM2)
164         CALL INDXB (I-I2-I1,IR-1,IM3,NM3)
165         CALL INDXB (I-I1,IR-1,IM1,NM1)
166         CALL PRDCT (NM2,B(IM2),NM3,B(IM3),NM1,B(IM1),0,DUM,W1,W1,M,AM,
167     1               BM,CM,WD,WW,WU)
168         CALL PRDCT (NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),W1,W1,M,AM,BM,
169     1               CM,WD,WW,WU)
170         DO 125 J=1,M
171            Y(J,I) = Y(J,I)-W1(J)
172  125    CONTINUE
173  126 CONTINUE
174C
175      IZR = NM
176      DO 131 L=1,KDO
177         IR = L-1
178         I2 = 2**IR
179         I1 = I2/2
180         I3 = I2+I1
181         I4 = I2+I2
182         IRM1 = IR-1
183         DO 130 I=I4,IF,I4
184            IPI1 = I+I1
185            IPI2 = I+I2
186            IPI3 = I+I3
187            IF (IPI2-IZR) 127,128,127
188  127       IF (I-IZR) 130,131,130
189  128       CALL INDXC (I,IR,IDXC,NC)
190            CALL INDXB (IPI2,IR,IP2,NP2)
191            CALL INDXB (IPI1,IRM1,IP1,NP1)
192            CALL INDXB (IPI3,IRM1,IP3,NP3)
193            CALL PRDCT (NP2,B(IP2),NP1,B(IP1),NP3,B(IP3),0,DUM,W2,W2,M,
194     1                  AM,BM,CM,WD,WW,WU)
195            CALL PRDCT (NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),W2,W2,M,AM,
196     1                  BM,CM,WD,WW,WU)
197            DO 129 J=1,M
198               Y(J,I) = Y(J,I)-W2(J)
199  129       CONTINUE
200            IZR = I
201            GO TO 131
202  130    CONTINUE
203  131 CONTINUE
204C
205C BEGIN BACK SUBSTITUTION PHASE
206C
207  132 DO 144 LL=1,K
208         L = K-LL+1
209         IR = L-1
210         IRM1 = IR-1
211         I2 = 2**IR
212         I1 = I2/2
213         I4 = I2+I2
214         IFD = IF-I2
215         DO 143 I=I2,IFD,I4
216            IF (I-NM) 133,133,143
217  133       IMI1 = I-I1
218            IMI2 = I-I2
219            IPI1 = I+I1
220            IPI2 = I+I2
221            CALL INDXA (I,IR,IDXA,NA)
222            CALL INDXC (I,IR,IDXC,NC)
223            CALL INDXB (I,IR,IZ,NZ)
224            CALL INDXB (IMI1,IRM1,IM1,NM1)
225            CALL INDXB (IPI1,IRM1,IP1,NP1)
226            IF (I-I2) 134,134,136
227  134       DO 135 J=1,M
228               W1(J) = 0.
229  135       CONTINUE
230            GO TO 137
231  136       CALL PRDCT (NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),Y(1,IMI2),
232     1                  W1,M,AM,BM,CM,WD,WW,WU)
233  137       IF (IPI2-NM) 140,140,138
234  138       DO 139 J=1,M
235               W2(J) = 0.
236  139       CONTINUE
237            GO TO 141
238  140       CALL PRDCT (NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),Y(1,IPI2),
239     1                  W2,M,AM,BM,CM,WD,WW,WU)
240  141       DO 142 J=1,M
241               W1(J) = Y(J,I)+W1(J)+W2(J)
242  142       CONTINUE
243            CALL PRDCT (NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,W1,Y(1,I),
244     1                  M,AM,BM,CM,WD,WW,WU)
245  143    CONTINUE
246  144 CONTINUE
247      RETURN
248      END
249