158
159
160
161
162
163
164 CHARACTER UPLO
165 INTEGER INFO, N
166
167
168 INTEGER IPIV( * )
169 COMPLEX*16 AP( * )
170
171
172
173
174
175 DOUBLE PRECISION ZERO, ONE
176 parameter( zero = 0.0d+0, one = 1.0d+0 )
177 DOUBLE PRECISION EIGHT, SEVTEN
178 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
179 COMPLEX*16 CONE
180 parameter( cone = ( 1.0d+0, 0.0d+0 ) )
181
182
183 LOGICAL UPPER
184 INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
185 $ KSTEP, KX, NPP
186 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX
187 COMPLEX*16 D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, ZDUM
188
189
190 LOGICAL LSAME
191 INTEGER IZAMAX
193
194
196
197
198 INTRINSIC abs, dble, dimag, max, sqrt
199
200
201 DOUBLE PRECISION CABS1
202
203
204 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
205
206
207
208
209
210 info = 0
211 upper =
lsame( uplo,
'U' )
212 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
213 info = -1
214 ELSE IF( n.LT.0 ) THEN
215 info = -2
216 END IF
217 IF( info.NE.0 ) THEN
218 CALL xerbla(
'ZSPTRF', -info )
219 RETURN
220 END IF
221
222
223
224 alpha = ( one+sqrt( sevten ) ) / eight
225
226 IF( upper ) THEN
227
228
229
230
231
232
233 k = n
234 kc = ( n-1 )*n / 2 + 1
235 10 CONTINUE
236 knc = kc
237
238
239
240 IF( k.LT.1 )
241 $ GO TO 110
242 kstep = 1
243
244
245
246
247 absakk = cabs1( ap( kc+k-1 ) )
248
249
250
251
252 IF( k.GT.1 ) THEN
253 imax =
izamax( k-1, ap( kc ), 1 )
254 colmax = cabs1( ap( kc+imax-1 ) )
255 ELSE
256 colmax = zero
257 END IF
258
259 IF( max( absakk, colmax ).EQ.zero ) THEN
260
261
262
263 IF( info.EQ.0 )
264 $ info = k
265 kp = k
266 ELSE
267 IF( absakk.GE.alpha*colmax ) THEN
268
269
270
271 kp = k
272 ELSE
273
274 rowmax = zero
275 jmax = imax
276 kx = imax*( imax+1 ) / 2 + imax
277 DO 20 j = imax + 1, k
278 IF( cabs1( ap( kx ) ).GT.rowmax ) THEN
279 rowmax = cabs1( ap( kx ) )
280 jmax = j
281 END IF
282 kx = kx + j
283 20 CONTINUE
284 kpc = ( imax-1 )*imax / 2 + 1
285 IF( imax.GT.1 ) THEN
286 jmax =
izamax( imax-1, ap( kpc ), 1 )
287 rowmax = max( rowmax, cabs1( ap( kpc+jmax-1 ) ) )
288 END IF
289
290 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
291
292
293
294 kp = k
295 ELSE IF( cabs1( ap( kpc+imax-1 ) ).GE.alpha*rowmax ) THEN
296
297
298
299
300 kp = imax
301 ELSE
302
303
304
305
306 kp = imax
307 kstep = 2
308 END IF
309 END IF
310
311 kk = k - kstep + 1
312 IF( kstep.EQ.2 )
313 $ knc = knc - k + 1
314 IF( kp.NE.kk ) THEN
315
316
317
318
319 CALL zswap( kp-1, ap( knc ), 1, ap( kpc ), 1 )
320 kx = kpc + kp - 1
321 DO 30 j = kp + 1, kk - 1
322 kx = kx + j - 1
323 t = ap( knc+j-1 )
324 ap( knc+j-1 ) = ap( kx )
325 ap( kx ) = t
326 30 CONTINUE
327 t = ap( knc+kk-1 )
328 ap( knc+kk-1 ) = ap( kpc+kp-1 )
329 ap( kpc+kp-1 ) = t
330 IF( kstep.EQ.2 ) THEN
331 t = ap( kc+k-2 )
332 ap( kc+k-2 ) = ap( kc+kp-1 )
333 ap( kc+kp-1 ) = t
334 END IF
335 END IF
336
337
338
339 IF( kstep.EQ.1 ) THEN
340
341
342
343
344
345
346
347
348
349
350
351 r1 = cone / ap( kc+k-1 )
352 CALL zspr( uplo, k-1, -r1, ap( kc ), 1, ap )
353
354
355
356 CALL zscal( k-1, r1, ap( kc ), 1 )
357 ELSE
358
359
360
361
362
363
364
365
366
367
368
369
370
371 IF( k.GT.2 ) THEN
372
373 d12 = ap( k-1+( k-1 )*k / 2 )
374 d22 = ap( k-1+( k-2 )*( k-1 ) / 2 ) / d12
375 d11 = ap( k+( k-1 )*k / 2 ) / d12
376 t = cone / ( d11*d22-cone )
377 d12 = t / d12
378
379 DO 50 j = k - 2, 1, -1
380 wkm1 = d12*( d11*ap( j+( k-2 )*( k-1 ) / 2 )-
381 $ ap( j+( k-1 )*k / 2 ) )
382 wk = d12*( d22*ap( j+( k-1 )*k / 2 )-
383 $ ap( j+( k-2 )*( k-1 ) / 2 ) )
384 DO 40 i = j, 1, -1
385 ap( i+( j-1 )*j / 2 ) = ap( i+( j-1 )*j / 2 ) -
386 $ ap( i+( k-1 )*k / 2 )*wk -
387 $ ap( i+( k-2 )*( k-1 ) / 2 )*wkm1
388 40 CONTINUE
389 ap( j+( k-1 )*k / 2 ) = wk
390 ap( j+( k-2 )*( k-1 ) / 2 ) = wkm1
391 50 CONTINUE
392
393 END IF
394 END IF
395 END IF
396
397
398
399 IF( kstep.EQ.1 ) THEN
400 ipiv( k ) = kp
401 ELSE
402 ipiv( k ) = -kp
403 ipiv( k-1 ) = -kp
404 END IF
405
406
407
408 k = k - kstep
409 kc = knc - k
410 GO TO 10
411
412 ELSE
413
414
415
416
417
418
419 k = 1
420 kc = 1
421 npp = n*( n+1 ) / 2
422 60 CONTINUE
423 knc = kc
424
425
426
427 IF( k.GT.n )
428 $ GO TO 110
429 kstep = 1
430
431
432
433
434 absakk = cabs1( ap( kc ) )
435
436
437
438
439 IF( k.LT.n ) THEN
440 imax = k +
izamax( n-k, ap( kc+1 ), 1 )
441 colmax = cabs1( ap( kc+imax-k ) )
442 ELSE
443 colmax = zero
444 END IF
445
446 IF( max( absakk, colmax ).EQ.zero ) THEN
447
448
449
450 IF( info.EQ.0 )
451 $ info = k
452 kp = k
453 ELSE
454 IF( absakk.GE.alpha*colmax ) THEN
455
456
457
458 kp = k
459 ELSE
460
461
462
463
464 rowmax = zero
465 kx = kc + imax - k
466 DO 70 j = k, imax - 1
467 IF( cabs1( ap( kx ) ).GT.rowmax ) THEN
468 rowmax = cabs1( ap( kx ) )
469 jmax = j
470 END IF
471 kx = kx + n - j
472 70 CONTINUE
473 kpc = npp - ( n-imax+1 )*( n-imax+2 ) / 2 + 1
474 IF( imax.LT.n ) THEN
475 jmax = imax +
izamax( n-imax, ap( kpc+1 ), 1 )
476 rowmax = max( rowmax, cabs1( ap( kpc+jmax-imax ) ) )
477 END IF
478
479 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
480
481
482
483 kp = k
484 ELSE IF( cabs1( ap( kpc ) ).GE.alpha*rowmax ) THEN
485
486
487
488
489 kp = imax
490 ELSE
491
492
493
494
495 kp = imax
496 kstep = 2
497 END IF
498 END IF
499
500 kk = k + kstep - 1
501 IF( kstep.EQ.2 )
502 $ knc = knc + n - k + 1
503 IF( kp.NE.kk ) THEN
504
505
506
507
508 IF( kp.LT.n )
509 $
CALL zswap( n-kp, ap( knc+kp-kk+1 ), 1, ap( kpc+1 ),
510 $ 1 )
511 kx = knc + kp - kk
512 DO 80 j = kk + 1, kp - 1
513 kx = kx + n - j + 1
514 t = ap( knc+j-kk )
515 ap( knc+j-kk ) = ap( kx )
516 ap( kx ) = t
517 80 CONTINUE
518 t = ap( knc )
519 ap( knc ) = ap( kpc )
520 ap( kpc ) = t
521 IF( kstep.EQ.2 ) THEN
522 t = ap( kc+1 )
523 ap( kc+1 ) = ap( kc+kp-k )
524 ap( kc+kp-k ) = t
525 END IF
526 END IF
527
528
529
530 IF( kstep.EQ.1 ) THEN
531
532
533
534
535
536
537
538 IF( k.LT.n ) THEN
539
540
541
542
543
544 r1 = cone / ap( kc )
545 CALL zspr( uplo, n-k, -r1, ap( kc+1 ), 1,
546 $ ap( kc+n-k+1 ) )
547
548
549
550 CALL zscal( n-k, r1, ap( kc+1 ), 1 )
551 END IF
552 ELSE
553
554
555
556
557
558
559
560
561 IF( k.LT.n-1 ) THEN
562
563
564
565
566
567
568
569
570
571 d21 = ap( k+1+( k-1 )*( 2*n-k ) / 2 )
572 d11 = ap( k+1+k*( 2*n-k-1 ) / 2 ) / d21
573 d22 = ap( k+( k-1 )*( 2*n-k ) / 2 ) / d21
574 t = cone / ( d11*d22-cone )
575 d21 = t / d21
576
577 DO 100 j = k + 2, n
578 wk = d21*( d11*ap( j+( k-1 )*( 2*n-k ) / 2 )-
579 $ ap( j+k*( 2*n-k-1 ) / 2 ) )
580 wkp1 = d21*( d22*ap( j+k*( 2*n-k-1 ) / 2 )-
581 $ ap( j+( k-1 )*( 2*n-k ) / 2 ) )
582 DO 90 i = j, n
583 ap( i+( j-1 )*( 2*n-j ) / 2 ) = ap( i+( j-1 )*
584 $ ( 2*n-j ) / 2 ) - ap( i+( k-1 )*( 2*n-k ) /
585 $ 2 )*wk - ap( i+k*( 2*n-k-1 ) / 2 )*wkp1
586 90 CONTINUE
587 ap( j+( k-1 )*( 2*n-k ) / 2 ) = wk
588 ap( j+k*( 2*n-k-1 ) / 2 ) = wkp1
589 100 CONTINUE
590 END IF
591 END IF
592 END IF
593
594
595
596 IF( kstep.EQ.1 ) THEN
597 ipiv( k ) = kp
598 ELSE
599 ipiv( k ) = -kp
600 ipiv( k+1 ) = -kp
601 END IF
602
603
604
605 k = k + kstep
606 kc = knc + n - k + 2
607 GO TO 60
608
609 END IF
610
611 110 CONTINUE
612 RETURN
613
614
615
subroutine xerbla(srname, info)
subroutine zspr(uplo, n, alpha, x, incx, ap)
ZSPR performs the symmetrical rank-1 update of a complex symmetric packed matrix.
integer function izamax(n, zx, incx)
IZAMAX
logical function lsame(ca, cb)
LSAME
subroutine zscal(n, za, zx, incx)
ZSCAL
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP