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