LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
zhetrf_rk.f
Go to the documentation of this file.
1 *> \brief \b ZHETRF_RK computes the factorization of a complex Hermitian indefinite matrix using the bounded Bunch-Kaufman (rook) diagonal pivoting method (BLAS3 blocked algorithm).
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZHETRF_RK + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrf_rk.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrf_rk.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrf_rk.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHETRF_RK( UPLO, N, A, LDA, E, IPIV, WORK, LWORK,
22 * INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER UPLO
26 * INTEGER INFO, LDA, LWORK, N
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IPIV( * )
30 * COMPLEX*16 A( LDA, * ), E ( * ), WORK( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *> ZHETRF_RK computes the factorization of a complex Hermitian matrix A
39 *> using the bounded Bunch-Kaufman (rook) diagonal pivoting method:
40 *>
41 *> A = P*U*D*(U**H)*(P**T) or A = P*L*D*(L**H)*(P**T),
42 *>
43 *> where U (or L) is unit upper (or lower) triangular matrix,
44 *> U**H (or L**H) is the conjugate of U (or L), P is a permutation
45 *> matrix, P**T is the transpose of P, and D is Hermitian and block
46 *> diagonal with 1-by-1 and 2-by-2 diagonal blocks.
47 *>
48 *> This is the blocked version of the algorithm, calling Level 3 BLAS.
49 *> For more information see Further Details section.
50 *> \endverbatim
51 *
52 * Arguments:
53 * ==========
54 *
55 *> \param[in] UPLO
56 *> \verbatim
57 *> UPLO is CHARACTER*1
58 *> Specifies whether the upper or lower triangular part of the
59 *> Hermitian matrix A is stored:
60 *> = 'U': Upper triangular
61 *> = 'L': Lower triangular
62 *> \endverbatim
63 *>
64 *> \param[in] N
65 *> \verbatim
66 *> N is INTEGER
67 *> The order of the matrix A. N >= 0.
68 *> \endverbatim
69 *>
70 *> \param[in,out] A
71 *> \verbatim
72 *> A is COMPLEX*16 array, dimension (LDA,N)
73 *> On entry, the Hermitian matrix A.
74 *> If UPLO = 'U': the leading N-by-N upper triangular part
75 *> of A contains the upper triangular part of the matrix A,
76 *> and the strictly lower triangular part of A is not
77 *> referenced.
78 *>
79 *> If UPLO = 'L': the leading N-by-N lower triangular part
80 *> of A contains the lower triangular part of the matrix A,
81 *> and the strictly upper triangular part of A is not
82 *> referenced.
83 *>
84 *> On exit, contains:
85 *> a) ONLY diagonal elements of the Hermitian block diagonal
86 *> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
87 *> (superdiagonal (or subdiagonal) elements of D
88 *> are stored on exit in array E), and
89 *> b) If UPLO = 'U': factor U in the superdiagonal part of A.
90 *> If UPLO = 'L': factor L in the subdiagonal part of A.
91 *> \endverbatim
92 *>
93 *> \param[in] LDA
94 *> \verbatim
95 *> LDA is INTEGER
96 *> The leading dimension of the array A. LDA >= max(1,N).
97 *> \endverbatim
98 *>
99 *> \param[out] E
100 *> \verbatim
101 *> E is COMPLEX*16 array, dimension (N)
102 *> On exit, contains the superdiagonal (or subdiagonal)
103 *> elements of the Hermitian block diagonal matrix D
104 *> with 1-by-1 or 2-by-2 diagonal blocks, where
105 *> If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
106 *> If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.
107 *>
108 *> NOTE: For 1-by-1 diagonal block D(k), where
109 *> 1 <= k <= N, the element E(k) is set to 0 in both
110 *> UPLO = 'U' or UPLO = 'L' cases.
111 *> \endverbatim
112 *>
113 *> \param[out] IPIV
114 *> \verbatim
115 *> IPIV is INTEGER array, dimension (N)
116 *> IPIV describes the permutation matrix P in the factorization
117 *> of matrix A as follows. The absolute value of IPIV(k)
118 *> represents the index of row and column that were
119 *> interchanged with the k-th row and column. The value of UPLO
120 *> describes the order in which the interchanges were applied.
121 *> Also, the sign of IPIV represents the block structure of
122 *> the Hermitian block diagonal matrix D with 1-by-1 or 2-by-2
123 *> diagonal blocks which correspond to 1 or 2 interchanges
124 *> at each factorization step. For more info see Further
125 *> Details section.
126 *>
127 *> If UPLO = 'U',
128 *> ( in factorization order, k decreases from N to 1 ):
129 *> a) A single positive entry IPIV(k) > 0 means:
130 *> D(k,k) is a 1-by-1 diagonal block.
131 *> If IPIV(k) != k, rows and columns k and IPIV(k) were
132 *> interchanged in the matrix A(1:N,1:N);
133 *> If IPIV(k) = k, no interchange occurred.
134 *>
135 *> b) A pair of consecutive negative entries
136 *> IPIV(k) < 0 and IPIV(k-1) < 0 means:
137 *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
138 *> (NOTE: negative entries in IPIV appear ONLY in pairs).
139 *> 1) If -IPIV(k) != k, rows and columns
140 *> k and -IPIV(k) were interchanged
141 *> in the matrix A(1:N,1:N).
142 *> If -IPIV(k) = k, no interchange occurred.
143 *> 2) If -IPIV(k-1) != k-1, rows and columns
144 *> k-1 and -IPIV(k-1) were interchanged
145 *> in the matrix A(1:N,1:N).
146 *> If -IPIV(k-1) = k-1, no interchange occurred.
147 *>
148 *> c) In both cases a) and b), always ABS( IPIV(k) ) <= k.
149 *>
150 *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
151 *>
152 *> If UPLO = 'L',
153 *> ( in factorization order, k increases from 1 to N ):
154 *> a) A single positive entry IPIV(k) > 0 means:
155 *> D(k,k) is a 1-by-1 diagonal block.
156 *> If IPIV(k) != k, rows and columns k and IPIV(k) were
157 *> interchanged in the matrix A(1:N,1:N).
158 *> If IPIV(k) = k, no interchange occurred.
159 *>
160 *> b) A pair of consecutive negative entries
161 *> IPIV(k) < 0 and IPIV(k+1) < 0 means:
162 *> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
163 *> (NOTE: negative entries in IPIV appear ONLY in pairs).
164 *> 1) If -IPIV(k) != k, rows and columns
165 *> k and -IPIV(k) were interchanged
166 *> in the matrix A(1:N,1:N).
167 *> If -IPIV(k) = k, no interchange occurred.
168 *> 2) If -IPIV(k+1) != k+1, rows and columns
169 *> k-1 and -IPIV(k-1) were interchanged
170 *> in the matrix A(1:N,1:N).
171 *> If -IPIV(k+1) = k+1, no interchange occurred.
172 *>
173 *> c) In both cases a) and b), always ABS( IPIV(k) ) >= k.
174 *>
175 *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
176 *> \endverbatim
177 *>
178 *> \param[out] WORK
179 *> \verbatim
180 *> WORK is COMPLEX*16 array, dimension ( MAX(1,LWORK) ).
181 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
182 *> \endverbatim
183 *>
184 *> \param[in] LWORK
185 *> \verbatim
186 *> LWORK is INTEGER
187 *> The length of WORK. LWORK >=1. For best performance
188 *> LWORK >= N*NB, where NB is the block size returned
189 *> by ILAENV.
190 *>
191 *> If LWORK = -1, then a workspace query is assumed;
192 *> the routine only calculates the optimal size of the WORK
193 *> array, returns this value as the first entry of the WORK
194 *> array, and no error message related to LWORK is issued
195 *> by XERBLA.
196 *> \endverbatim
197 *>
198 *> \param[out] INFO
199 *> \verbatim
200 *> INFO is INTEGER
201 *> = 0: successful exit
202 *>
203 *> < 0: If INFO = -k, the k-th argument had an illegal value
204 *>
205 *> > 0: If INFO = k, the matrix A is singular, because:
206 *> If UPLO = 'U': column k in the upper
207 *> triangular part of A contains all zeros.
208 *> If UPLO = 'L': column k in the lower
209 *> triangular part of A contains all zeros.
210 *>
211 *> Therefore D(k,k) is exactly zero, and superdiagonal
212 *> elements of column k of U (or subdiagonal elements of
213 *> column k of L ) are all zeros. The factorization has
214 *> been completed, but the block diagonal matrix D is
215 *> exactly singular, and division by zero will occur if
216 *> it is used to solve a system of equations.
217 *>
218 *> NOTE: INFO only stores the first occurrence of
219 *> a singularity, any subsequent occurrence of singularity
220 *> is not stored in INFO even though the factorization
221 *> always completes.
222 *> \endverbatim
223 *
224 * Authors:
225 * ========
226 *
227 *> \author Univ. of Tennessee
228 *> \author Univ. of California Berkeley
229 *> \author Univ. of Colorado Denver
230 *> \author NAG Ltd.
231 *
232 *> \ingroup complex16HEcomputational
233 *
234 *> \par Further Details:
235 * =====================
236 *>
237 *> \verbatim
238 *> TODO: put correct description
239 *> \endverbatim
240 *
241 *> \par Contributors:
242 * ==================
243 *>
244 *> \verbatim
245 *>
246 *> December 2016, Igor Kozachenko,
247 *> Computer Science Division,
248 *> University of California, Berkeley
249 *>
250 *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
251 *> School of Mathematics,
252 *> University of Manchester
253 *>
254 *> \endverbatim
255 *
256 * =====================================================================
257  SUBROUTINE zhetrf_rk( UPLO, N, A, LDA, E, IPIV, WORK, LWORK,
258  $ INFO )
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 zlahef_rk, zhetf2_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, 'ZHETRF_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( 'ZHETRF_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, 'ZHETRF_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 ZLAHEF_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 zlahef_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 zhetf2_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 ZLAHEF_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 zlahef_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 zhetf2_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 ZHETRF_RK
494 *
495  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zlahef_rk(UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW, INFO)
ZLAHEF_RK computes a partial factorization of a complex Hermitian indefinite matrix using bounded Bun...
Definition: zlahef_rk.f:262
subroutine zhetrf_rk(UPLO, N, A, LDA, E, IPIV, WORK, LWORK, INFO)
ZHETRF_RK computes the factorization of a complex Hermitian indefinite matrix using the bounded Bunch...
Definition: zhetrf_rk.f:259
subroutine zhetf2_rk(UPLO, N, A, LDA, E, IPIV, INFO)
ZHETF2_RK computes the factorization of a complex Hermitian indefinite matrix using the bounded Bunch...
Definition: zhetf2_rk.f:241