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