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

◆ 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).

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

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.
 For more information see Further Details section.
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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 REAL SROUNDUP_LWORK
284 EXTERNAL lsame, ilaenv, sroundup_lwork
285* ..
286* .. External Subroutines ..
287 EXTERNAL slasyf_rk, ssytf2_rk, sswap, xerbla
288* ..
289* .. Intrinsic Functions ..
290 INTRINSIC abs, max
291* ..
292* .. Executable Statements ..
293*
294* Test the input parameters.
295*
296 info = 0
297 upper = lsame( uplo, 'U' )
298 lquery = ( lwork.EQ.-1 )
299 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
300 info = -1
301 ELSE IF( n.LT.0 ) THEN
302 info = -2
303 ELSE IF( lda.LT.max( 1, n ) ) THEN
304 info = -4
305 ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
306 info = -8
307 END IF
308*
309 IF( info.EQ.0 ) THEN
310*
311* Determine the block size
312*
313 nb = ilaenv( 1, 'SSYTRF_RK', uplo, n, -1, -1, -1 )
314 lwkopt = max( 1, n*nb )
315 work( 1 ) = sroundup_lwork(lwkopt)
316 END IF
317*
318 IF( info.NE.0 ) THEN
319 CALL xerbla( 'SSYTRF_RK', -info )
320 RETURN
321 ELSE IF( lquery ) THEN
322 RETURN
323 END IF
324*
325 nbmin = 2
326 ldwork = n
327 IF( nb.GT.1 .AND. nb.LT.n ) THEN
328 iws = ldwork*nb
329 IF( lwork.LT.iws ) THEN
330 nb = max( lwork / ldwork, 1 )
331 nbmin = max( 2, ilaenv( 2, 'SSYTRF_RK',
332 $ uplo, n, -1, -1, -1 ) )
333 END IF
334 ELSE
335 iws = 1
336 END IF
337 IF( nb.LT.nbmin )
338 $ nb = n
339*
340 IF( upper ) THEN
341*
342* Factorize A as U*D*U**T using the upper triangle of A
343*
344* K is the main loop index, decreasing from N to 1 in steps of
345* KB, where KB is the number of columns factorized by SLASYF_RK;
346* KB is either NB or NB-1, or K for the last block
347*
348 k = n
349 10 CONTINUE
350*
351* If K < 1, exit from loop
352*
353 IF( k.LT.1 )
354 $ GO TO 15
355*
356 IF( k.GT.nb ) THEN
357*
358* Factorize columns k-kb+1:k of A and use blocked code to
359* update columns 1:k-kb
360*
361 CALL slasyf_rk( uplo, k, nb, kb, a, lda, e,
362 $ ipiv, work, ldwork, iinfo )
363 ELSE
364*
365* Use unblocked code to factorize columns 1:k of A
366*
367 CALL ssytf2_rk( uplo, k, a, lda, e, ipiv, iinfo )
368 kb = k
369 END IF
370*
371* Set INFO on the first occurrence of a zero pivot
372*
373 IF( info.EQ.0 .AND. iinfo.GT.0 )
374 $ info = iinfo
375*
376* No need to adjust IPIV
377*
378*
379* Apply permutations to the leading panel 1:k-1
380*
381* Read IPIV from the last block factored, i.e.
382* indices k-kb+1:k and apply row permutations to the
383* last k+1 colunms k+1:N after that block
384* (We can do the simple loop over IPIV with decrement -1,
385* since the ABS value of IPIV( I ) represents the row index
386* of the interchange with row i in both 1x1 and 2x2 pivot cases)
387*
388 IF( k.LT.n ) THEN
389 DO i = k, ( k - kb + 1 ), -1
390 ip = abs( ipiv( i ) )
391 IF( ip.NE.i ) THEN
392 CALL sswap( n-k, a( i, k+1 ), lda,
393 $ a( ip, k+1 ), lda )
394 END IF
395 END DO
396 END IF
397*
398* Decrease K and return to the start of the main loop
399*
400 k = k - kb
401 GO TO 10
402*
403* This label is the exit from main loop over K decreasing
404* from N to 1 in steps of KB
405*
406 15 CONTINUE
407*
408 ELSE
409*
410* Factorize A as L*D*L**T using the lower triangle of A
411*
412* K is the main loop index, increasing from 1 to N in steps of
413* KB, where KB is the number of columns factorized by SLASYF_RK;
414* KB is either NB or NB-1, or N-K+1 for the last block
415*
416 k = 1
417 20 CONTINUE
418*
419* If K > N, exit from loop
420*
421 IF( k.GT.n )
422 $ GO TO 35
423*
424 IF( k.LE.n-nb ) THEN
425*
426* Factorize columns k:k+kb-1 of A and use blocked code to
427* update columns k+kb:n
428*
429 CALL slasyf_rk( uplo, n-k+1, nb, kb, a( k, k ), lda, e( k ),
430 $ ipiv( k ), work, ldwork, iinfo )
431
432
433 ELSE
434*
435* Use unblocked code to factorize columns k:n of A
436*
437 CALL ssytf2_rk( uplo, n-k+1, a( k, k ), lda, e( k ),
438 $ ipiv( k ), iinfo )
439 kb = n - k + 1
440*
441 END IF
442*
443* Set INFO on the first occurrence of a zero pivot
444*
445 IF( info.EQ.0 .AND. iinfo.GT.0 )
446 $ info = iinfo + k - 1
447*
448* Adjust IPIV
449*
450 DO i = k, k + kb - 1
451 IF( ipiv( i ).GT.0 ) THEN
452 ipiv( i ) = ipiv( i ) + k - 1
453 ELSE
454 ipiv( i ) = ipiv( i ) - k + 1
455 END IF
456 END DO
457*
458* Apply permutations to the leading panel 1:k-1
459*
460* Read IPIV from the last block factored, i.e.
461* indices k:k+kb-1 and apply row permutations to the
462* first k-1 colunms 1:k-1 before that block
463* (We can do the simple loop over IPIV with increment 1,
464* since the ABS value of IPIV( I ) represents the row index
465* of the interchange with row i in both 1x1 and 2x2 pivot cases)
466*
467 IF( k.GT.1 ) THEN
468 DO i = k, ( k + kb - 1 ), 1
469 ip = abs( ipiv( i ) )
470 IF( ip.NE.i ) THEN
471 CALL sswap( k-1, a( i, 1 ), lda,
472 $ a( ip, 1 ), lda )
473 END IF
474 END DO
475 END IF
476*
477* Increase K and return to the start of the main loop
478*
479 k = k + kb
480 GO TO 20
481*
482* This label is the exit from main loop over K increasing
483* from 1 to N in steps of KB
484*
485 35 CONTINUE
486*
487* End Lower
488*
489 END IF
490*
491 work( 1 ) = sroundup_lwork(lwkopt)
492 RETURN
493*
494* End of SSYTRF_RK
495*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
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: