 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

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

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.
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.```
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*
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: