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