LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zsytrf_rk()

subroutine zsytrf_rk ( character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  E,
integer, dimension( * )  IPIV,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

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

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

Purpose:
 ZSYTRF_RK computes the factorization of a complex 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 zsytrf_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  COMPLEX*16 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 zlasyf_rk, zsytf2_rk, zswap, 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, 'ZSYTRF_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( 'ZSYTRF_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, 'ZSYTRF_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 ZLASYF_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 zlasyf_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 zsytf2_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 zswap( 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 ZLASYF_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 zlasyf_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 zsytf2_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 *
447 * Adjust IPIV
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 zswap( 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 ZSYTRF_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 zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zlasyf_rk(UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW, INFO)
ZLASYF_RK computes a partial factorization of a complex symmetric indefinite matrix using bounded Bun...
Definition: zlasyf_rk.f:262
subroutine zsytf2_rk(UPLO, N, A, LDA, E, IPIV, INFO)
ZSYTF2_RK computes the factorization of a complex symmetric indefinite matrix using the bounded Bunch...
Definition: zsytf2_rk.f:241
Here is the call graph for this function:
Here is the caller graph for this function: