 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ ssytrf_rk()

 subroutine ssytrf_rk ( character UPLO, integer N, real, dimension( lda, * ) A, integer LDA, real, dimension( * ) E, integer, dimension( * ) IPIV, real, dimension( * ) WORK, integer LWORK, integer INFO )

SSYTRF_RK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Kaufman (rook) diagonal pivoting method (BLAS3 blocked algorithm).

Purpose:
``` SSYTRF_RK computes the factorization of a real symmetric matrix A
using the bounded Bunch-Kaufman (rook) diagonal pivoting method:

A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T),

where U (or L) is unit upper (or lower) triangular matrix,
U**T (or L**T) is the transpose of U (or L), P is a permutation
matrix, P**T is the transpose of P, and D is symmetric and block
diagonal with 1-by-1 and 2-by-2 diagonal blocks.

This is the blocked version of the algorithm, calling Level 3 BLAS.
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the symmetric matrix A is stored: = 'U': Upper triangular = 'L': Lower triangular``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in,out] A ``` A is REAL array, dimension (LDA,N) On entry, the symmetric matrix A. If UPLO = 'U': the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L': the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. On exit, contains: a) ONLY diagonal elements of the symmetric block diagonal matrix D on the diagonal of A, i.e. D(k,k) = A(k,k); (superdiagonal (or subdiagonal) elements of D are stored on exit in array E), and b) If UPLO = 'U': factor U in the superdiagonal part of A. If UPLO = 'L': factor L in the subdiagonal part of A.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).``` [out] E ``` E is REAL array, dimension (N) On exit, contains the superdiagonal (or subdiagonal) elements of the symmetric block diagonal matrix D with 1-by-1 or 2-by-2 diagonal blocks, where If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0; If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0. NOTE: For 1-by-1 diagonal block D(k), where 1 <= k <= N, the element E(k) is set to 0 in both UPLO = 'U' or UPLO = 'L' cases.``` [out] IPIV ``` IPIV is INTEGER array, dimension (N) IPIV describes the permutation matrix P in the factorization of matrix A as follows. The absolute value of IPIV(k) represents the index of row and column that were interchanged with the k-th row and column. The value of UPLO describes the order in which the interchanges were applied. Also, the sign of IPIV represents the block structure of the symmetric block diagonal matrix D with 1-by-1 or 2-by-2 diagonal blocks which correspond to 1 or 2 interchanges at each factorization step. For more info see Further Details section. If UPLO = 'U', ( in factorization order, k decreases from N to 1 ): a) A single positive entry IPIV(k) > 0 means: D(k,k) is a 1-by-1 diagonal block. If IPIV(k) != k, rows and columns k and IPIV(k) were interchanged in the matrix A(1:N,1:N); If IPIV(k) = k, no interchange occurred. b) A pair of consecutive negative entries IPIV(k) < 0 and IPIV(k-1) < 0 means: D(k-1:k,k-1:k) is a 2-by-2 diagonal block. (NOTE: negative entries in IPIV appear ONLY in pairs). 1) If -IPIV(k) != k, rows and columns k and -IPIV(k) were interchanged in the matrix A(1:N,1:N). If -IPIV(k) = k, no interchange occurred. 2) If -IPIV(k-1) != k-1, rows and columns k-1 and -IPIV(k-1) were interchanged in the matrix A(1:N,1:N). If -IPIV(k-1) = k-1, no interchange occurred. c) In both cases a) and b), always ABS( IPIV(k) ) <= k. d) NOTE: Any entry IPIV(k) is always NONZERO on output. If UPLO = 'L', ( in factorization order, k increases from 1 to N ): a) A single positive entry IPIV(k) > 0 means: D(k,k) is a 1-by-1 diagonal block. If IPIV(k) != k, rows and columns k and IPIV(k) were interchanged in the matrix A(1:N,1:N). If IPIV(k) = k, no interchange occurred. b) A pair of consecutive negative entries IPIV(k) < 0 and IPIV(k+1) < 0 means: D(k:k+1,k:k+1) is a 2-by-2 diagonal block. (NOTE: negative entries in IPIV appear ONLY in pairs). 1) If -IPIV(k) != k, rows and columns k and -IPIV(k) were interchanged in the matrix A(1:N,1:N). If -IPIV(k) = k, no interchange occurred. 2) If -IPIV(k+1) != k+1, rows and columns k-1 and -IPIV(k-1) were interchanged in the matrix A(1:N,1:N). If -IPIV(k+1) = k+1, no interchange occurred. c) In both cases a) and b), always ABS( IPIV(k) ) >= k. d) NOTE: Any entry IPIV(k) is always NONZERO on output.``` [out] WORK ``` WORK is REAL array, dimension ( MAX(1,LWORK) ). On exit, if INFO = 0, WORK(1) returns the optimal LWORK.``` [in] LWORK ``` LWORK is INTEGER The length of WORK. LWORK >=1. For best performance LWORK >= N*NB, where NB is the block size returned by ILAENV. If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: If INFO = -k, the k-th argument had an illegal value > 0: If INFO = k, the matrix A is singular, because: If UPLO = 'U': column k in the upper triangular part of A contains all zeros. If UPLO = 'L': column k in the lower triangular part of A contains all zeros. Therefore D(k,k) is exactly zero, and superdiagonal elements of column k of U (or subdiagonal elements of column k of L ) are all zeros. 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. NOTE: INFO only stores the first occurrence of a singularity, any subsequent occurrence of singularity is not stored in INFO even though the factorization always completes.```
Further Details:
` TODO: put correct description`
Contributors:
```  December 2016,  Igor Kozachenko,
Computer Science Division,
University of California, Berkeley

September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
School of Mathematics,
University of Manchester```

