LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ ssptrf()

subroutine ssptrf ( character uplo,
integer n,
real, dimension( * ) ap,
integer, dimension( * ) ipiv,
integer info )

SSPTRF

Download SSPTRF + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SSPTRF computes the factorization of a real symmetric matrix A stored
!> in packed format using the Bunch-Kaufman diagonal pivoting method:
!>
!>    A = U*D*U**T  or  A = L*D*L**T
!>
!> where U (or L) is a product of permutation and unit upper (lower)
!> triangular matrices, and D is symmetric and block diagonal with
!> 1-by-1 and 2-by-2 diagonal blocks.
!> 
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!> 
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!> 
[in,out]AP
!>          AP is REAL array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the symmetric matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
!>
!>          On exit, the block diagonal matrix D and the multipliers used
!>          to obtain the factor U or L, stored as a packed triangular
!>          matrix overwriting A (see below for further details).
!> 
[out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          Details of the interchanges and the block structure of D.
!>          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
!>          interchanged and D(k,k) is a 1-by-1 diagonal block.
!>          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
!>          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
!>          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
!>          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
!>          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
!> 
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>          > 0: if INFO = i, D(i,i) is exactly zero.  The factorization
!>               has been completed, but the block diagonal matrix D is
!>               exactly singular, and division by zero will occur if it
!>               is used to solve a system of equations.
!> 
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  5-96 - Based on modifications by J. Lewis, Boeing Computer Services
!>         Company
!>
!>  If UPLO = 'U', then A = U*D*U**T, where
!>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
!>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
!>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
!>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
!>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
!>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
!>
!>             (   I    v    0   )   k-s
!>     U(k) =  (   0    I    0   )   s
!>             (   0    0    I   )   n-k
!>                k-s   s   n-k
!>
!>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
!>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
!>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
!>
!>  If UPLO = 'L', then A = L*D*L**T, where
!>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
!>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
!>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
!>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
!>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
!>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
!>
!>             (   I    0     0   )  k-1
!>     L(k) =  (   0    I     0   )  s
!>             (   0    v     I   )  n-k-s+1
!>                k-1   s  n-k-s+1
!>
!>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
!>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
!>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
!> 

Definition at line 154 of file ssptrf.f.

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