LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dsytrf_rk()

subroutine dsytrf_rk ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  E,
integer, dimension( * )  IPIV,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

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

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

Purpose:
 DSYTRF_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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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.
Date
December 2016
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 261 of file dsytrf_rk.f.

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