Definition at line 257 of file ssytrf_rk.f.

259 *
260 * -- LAPACK computational routine --
261 * -- LAPACK is a software package provided by Univ. of Tennessee, --
262 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
263 *
264 * .. Scalar Arguments ..
265  CHARACTER UPLO
266  INTEGER INFO, LDA, LWORK, N
267 * ..
268 * .. Array Arguments ..
269  INTEGER IPIV( * )
270  REAL A( LDA, * ), E( * ), WORK( * )
271 * ..
272 *
273 * =====================================================================
274 *
275 * .. Local Scalars ..
276  LOGICAL LQUERY, UPPER
277  INTEGER I, IINFO, IP, IWS, K, KB, LDWORK, LWKOPT,
278  \$ NB, NBMIN
279 * ..
280 * .. External Functions ..
281  LOGICAL LSAME
282  INTEGER ILAENV
283  EXTERNAL lsame, ilaenv
284 * ..
285 * .. External Subroutines ..
286  EXTERNAL slasyf_rk, ssytf2_rk, sswap, xerbla
287 * ..
288 * .. Intrinsic Functions ..
289  INTRINSIC abs, max
290 * ..
291 * .. Executable Statements ..
292 *
293 * Test the input parameters.
294 *
295  info = 0
296  upper = lsame( uplo, 'U' )
297  lquery = ( lwork.EQ.-1 )
298  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
299  info = -1
300  ELSE IF( n.LT.0 ) THEN
301  info = -2
302  ELSE IF( lda.LT.max( 1, n ) ) THEN
303  info = -4
304  ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
305  info = -8
306  END IF
307 *
308  IF( info.EQ.0 ) THEN
309 *
310 * Determine the block size
311 *
312  nb = ilaenv( 1, 'SSYTRF_RK', uplo, n, -1, -1, -1 )
313  lwkopt = n*nb
314  work( 1 ) = lwkopt
315  END IF
316 *
317  IF( info.NE.0 ) THEN
318  CALL xerbla( 'SSYTRF_RK', -info )
319  RETURN
320  ELSE IF( lquery ) THEN
321  RETURN
322  END IF
323 *
324  nbmin = 2
325  ldwork = n
326  IF( nb.GT.1 .AND. nb.LT.n ) THEN
327  iws = ldwork*nb
328  IF( lwork.LT.iws ) THEN
329  nb = max( lwork / ldwork, 1 )
330  nbmin = max( 2, ilaenv( 2, 'SSYTRF_RK',
331  \$ uplo, n, -1, -1, -1 ) )
332  END IF
333  ELSE
334  iws = 1
335  END IF
336  IF( nb.LT.nbmin )
337  \$ nb = n
338 *
339  IF( upper ) THEN
340 *
341 * Factorize A as U*D*U**T using the upper triangle of A
342 *
343 * K is the main loop index, decreasing from N to 1 in steps of
344 * KB, where KB is the number of columns factorized by SLASYF_RK;
345 * KB is either NB or NB-1, or K for the last block
346 *
347  k = n
348  10 CONTINUE
349 *
350 * If K < 1, exit from loop
351 *
352  IF( k.LT.1 )
353  \$ GO TO 15
354 *
355  IF( k.GT.nb ) THEN
356 *
357 * Factorize columns k-kb+1:k of A and use blocked code to
358 * update columns 1:k-kb
359 *
360  CALL slasyf_rk( uplo, k, nb, kb, a, lda, e,
361  \$ ipiv, work, ldwork, iinfo )
362  ELSE
363 *
364 * Use unblocked code to factorize columns 1:k of A
365 *
366  CALL ssytf2_rk( uplo, k, a, lda, e, ipiv, iinfo )
367  kb = k
368  END IF
369 *
370 * Set INFO on the first occurrence of a zero pivot
371 *
372  IF( info.EQ.0 .AND. iinfo.GT.0 )
373  \$ info = iinfo
374 *
375 * No need to adjust IPIV
376 *
377 *
378 * Apply permutations to the leading panel 1:k-1
379 *
380 * Read IPIV from the last block factored, i.e.
381 * indices k-kb+1:k and apply row permutations to the
382 * last k+1 colunms k+1:N after that block
383 * (We can do the simple loop over IPIV with decrement -1,
384 * since the ABS value of IPIV( I ) represents the row index
385 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
386 *
387  IF( k.LT.n ) THEN
388  DO i = k, ( k - kb + 1 ), -1
389  ip = abs( ipiv( i ) )
390  IF( ip.NE.i ) THEN
391  CALL sswap( n-k, a( i, k+1 ), lda,
392  \$ a( ip, k+1 ), lda )
393  END IF
394  END DO
395  END IF
396 *
397 * Decrease K and return to the start of the main loop
398 *
399  k = k - kb
400  GO TO 10
401 *
402 * This label is the exit from main loop over K decreasing
403 * from N to 1 in steps of KB
404 *
405  15 CONTINUE
406 *
407  ELSE
408 *
409 * Factorize A as L*D*L**T using the lower triangle of A
410 *
411 * K is the main loop index, increasing from 1 to N in steps of
412 * KB, where KB is the number of columns factorized by SLASYF_RK;
413 * KB is either NB or NB-1, or N-K+1 for the last block
414 *
415  k = 1
416  20 CONTINUE
417 *
418 * If K > N, exit from loop
419 *
420  IF( k.GT.n )
421  \$ GO TO 35
422 *
423  IF( k.LE.n-nb ) THEN
424 *
425 * Factorize columns k:k+kb-1 of A and use blocked code to
426 * update columns k+kb:n
427 *
428  CALL slasyf_rk( uplo, n-k+1, nb, kb, a( k, k ), lda, e( k ),
429  \$ ipiv( k ), work, ldwork, iinfo )
430
431
432  ELSE
433 *
434 * Use unblocked code to factorize columns k:n of A
435 *
436  CALL ssytf2_rk( uplo, n-k+1, a( k, k ), lda, e( k ),
437  \$ ipiv( k ), iinfo )
438  kb = n - k + 1
439 *
440  END IF
441 *
442 * Set INFO on the first occurrence of a zero pivot
443 *
444  IF( info.EQ.0 .AND. iinfo.GT.0 )
445  \$ info = iinfo + k - 1
446 *
448 *
449  DO i = k, k + kb - 1
450  IF( ipiv( i ).GT.0 ) THEN
451  ipiv( i ) = ipiv( i ) + k - 1
452  ELSE
453  ipiv( i ) = ipiv( i ) - k + 1
454  END IF
455  END DO
456 *
457 * Apply permutations to the leading panel 1:k-1
458 *
459 * Read IPIV from the last block factored, i.e.
460 * indices k:k+kb-1 and apply row permutations to the
461 * first k-1 colunms 1:k-1 before that block
462 * (We can do the simple loop over IPIV with increment 1,
463 * since the ABS value of IPIV( I ) represents the row index
464 * of the interchange with row i in both 1x1 and 2x2 pivot cases)
465 *
466  IF( k.GT.1 ) THEN
467  DO i = k, ( k + kb - 1 ), 1
468  ip = abs( ipiv( i ) )
469  IF( ip.NE.i ) THEN
470  CALL sswap( k-1, a( i, 1 ), lda,
471  \$ a( ip, 1 ), lda )
472  END IF
473  END DO
474  END IF
475 *
476 * Increase K and return to the start of the main loop
477 *
478  k = k + kb
479  GO TO 20
480 *
481 * This label is the exit from main loop over K increasing
482 * from 1 to N in steps of KB
483 *
484  35 CONTINUE
485 *
486 * End Lower
487 *
488  END IF
489 *
490  work( 1 ) = lwkopt
491  RETURN
492 *
493 * End of SSYTRF_RK
494 *
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine slasyf_rk(UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW, INFO)
SLASYF_RK computes a partial factorization of a real symmetric indefinite matrix using bounded Bunch-...
Definition: slasyf_rk.f:262
subroutine ssytf2_rk(UPLO, N, A, LDA, E, IPIV, INFO)
SSYTF2_RK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Ka...
Definition: ssytf2_rk.f:241
Here is the call graph for this function:
Here is the caller graph for this function: