113
114
115
116
117
118
119 CHARACTER UPLO
120 INTEGER INFO, LDB, N, NRHS
121
122
123 INTEGER IPIV( * )
124 COMPLEX*16 AP( * ), B( LDB, * )
125
126
127
128
129
130 COMPLEX*16 ONE
131 parameter( one = ( 1.0d+0, 0.0d+0 ) )
132
133
134 LOGICAL UPPER
135 INTEGER J, K, KC, KP
136 DOUBLE PRECISION S
137 COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
138
139
140 LOGICAL LSAME
142
143
146
147
148 INTRINSIC dble, dconjg, max
149
150
151
152 info = 0
153 upper =
lsame( uplo,
'U' )
154 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
155 info = -1
156 ELSE IF( n.LT.0 ) THEN
157 info = -2
158 ELSE IF( nrhs.LT.0 ) THEN
159 info = -3
160 ELSE IF( ldb.LT.max( 1, n ) ) THEN
161 info = -7
162 END IF
163 IF( info.NE.0 ) THEN
164 CALL xerbla(
'ZHPTRS', -info )
165 RETURN
166 END IF
167
168
169
170 IF( n.EQ.0 .OR. nrhs.EQ.0 )
171 $ RETURN
172
173 IF( upper ) THEN
174
175
176
177
178
179
180
181
182 k = n
183 kc = n*( n+1 ) / 2 + 1
184 10 CONTINUE
185
186
187
188 IF( k.LT.1 )
189 $ GO TO 30
190
191 kc = kc - k
192 IF( ipiv( k ).GT.0 ) THEN
193
194
195
196
197
198 kp = ipiv( k )
199 IF( kp.NE.k )
200 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
201
202
203
204
205 CALL zgeru( k-1, nrhs, -one, ap( kc ), 1, b( k, 1 ), ldb,
206 $ b( 1, 1 ), ldb )
207
208
209
210 s = dble( one ) / dble( ap( kc+k-1 ) )
211 CALL zdscal( nrhs, s, b( k, 1 ), ldb )
212 k = k - 1
213 ELSE
214
215
216
217
218
219 kp = -ipiv( k )
220 IF( kp.NE.k-1 )
221 $
CALL zswap( nrhs, b( k-1, 1 ), ldb, b( kp, 1 ), ldb )
222
223
224
225
226 CALL zgeru( k-2, nrhs, -one, ap( kc ), 1, b( k, 1 ), ldb,
227 $ b( 1, 1 ), ldb )
228 CALL zgeru( k-2, nrhs, -one, ap( kc-( k-1 ) ), 1,
229 $ b( k-1, 1 ), ldb, b( 1, 1 ), ldb )
230
231
232
233 akm1k = ap( kc+k-2 )
234 akm1 = ap( kc-1 ) / akm1k
235 ak = ap( kc+k-1 ) / dconjg( akm1k )
236 denom = akm1*ak - one
237 DO 20 j = 1, nrhs
238 bkm1 = b( k-1, j ) / akm1k
239 bk = b( k, j ) / dconjg( akm1k )
240 b( k-1, j ) = ( ak*bkm1-bk ) / denom
241 b( k, j ) = ( akm1*bk-bkm1 ) / denom
242 20 CONTINUE
243 kc = kc - k + 1
244 k = k - 2
245 END IF
246
247 GO TO 10
248 30 CONTINUE
249
250
251
252
253
254
255 k = 1
256 kc = 1
257 40 CONTINUE
258
259
260
261 IF( k.GT.n )
262 $ GO TO 50
263
264 IF( ipiv( k ).GT.0 ) THEN
265
266
267
268
269
270
271 IF( k.GT.1 ) THEN
272 CALL zlacgv( nrhs, b( k, 1 ), ldb )
273 CALL zgemv(
'Conjugate transpose', k-1, nrhs, -one, b,
274 $ ldb, ap( kc ), 1, one, b( k, 1 ), ldb )
275 CALL zlacgv( nrhs, b( k, 1 ), ldb )
276 END IF
277
278
279
280 kp = ipiv( k )
281 IF( kp.NE.k )
282 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
283 kc = kc + k
284 k = k + 1
285 ELSE
286
287
288
289
290
291
292 IF( k.GT.1 ) THEN
293 CALL zlacgv( nrhs, b( k, 1 ), ldb )
294 CALL zgemv(
'Conjugate transpose', k-1, nrhs, -one, b,
295 $ ldb, ap( kc ), 1, one, b( k, 1 ), ldb )
296 CALL zlacgv( nrhs, b( k, 1 ), ldb )
297
298 CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
299 CALL zgemv(
'Conjugate transpose', k-1, nrhs, -one, b,
300 $ ldb, ap( kc+k ), 1, one, b( k+1, 1 ), ldb )
301 CALL zlacgv( nrhs, b( k+1, 1 ), ldb )
302 END IF
303
304
305
306 kp = -ipiv( k )
307 IF( kp.NE.k )
308 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
309 kc = kc + 2*k + 1
310 k = k + 2
311 END IF
312
313 GO TO 40
314 50 CONTINUE
315
316 ELSE
317
318
319
320
321
322
323
324
325 k = 1
326 kc = 1
327 60 CONTINUE
328
329
330
331 IF( k.GT.n )
332 $ GO TO 80
333
334 IF( ipiv( k ).GT.0 ) THEN
335
336
337
338
339
340 kp = ipiv( k )
341 IF( kp.NE.k )
342 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
343
344
345
346
347 IF( k.LT.n )
348 $
CALL zgeru( n-k, nrhs, -one, ap( kc+1 ), 1, b( k, 1 ),
349 $ ldb, b( k+1, 1 ), ldb )
350
351
352
353 s = dble( one ) / dble( ap( kc ) )
354 CALL zdscal( nrhs, s, b( k, 1 ), ldb )
355 kc = kc + n - k + 1
356 k = k + 1
357 ELSE
358
359
360
361
362
363 kp = -ipiv( k )
364 IF( kp.NE.k+1 )
365 $
CALL zswap( nrhs, b( k+1, 1 ), ldb, b( kp, 1 ), ldb )
366
367
368
369
370 IF( k.LT.n-1 ) THEN
371 CALL zgeru( n-k-1, nrhs, -one, ap( kc+2 ), 1, b( k,
372 $ 1 ),
373 $ ldb, b( k+2, 1 ), ldb )
374 CALL zgeru( n-k-1, nrhs, -one, ap( kc+n-k+2 ), 1,
375 $ b( k+1, 1 ), ldb, b( k+2, 1 ), ldb )
376 END IF
377
378
379
380 akm1k = ap( kc+1 )
381 akm1 = ap( kc ) / dconjg( akm1k )
382 ak = ap( kc+n-k+1 ) / akm1k
383 denom = akm1*ak - one
384 DO 70 j = 1, nrhs
385 bkm1 = b( k, j ) / dconjg( akm1k )
386 bk = b( k+1, j ) / akm1k
387 b( k, j ) = ( ak*bkm1-bk ) / denom
388 b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
389 70 CONTINUE
390 kc = kc + 2*( n-k ) + 1
391 k = k + 2
392 END IF
393
394 GO TO 60
395 80 CONTINUE
396
397
398
399
400
401
402 k = n
403 kc = n*( n+1 ) / 2 + 1
404 90 CONTINUE
405
406
407
408 IF( k.LT.1 )
409 $ GO TO 100
410
411 kc = kc - ( n-k+1 )
412 IF( ipiv( k ).GT.0 ) THEN
413
414
415
416
417
418
419 IF( k.LT.n ) THEN
420 CALL zlacgv( nrhs, b( k, 1 ), ldb )
421 CALL zgemv(
'Conjugate transpose', n-k, nrhs, -one,
422 $ b( k+1, 1 ), ldb, ap( kc+1 ), 1, one,
423 $ b( k, 1 ), ldb )
424 CALL zlacgv( nrhs, b( k, 1 ), ldb )
425 END IF
426
427
428
429 kp = ipiv( k )
430 IF( kp.NE.k )
431 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
432 k = k - 1
433 ELSE
434
435
436
437
438
439
440 IF( k.LT.n ) THEN
441 CALL zlacgv( nrhs, b( k, 1 ), ldb )
442 CALL zgemv(
'Conjugate transpose', n-k, nrhs, -one,
443 $ b( k+1, 1 ), ldb, ap( kc+1 ), 1, one,
444 $ b( k, 1 ), ldb )
445 CALL zlacgv( nrhs, b( k, 1 ), ldb )
446
447 CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
448 CALL zgemv(
'Conjugate transpose', n-k, nrhs, -one,
449 $ b( k+1, 1 ), ldb, ap( kc-( n-k ) ), 1, one,
450 $ b( k-1, 1 ), ldb )
451 CALL zlacgv( nrhs, b( k-1, 1 ), ldb )
452 END IF
453
454
455
456 kp = -ipiv( k )
457 IF( kp.NE.k )
458 $
CALL zswap( nrhs, b( k, 1 ), ldb, b( kp, 1 ), ldb )
459 kc = kc - ( n-k+2 )
460 k = k - 2
461 END IF
462
463 GO TO 90
464 100 CONTINUE
465 END IF
466
467 RETURN
468
469
470
subroutine xerbla(srname, info)
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
subroutine zgeru(m, n, alpha, x, incx, y, incy, a, lda)
ZGERU
subroutine zlacgv(n, x, incx)
ZLACGV conjugates a complex vector.
logical function lsame(ca, cb)
LSAME
subroutine zdscal(n, da, zx, incx)
ZDSCAL
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP