1*DECK PPADD 2 SUBROUTINE PPADD (N, IERROR, A, C, CBP, BP, BH) 3C***BEGIN PROLOGUE PPADD 4C***SUBSIDIARY 5C***PURPOSE Subsidiary to BLKTRI 6C***LIBRARY SLATEC 7C***TYPE SINGLE PRECISION (PPADD-S) 8C***AUTHOR (UNKNOWN) 9C***DESCRIPTION 10C 11C PPADD computes the eigenvalues of the periodic tridiagonal matrix 12C with coefficients AN,BN,CN. 13C 14C N is the order of the BH and BP polynomials. 15C BP contains the eigenvalues on output. 16C CBP is the same as BP except type complex. 17C BH is used to temporarily store the roots of the B HAT polynomial 18C which enters through BP. 19C 20C***SEE ALSO BLKTRI 21C***ROUTINES CALLED BSRH, PPSGF, PPSPF, PSGF 22C***COMMON BLOCKS CBLKT 23C***REVISION HISTORY (YYMMDD) 24C 801001 DATE WRITTEN 25C 890531 Changed all specific intrinsics to generic. (WRB) 26C 891214 Prologue converted to Version 4.0 format. (BAB) 27C 900402 Added TYPE section. (WRB) 28C***END PROLOGUE PPADD 29C 30 COMPLEX CX ,FSG ,HSG , 31 1 DD ,F ,FP ,FPP , 32 2 CDIS ,R1 ,R2 ,R3 , 33 3 CBP 34 DIMENSION A(*) ,C(*) ,BP(*) ,BH(*) , 35 1 CBP(*) 36 COMMON /CBLKT/ NPP ,K ,EPS ,CNV , 37 1 NM ,NCMPLX ,IK 38 EXTERNAL PSGF ,PPSPF ,PPSGF 39C***FIRST EXECUTABLE STATEMENT PPADD 40 SCNV = SQRT(CNV) 41 IZ = N 42 IF (BP(N)-BP(1)) 101,142,103 43 101 DO 102 J=1,N 44 NT = N-J 45 BH(J) = BP(NT+1) 46 102 CONTINUE 47 GO TO 105 48 103 DO 104 J=1,N 49 BH(J) = BP(J) 50 104 CONTINUE 51 105 NCMPLX = 0 52 MODIZ = MOD(IZ,2) 53 IS = 1 54 IF (MODIZ) 106,107,106 55 106 IF (A(1)) 110,142,107 56 107 XL = BH(1) 57 DB = BH(3)-BH(1) 58 108 XL = XL-DB 59 IF (PSGF(XL,IZ,C,A,BH)) 108,108,109 60 109 SGN = -1. 61 CBP(1) = CMPLX(BSRH(XL,BH(1),IZ,C,A,BH,PSGF,SGN),0.) 62 IS = 2 63 110 IF = IZ-1 64 IF (MODIZ) 111,112,111 65 111 IF (A(1)) 112,142,115 66 112 XR = BH(IZ) 67 DB = BH(IZ)-BH(IZ-2) 68 113 XR = XR+DB 69 IF (PSGF(XR,IZ,C,A,BH)) 113,114,114 70 114 SGN = 1. 71 CBP(IZ) = CMPLX(BSRH(BH(IZ),XR,IZ,C,A,BH,PSGF,SGN),0.) 72 IF = IZ-2 73 115 DO 136 IG=IS,IF,2 74 XL = BH(IG) 75 XR = BH(IG+1) 76 SGN = -1. 77 XM = BSRH(XL,XR,IZ,C,A,BH,PPSPF,SGN) 78 PSG = PSGF(XM,IZ,C,A,BH) 79 IF (ABS(PSG)-EPS) 118,118,116 80 116 IF (PSG*PPSGF(XM,IZ,C,A,BH)) 117,118,119 81C 82C CASE OF A REAL ZERO 83C 84 117 SGN = 1. 85 CBP(IG) = CMPLX(BSRH(BH(IG),XM,IZ,C,A,BH,PSGF,SGN),0.) 86 SGN = -1. 87 CBP(IG+1) = CMPLX(BSRH(XM,BH(IG+1),IZ,C,A,BH,PSGF,SGN),0.) 88 GO TO 136 89C 90C CASE OF A MULTIPLE ZERO 91C 92 118 CBP(IG) = CMPLX(XM,0.) 93 CBP(IG+1) = CMPLX(XM,0.) 94 GO TO 136 95C 96C CASE OF A COMPLEX ZERO 97C 98 119 IT = 0 99 ICV = 0 100 CX = CMPLX(XM,0.) 101 120 FSG = (1.,0.) 102 HSG = (1.,0.) 103 FP = (0.,0.) 104 FPP = (0.,0.) 105 DO 121 J=1,IZ 106 DD = 1./(CX-BH(J)) 107 FSG = FSG*A(J)*DD 108 HSG = HSG*C(J)*DD 109 FP = FP+DD 110 FPP = FPP-DD*DD 111 121 CONTINUE 112 IF (MODIZ) 123,122,123 113 122 F = (1.,0.)-FSG-HSG 114 GO TO 124 115 123 F = (1.,0.)+FSG+HSG 116 124 I3 = 0 117 IF (ABS(FP)) 126,126,125 118 125 I3 = 1 119 R3 = -F/FP 120 126 IF (ABS(FPP)) 132,132,127 121 127 CDIS = SQRT(FP**2-2.*F*FPP) 122 R1 = CDIS-FP 123 R2 = -FP-CDIS 124 IF (ABS(R1)-ABS(R2)) 129,129,128 125 128 R1 = R1/FPP 126 GO TO 130 127 129 R1 = R2/FPP 128 130 R2 = 2.*F/FPP/R1 129 IF (ABS(R2) .LT. ABS(R1)) R1 = R2 130 IF (I3) 133,133,131 131 131 IF (ABS(R3) .LT. ABS(R1)) R1 = R3 132 GO TO 133 133 132 R1 = R3 134 133 CX = CX+R1 135 IT = IT+1 136 IF (IT .GT. 50) GO TO 142 137 IF (ABS(R1) .GT. SCNV) GO TO 120 138 IF (ICV) 134,134,135 139 134 ICV = 1 140 GO TO 120 141 135 CBP(IG) = CX 142 CBP(IG+1) = CONJG(CX) 143 136 CONTINUE 144 IF (ABS(CBP(N))-ABS(CBP(1))) 137,142,139 145 137 NHALF = N/2 146 DO 138 J=1,NHALF 147 NT = N-J 148 CX = CBP(J) 149 CBP(J) = CBP(NT+1) 150 CBP(NT+1) = CX 151 138 CONTINUE 152 139 NCMPLX = 1 153 DO 140 J=2,IZ 154 IF (AIMAG(CBP(J))) 143,140,143 155 140 CONTINUE 156 NCMPLX = 0 157 DO 141 J=2,IZ 158 BP(J) = REAL(CBP(J)) 159 141 CONTINUE 160 GO TO 143 161 142 IERROR = 4 162 143 CONTINUE 163 RETURN 164 END 165