127
128
129
130
131
132
133 CHARACTER UPLO
134 INTEGER INFO, LDA, LDB, N, NRHS
135
136
137 INTEGER IPIV( * )
138 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
139
140
141
142
143
144 COMPLEX*16 ONE
145 parameter( one = (1.0d+0,0.0d+0) )
146
147
148 LOGICAL UPPER
149 INTEGER I, IINFO, J, K, KP
150 DOUBLE PRECISION S
151 COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
152
153
154 LOGICAL LSAME
156
157
159
160
161 INTRINSIC dble, dconjg, max
162
163
164
165 info = 0
166 upper =
lsame( uplo,
'U' )
167 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
168 info = -1
169 ELSE IF( n.LT.0 ) THEN
170 info = -2
171 ELSE IF( nrhs.LT.0 ) THEN
172 info = -3
173 ELSE IF( lda.LT.max( 1, n ) ) THEN
174 info = -5
175 ELSE IF( ldb.LT.max( 1, n ) ) THEN
176 info = -8
177 END IF
178 IF( info.NE.0 ) THEN
179 CALL xerbla(
'ZHETRS2', -info )
180 RETURN
181 END IF
182
183
184
185 IF( n.EQ.0 .OR. nrhs.EQ.0 )
186 $ RETURN
187
188
189
190 CALL zsyconv( uplo,
'C', n, a, lda, ipiv, work, iinfo )
191
192 IF( upper ) THEN
193
194
195
196
197 k=n
198 DO WHILE ( k .GE. 1 )
199 IF( ipiv( k ).GT.0 ) THEN
200
201
202 kp = ipiv( k )
203 IF( kp.NE.k )
204 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
205 k=k-1
206 ELSE
207
208
209 kp = -ipiv( k )
210 IF( kp.EQ.-ipiv( k-1 ) )
211 $
CALL zswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
212 k=k-2
213 END IF
214 END DO
215
216
217
218 CALL ztrsm(
'L',
'U',
'N',
'U',n,nrhs,one,a,lda,b,ldb)
219
220
221
222 i=n
223 DO WHILE ( i .GE. 1 )
224 IF( ipiv(i) .GT. 0 ) THEN
225 s = dble( one ) / dble( a( i, i ) )
226 CALL zdscal( nrhs, s, b( i, 1 ), ldb )
227 ELSEIF ( i .GT. 1) THEN
228 IF ( ipiv(i-1) .EQ. ipiv(i) ) THEN
229 akm1k = work(i)
230 akm1 = a( i-1, i-1 ) / akm1k
231 ak = a( i, i ) / dconjg( akm1k )
232 denom = akm1*ak - one
233 DO 15 j = 1, nrhs
234 bkm1 = b( i-1, j ) / akm1k
235 bk = b( i, j ) / dconjg( akm1k )
236 b( i-1, j ) = ( ak*bkm1-bk ) / denom
237 b( i, j ) = ( akm1*bk-bkm1 ) / denom
238 15 CONTINUE
239 i = i - 1
240 ENDIF
241 ENDIF
242 i = i - 1
243 END DO
244
245
246
247 CALL ztrsm(
'L',
'U',
'C',
'U',n,nrhs,one,a,lda,b,ldb)
248
249
250
251 k=1
252 DO WHILE ( k .LE. n )
253 IF( ipiv( k ).GT.0 ) THEN
254
255
256 kp = ipiv( k )
257 IF( kp.NE.k )
258 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
259 k=k+1
260 ELSE
261
262
263 kp = -ipiv( k )
264 IF( k .LT. n .AND. kp.EQ.-ipiv( k+1 ) )
265 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
266 k=k+2
267 ENDIF
268 END DO
269
270 ELSE
271
272
273
274
275 k=1
276 DO WHILE ( k .LE. n )
277 IF( ipiv( k ).GT.0 ) THEN
278
279
280 kp = ipiv( k )
281 IF( kp.NE.k )
282 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
283 k=k+1
284 ELSE
285
286
287 kp = -ipiv( k+1 )
288 IF( kp.EQ.-ipiv( k ) )
289 $
CALL zswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
290 k=k+2
291 ENDIF
292 END DO
293
294
295
296 CALL ztrsm(
'L',
'L',
'N',
'U',n,nrhs,one,a,lda,b,ldb)
297
298
299
300 i=1
301 DO WHILE ( i .LE. n )
302 IF( ipiv(i) .GT. 0 ) THEN
303 s = dble( one ) / dble( a( i, i ) )
304 CALL zdscal( nrhs, s, b( i, 1 ), ldb )
305 ELSE
306 akm1k = work(i)
307 akm1 = a( i, i ) / dconjg( akm1k )
308 ak = a( i+1, i+1 ) / akm1k
309 denom = akm1*ak - one
310 DO 25 j = 1, nrhs
311 bkm1 = b( i, j ) / dconjg( akm1k )
312 bk = b( i+1, j ) / akm1k
313 b( i, j ) = ( ak*bkm1-bk ) / denom
314 b( i+1, j ) = ( akm1*bk-bkm1 ) / denom
315 25 CONTINUE
316 i = i + 1
317 ENDIF
318 i = i + 1
319 END DO
320
321
322
323 CALL ztrsm(
'L',
'L',
'C',
'U',n,nrhs,one,a,lda,b,ldb)
324
325
326
327 k=n
328 DO WHILE ( k .GE. 1 )
329 IF( ipiv( k ).GT.0 ) THEN
330
331
332 kp = ipiv( k )
333 IF( kp.NE.k )
334 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
335 k=k-1
336 ELSE
337
338
339 kp = -ipiv( k )
340 IF( k.GT.1 .AND. kp.EQ.-ipiv( k-1 ) )
341 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
342 k=k-2
343 ENDIF
344 END DO
345
346 END IF
347
348
349
350 CALL zsyconv( uplo,
'R', n, a, lda, ipiv, work, iinfo )
351
352 RETURN
353
354
355
subroutine xerbla(srname, info)
logical function lsame(ca, cb)
LSAME
subroutine zdscal(n, da, zx, incx)
ZDSCAL
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
subroutine zsyconv(uplo, way, n, a, lda, ipiv, e, info)
ZSYCONV
subroutine ztrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRSM