1 subroutine zhpsl(ap,n,kpvt,b) 2 integer n,kpvt(1) 3 complex*16 ap(1),b(1) 4c 5c zhisl solves the complex*16 hermitian system 6c a * x = b 7c using the factors computed by zhpfa. 8c 9c on entry 10c 11c ap complex*16(n*(n+1)/2) 12c the output from zhpfa. 13c 14c n integer 15c the order of the matrix a . 16c 17c kpvt integer(n) 18c the pivot vector from zhpfa. 19c 20c b complex*16(n) 21c the right hand side vector. 22c 23c on return 24c 25c b the solution vector x . 26c 27c error condition 28c 29c a division by zero may occur if zhpco has set rcond .eq. 0.0 30c or zhpfa has set info .ne. 0 . 31c 32c to compute inverse(a) * c where c is a matrix 33c with p columns 34c call zhpfa(ap,n,kpvt,info) 35c if (info .ne. 0) go to ... 36c do 10 j = 1, p 37c call zhpsl(ap,n,kpvt,c(1,j)) 38c 10 continue 39c 40c linpack. this version dated 08/14/78 . 41c james bunch, univ. calif. san diego, argonne nat. lab. 42c 43c subroutines and functions 44c 45c blas zaxpy,zdotc 46c fortran dconjg,iabs 47c 48c internal variables. 49c 50 complex*16 ak,akm1,bk,bkm1,zdotc,denom,temp 51 integer ik,ikm1,ikp1,k,kk,km1k,km1km1,kp 52c 53c loop backward applying the transformations and 54c d inverse to b. 55c 56 k = n 57 ik = (n*(n - 1))/2 58 10 if (k .eq. 0) go to 80 59 kk = ik + k 60 if (kpvt(k) .lt. 0) go to 40 61c 62c 1 x 1 pivot block. 63c 64 if (k .eq. 1) go to 30 65 kp = kpvt(k) 66 if (kp .eq. k) go to 20 67c 68c interchange. 69c 70 temp = b(k) 71 b(k) = b(kp) 72 b(kp) = temp 73 20 continue 74c 75c apply the transformation. 76c 77 call zaxpy(k-1,b(k),ap(ik+1),1,b(1),1) 78 30 continue 79c 80c apply d inverse. 81c 82 b(k) = b(k)/ap(kk) 83 k = k - 1 84 ik = ik - k 85 go to 70 86 40 continue 87c 88c 2 x 2 pivot block. 89c 90 ikm1 = ik - (k - 1) 91 if (k .eq. 2) go to 60 92 kp = iabs(kpvt(k)) 93 if (kp .eq. k - 1) go to 50 94c 95c interchange. 96c 97 temp = b(k-1) 98 b(k-1) = b(kp) 99 b(kp) = temp 100 50 continue 101c 102c apply the transformation. 103c 104 call zaxpy(k-2,b(k),ap(ik+1),1,b(1),1) 105 call zaxpy(k-2,b(k-1),ap(ikm1+1),1,b(1),1) 106 60 continue 107c 108c apply d inverse. 109c 110 km1k = ik + k - 1 111 kk = ik + k 112 ak = ap(kk)/dconjg(ap(km1k)) 113 km1km1 = ikm1 + k - 1 114 akm1 = ap(km1km1)/ap(km1k) 115 bk = b(k)/dconjg(ap(km1k)) 116 bkm1 = b(k-1)/ap(km1k) 117 denom = ak*akm1 - 1.0d0 118 b(k) = (akm1*bk - bkm1)/denom 119 b(k-1) = (ak*bkm1 - bk)/denom 120 k = k - 2 121 ik = ik - (k + 1) - k 122 70 continue 123 go to 10 124 80 continue 125c 126c loop forward applying the transformations. 127c 128 k = 1 129 ik = 0 130 90 if (k .gt. n) go to 160 131 if (kpvt(k) .lt. 0) go to 120 132c 133c 1 x 1 pivot block. 134c 135 if (k .eq. 1) go to 110 136c 137c apply the transformation. 138c 139 b(k) = b(k) + zdotc(k-1,ap(ik+1),1,b(1),1) 140 kp = kpvt(k) 141 if (kp .eq. k) go to 100 142c 143c interchange. 144c 145 temp = b(k) 146 b(k) = b(kp) 147 b(kp) = temp 148 100 continue 149 110 continue 150 ik = ik + k 151 k = k + 1 152 go to 150 153 120 continue 154c 155c 2 x 2 pivot block. 156c 157 if (k .eq. 1) go to 140 158c 159c apply the transformation. 160c 161 b(k) = b(k) + zdotc(k-1,ap(ik+1),1,b(1),1) 162 ikp1 = ik + k 163 b(k+1) = b(k+1) + zdotc(k-1,ap(ikp1+1),1,b(1),1) 164 kp = iabs(kpvt(k)) 165 if (kp .eq. k) go to 130 166c 167c interchange. 168c 169 temp = b(k) 170 b(k) = b(kp) 171 b(kp) = temp 172 130 continue 173 140 continue 174 ik = ik + k + k + 1 175 k = k + 2 176 150 continue 177 go to 90 178 160 continue 179 return 180 end 181