LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
zhetf2_rk.f
Go to the documentation of this file.
1 *> \brief \b ZHETF2_RK computes the factorization of a complex Hermitian indefinite matrix using the bounded Bunch-Kaufman (rook) diagonal pivoting method (BLAS2 unblocked algorithm).
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZHETF2_RK + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetf2_rk.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetf2_rk.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetf2_rk.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHETF2_RK( UPLO, N, A, LDA, E, IPIV, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER UPLO
25 * INTEGER INFO, LDA, N
26 * ..
27 * .. Array Arguments ..
28 * INTEGER IPIV( * )
29 * COMPLEX*16 A( LDA, * ), E ( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *> ZHETF2_RK computes the factorization of a complex Hermitian matrix A
38 *> using the bounded Bunch-Kaufman (rook) diagonal pivoting method:
39 *>
40 *> A = P*U*D*(U**H)*(P**T) or A = P*L*D*(L**H)*(P**T),
41 *>
42 *> where U (or L) is unit upper (or lower) triangular matrix,
43 *> U**H (or L**H) is the conjugate of U (or L), P is a permutation
44 *> matrix, P**T is the transpose of P, and D is Hermitian and block
45 *> diagonal with 1-by-1 and 2-by-2 diagonal blocks.
46 *>
47 *> This is the unblocked version of the algorithm, calling Level 2 BLAS.
48 *> For more information see Further Details section.
49 *> \endverbatim
50 *
51 * Arguments:
52 * ==========
53 *
54 *> \param[in] UPLO
55 *> \verbatim
56 *> UPLO is CHARACTER*1
57 *> Specifies whether the upper or lower triangular part of the
58 *> Hermitian matrix A is stored:
59 *> = 'U': Upper triangular
60 *> = 'L': Lower triangular
61 *> \endverbatim
62 *>
63 *> \param[in] N
64 *> \verbatim
65 *> N is INTEGER
66 *> The order of the matrix A. N >= 0.
67 *> \endverbatim
68 *>
69 *> \param[in,out] A
70 *> \verbatim
71 *> A is COMPLEX*16 array, dimension (LDA,N)
72 *> On entry, the Hermitian matrix A.
73 *> If UPLO = 'U': the leading N-by-N upper triangular part
74 *> of A contains the upper triangular part of the matrix A,
75 *> and the strictly lower triangular part of A is not
76 *> referenced.
77 *>
78 *> If UPLO = 'L': the leading N-by-N lower triangular part
79 *> of A contains the lower triangular part of the matrix A,
80 *> and the strictly upper triangular part of A is not
81 *> referenced.
82 *>
83 *> On exit, contains:
84 *> a) ONLY diagonal elements of the Hermitian block diagonal
85 *> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
86 *> (superdiagonal (or subdiagonal) elements of D
87 *> are stored on exit in array E), and
88 *> b) If UPLO = 'U': factor U in the superdiagonal part of A.
89 *> If UPLO = 'L': factor L in the subdiagonal part of A.
90 *> \endverbatim
91 *>
92 *> \param[in] LDA
93 *> \verbatim
94 *> LDA is INTEGER
95 *> The leading dimension of the array A. LDA >= max(1,N).
96 *> \endverbatim
97 *>
98 *> \param[out] E
99 *> \verbatim
100 *> E is COMPLEX*16 array, dimension (N)
101 *> On exit, contains the superdiagonal (or subdiagonal)
102 *> elements of the Hermitian block diagonal matrix D
103 *> with 1-by-1 or 2-by-2 diagonal blocks, where
104 *> If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
105 *> If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.
106 *>
107 *> NOTE: For 1-by-1 diagonal block D(k), where
108 *> 1 <= k <= N, the element E(k) is set to 0 in both
109 *> UPLO = 'U' or UPLO = 'L' cases.
110 *> \endverbatim
111 *>
112 *> \param[out] IPIV
113 *> \verbatim
114 *> IPIV is INTEGER array, dimension (N)
115 *> IPIV describes the permutation matrix P in the factorization
116 *> of matrix A as follows. The absolute value of IPIV(k)
117 *> represents the index of row and column that were
118 *> interchanged with the k-th row and column. The value of UPLO
119 *> describes the order in which the interchanges were applied.
120 *> Also, the sign of IPIV represents the block structure of
121 *> the Hermitian block diagonal matrix D with 1-by-1 or 2-by-2
122 *> diagonal blocks which correspond to 1 or 2 interchanges
123 *> at each factorization step. For more info see Further
124 *> Details section.
125 *>
126 *> If UPLO = 'U',
127 *> ( in factorization order, k decreases from N to 1 ):
128 *> a) A single positive entry IPIV(k) > 0 means:
129 *> D(k,k) is a 1-by-1 diagonal block.
130 *> If IPIV(k) != k, rows and columns k and IPIV(k) were
131 *> interchanged in the matrix A(1:N,1:N);
132 *> If IPIV(k) = k, no interchange occurred.
133 *>
134 *> b) A pair of consecutive negative entries
135 *> IPIV(k) < 0 and IPIV(k-1) < 0 means:
136 *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
137 *> (NOTE: negative entries in IPIV appear ONLY in pairs).
138 *> 1) If -IPIV(k) != k, rows and columns
139 *> k and -IPIV(k) were interchanged
140 *> in the matrix A(1:N,1:N).
141 *> If -IPIV(k) = k, no interchange occurred.
142 *> 2) If -IPIV(k-1) != k-1, rows and columns
143 *> k-1 and -IPIV(k-1) were interchanged
144 *> in the matrix A(1:N,1:N).
145 *> If -IPIV(k-1) = k-1, no interchange occurred.
146 *>
147 *> c) In both cases a) and b), always ABS( IPIV(k) ) <= k.
148 *>
149 *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
150 *>
151 *> If UPLO = 'L',
152 *> ( in factorization order, k increases from 1 to N ):
153 *> a) A single positive entry IPIV(k) > 0 means:
154 *> D(k,k) is a 1-by-1 diagonal block.
155 *> If IPIV(k) != k, rows and columns k and IPIV(k) were
156 *> interchanged in the matrix A(1:N,1:N).
157 *> If IPIV(k) = k, no interchange occurred.
158 *>
159 *> b) A pair of consecutive negative entries
160 *> IPIV(k) < 0 and IPIV(k+1) < 0 means:
161 *> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
162 *> (NOTE: negative entries in IPIV appear ONLY in pairs).
163 *> 1) If -IPIV(k) != k, rows and columns
164 *> k and -IPIV(k) were interchanged
165 *> in the matrix A(1:N,1:N).
166 *> If -IPIV(k) = k, no interchange occurred.
167 *> 2) If -IPIV(k+1) != k+1, rows and columns
168 *> k-1 and -IPIV(k-1) were interchanged
169 *> in the matrix A(1:N,1:N).
170 *> If -IPIV(k+1) = k+1, no interchange occurred.
171 *>
172 *> c) In both cases a) and b), always ABS( IPIV(k) ) >= k.
173 *>
174 *> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
175 *> \endverbatim
176 *>
177 *> \param[out] INFO
178 *> \verbatim
179 *> INFO is INTEGER
180 *> = 0: successful exit
181 *>
182 *> < 0: If INFO = -k, the k-th argument had an illegal value
183 *>
184 *> > 0: If INFO = k, the matrix A is singular, because:
185 *> If UPLO = 'U': column k in the upper
186 *> triangular part of A contains all zeros.
187 *> If UPLO = 'L': column k in the lower
188 *> triangular part of A contains all zeros.
189 *>
190 *> Therefore D(k,k) is exactly zero, and superdiagonal
191 *> elements of column k of U (or subdiagonal elements of
192 *> column k of L ) are all zeros. The factorization has
193 *> been completed, but the block diagonal matrix D is
194 *> exactly singular, and division by zero will occur if
195 *> it is used to solve a system of equations.
196 *>
197 *> NOTE: INFO only stores the first occurrence of
198 *> a singularity, any subsequent occurrence of singularity
199 *> is not stored in INFO even though the factorization
200 *> always completes.
201 *> \endverbatim
202 *
203 * Authors:
204 * ========
205 *
206 *> \author Univ. of Tennessee
207 *> \author Univ. of California Berkeley
208 *> \author Univ. of Colorado Denver
209 *> \author NAG Ltd.
210 *
211 *> \ingroup complex16HEcomputational
212 *
213 *> \par Further Details:
214 * =====================
215 *>
216 *> \verbatim
217 *> TODO: put further details
218 *> \endverbatim
219 *
220 *> \par Contributors:
221 * ==================
222 *>
223 *> \verbatim
224 *>
225 *> December 2016, Igor Kozachenko,
226 *> Computer Science Division,
227 *> University of California, Berkeley
228 *>
229 *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
230 *> School of Mathematics,
231 *> University of Manchester
232 *>
233 *> 01-01-96 - Based on modifications by
234 *> J. Lewis, Boeing Computer Services Company
235 *> A. Petitet, Computer Science Dept.,
236 *> Univ. of Tenn., Knoxville abd , USA
237 *> \endverbatim
238 *
239 * =====================================================================
240  SUBROUTINE zhetf2_rk( UPLO, N, A, LDA, E, IPIV, INFO )
241 *
242 * -- LAPACK computational routine --
243 * -- LAPACK is a software package provided by Univ. of Tennessee, --
244 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
245 *
246 * .. Scalar Arguments ..
247  CHARACTER UPLO
248  INTEGER INFO, LDA, N
249 * ..
250 * .. Array Arguments ..
251  INTEGER IPIV( * )
252  COMPLEX*16 A( LDA, * ), E( * )
253 * ..
254 *
255 * ======================================================================
256 *
257 * .. Parameters ..
258  DOUBLE PRECISION ZERO, ONE
259  parameter( zero = 0.0d+0, one = 1.0d+0 )
260  DOUBLE PRECISION EIGHT, SEVTEN
261  parameter( eight = 8.0d+0, sevten = 17.0d+0 )
262  COMPLEX*16 CZERO
263  parameter( czero = ( 0.0d+0, 0.0d+0 ) )
264 * ..
265 * .. Local Scalars ..
266  LOGICAL DONE, UPPER
267  INTEGER I, II, IMAX, ITEMP, J, JMAX, K, KK, KP, KSTEP,
268  $ P
269  DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, DTEMP,
270  $ ROWMAX, TT, SFMIN
271  COMPLEX*16 D12, D21, T, WK, WKM1, WKP1, Z
272 * ..
273 * .. External Functions ..
274 *
275  LOGICAL LSAME
276  INTEGER IZAMAX
277  DOUBLE PRECISION DLAMCH, DLAPY2
278  EXTERNAL lsame, izamax, dlamch, dlapy2
279 * ..
280 * .. External Subroutines ..
281  EXTERNAL xerbla, zdscal, zher, zswap
282 * ..
283 * .. Intrinsic Functions ..
284  INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, sqrt
285 * ..
286 * .. Statement Functions ..
287  DOUBLE PRECISION CABS1
288 * ..
289 * .. Statement Function definitions ..
290  cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
291 * ..
292 * .. Executable Statements ..
293 *
294 * Test the input parameters.
295 *
296  info = 0
297  upper = lsame( uplo, 'U' )
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  END IF
305  IF( info.NE.0 ) THEN
306  CALL xerbla( 'ZHETF2_RK', -info )
307  RETURN
308  END IF
309 *
310 * Initialize ALPHA for use in choosing pivot block size.
311 *
312  alpha = ( one+sqrt( sevten ) ) / eight
313 *
314 * Compute machine safe minimum
315 *
316  sfmin = dlamch( 'S' )
317 *
318  IF( upper ) THEN
319 *
320 * Factorize A as U*D*U**H using the upper triangle of A
321 *
322 * Initialize the first entry of array E, where superdiagonal
323 * elements of D are stored
324 *
325  e( 1 ) = czero
326 *
327 * K is the main loop index, decreasing from N to 1 in steps of
328 * 1 or 2
329 *
330  k = n
331  10 CONTINUE
332 *
333 * If K < 1, exit from loop
334 *
335  IF( k.LT.1 )
336  $ GO TO 34
337  kstep = 1
338  p = k
339 *
340 * Determine rows and columns to be interchanged and whether
341 * a 1-by-1 or 2-by-2 pivot block will be used
342 *
343  absakk = abs( dble( a( k, k ) ) )
344 *
345 * IMAX is the row-index of the largest off-diagonal element in
346 * column K, and COLMAX is its absolute value.
347 * Determine both COLMAX and IMAX.
348 *
349  IF( k.GT.1 ) THEN
350  imax = izamax( k-1, a( 1, k ), 1 )
351  colmax = cabs1( a( imax, k ) )
352  ELSE
353  colmax = zero
354  END IF
355 *
356  IF( ( max( absakk, colmax ).EQ.zero ) ) THEN
357 *
358 * Column K is zero or underflow: set INFO and continue
359 *
360  IF( info.EQ.0 )
361  $ info = k
362  kp = k
363  a( k, k ) = dble( a( k, k ) )
364 *
365 * Set E( K ) to zero
366 *
367  IF( k.GT.1 )
368  $ e( k ) = czero
369 *
370  ELSE
371 *
372 * ============================================================
373 *
374 * BEGIN pivot search
375 *
376 * Case(1)
377 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
378 * (used to handle NaN and Inf)
379 *
380  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
381 *
382 * no interchange, use 1-by-1 pivot block
383 *
384  kp = k
385 *
386  ELSE
387 *
388  done = .false.
389 *
390 * Loop until pivot found
391 *
392  12 CONTINUE
393 *
394 * BEGIN pivot search loop body
395 *
396 *
397 * JMAX is the column-index of the largest off-diagonal
398 * element in row IMAX, and ROWMAX is its absolute value.
399 * Determine both ROWMAX and JMAX.
400 *
401  IF( imax.NE.k ) THEN
402  jmax = imax + izamax( k-imax, a( imax, imax+1 ),
403  $ lda )
404  rowmax = cabs1( a( imax, jmax ) )
405  ELSE
406  rowmax = zero
407  END IF
408 *
409  IF( imax.GT.1 ) THEN
410  itemp = izamax( imax-1, a( 1, imax ), 1 )
411  dtemp = cabs1( a( itemp, imax ) )
412  IF( dtemp.GT.rowmax ) THEN
413  rowmax = dtemp
414  jmax = itemp
415  END IF
416  END IF
417 *
418 * Case(2)
419 * Equivalent to testing for
420 * ABS( DBLE( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
421 * (used to handle NaN and Inf)
422 *
423  IF( .NOT.( abs( dble( a( imax, imax ) ) )
424  $ .LT.alpha*rowmax ) ) THEN
425 *
426 * interchange rows and columns K and IMAX,
427 * use 1-by-1 pivot block
428 *
429  kp = imax
430  done = .true.
431 *
432 * Case(3)
433 * Equivalent to testing for ROWMAX.EQ.COLMAX,
434 * (used to handle NaN and Inf)
435 *
436  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
437  $ THEN
438 *
439 * interchange rows and columns K-1 and IMAX,
440 * use 2-by-2 pivot block
441 *
442  kp = imax
443  kstep = 2
444  done = .true.
445 *
446 * Case(4)
447  ELSE
448 *
449 * Pivot not found: set params and repeat
450 *
451  p = imax
452  colmax = rowmax
453  imax = jmax
454  END IF
455 *
456 * END pivot search loop body
457 *
458  IF( .NOT.done ) GOTO 12
459 *
460  END IF
461 *
462 * END pivot search
463 *
464 * ============================================================
465 *
466 * KK is the column of A where pivoting step stopped
467 *
468  kk = k - kstep + 1
469 *
470 * For only a 2x2 pivot, interchange rows and columns K and P
471 * in the leading submatrix A(1:k,1:k)
472 *
473  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
474 * (1) Swap columnar parts
475  IF( p.GT.1 )
476  $ CALL zswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
477 * (2) Swap and conjugate middle parts
478  DO 14 j = p + 1, k - 1
479  t = dconjg( a( j, k ) )
480  a( j, k ) = dconjg( a( p, j ) )
481  a( p, j ) = t
482  14 CONTINUE
483 * (3) Swap and conjugate corner elements at row-col interserction
484  a( p, k ) = dconjg( a( p, k ) )
485 * (4) Swap diagonal elements at row-col intersection
486  r1 = dble( a( k, k ) )
487  a( k, k ) = dble( a( p, p ) )
488  a( p, p ) = r1
489 *
490 * Convert upper triangle of A into U form by applying
491 * the interchanges in columns k+1:N.
492 *
493  IF( k.LT.n )
494  $ CALL zswap( n-k, a( k, k+1 ), lda, a( p, k+1 ), lda )
495 *
496  END IF
497 *
498 * For both 1x1 and 2x2 pivots, interchange rows and
499 * columns KK and KP in the leading submatrix A(1:k,1:k)
500 *
501  IF( kp.NE.kk ) THEN
502 * (1) Swap columnar parts
503  IF( kp.GT.1 )
504  $ CALL zswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
505 * (2) Swap and conjugate middle parts
506  DO 15 j = kp + 1, kk - 1
507  t = dconjg( a( j, kk ) )
508  a( j, kk ) = dconjg( a( kp, j ) )
509  a( kp, j ) = t
510  15 CONTINUE
511 * (3) Swap and conjugate corner elements at row-col interserction
512  a( kp, kk ) = dconjg( a( kp, kk ) )
513 * (4) Swap diagonal elements at row-col intersection
514  r1 = dble( a( kk, kk ) )
515  a( kk, kk ) = dble( a( kp, kp ) )
516  a( kp, kp ) = r1
517 *
518  IF( kstep.EQ.2 ) THEN
519 * (*) Make sure that diagonal element of pivot is real
520  a( k, k ) = dble( a( k, k ) )
521 * (5) Swap row elements
522  t = a( k-1, k )
523  a( k-1, k ) = a( kp, k )
524  a( kp, k ) = t
525  END IF
526 *
527 * Convert upper triangle of A into U form by applying
528 * the interchanges in columns k+1:N.
529 *
530  IF( k.LT.n )
531  $ CALL zswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
532  $ lda )
533 *
534  ELSE
535 * (*) Make sure that diagonal element of pivot is real
536  a( k, k ) = dble( a( k, k ) )
537  IF( kstep.EQ.2 )
538  $ a( k-1, k-1 ) = dble( a( k-1, k-1 ) )
539  END IF
540 *
541 * Update the leading submatrix
542 *
543  IF( kstep.EQ.1 ) THEN
544 *
545 * 1-by-1 pivot block D(k): column k now holds
546 *
547 * W(k) = U(k)*D(k)
548 *
549 * where U(k) is the k-th column of U
550 *
551  IF( k.GT.1 ) THEN
552 *
553 * Perform a rank-1 update of A(1:k-1,1:k-1) and
554 * store U(k) in column k
555 *
556  IF( abs( dble( a( k, k ) ) ).GE.sfmin ) THEN
557 *
558 * Perform a rank-1 update of A(1:k-1,1:k-1) as
559 * A := A - U(k)*D(k)*U(k)**T
560 * = A - W(k)*1/D(k)*W(k)**T
561 *
562  d11 = one / dble( a( k, k ) )
563  CALL zher( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
564 *
565 * Store U(k) in column k
566 *
567  CALL zdscal( k-1, d11, a( 1, k ), 1 )
568  ELSE
569 *
570 * Store L(k) in column K
571 *
572  d11 = dble( a( k, k ) )
573  DO 16 ii = 1, k - 1
574  a( ii, k ) = a( ii, k ) / d11
575  16 CONTINUE
576 *
577 * Perform a rank-1 update of A(k+1:n,k+1:n) as
578 * A := A - U(k)*D(k)*U(k)**T
579 * = A - W(k)*(1/D(k))*W(k)**T
580 * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
581 *
582  CALL zher( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
583  END IF
584 *
585 * Store the superdiagonal element of D in array E
586 *
587  e( k ) = czero
588 *
589  END IF
590 *
591  ELSE
592 *
593 * 2-by-2 pivot block D(k): columns k and k-1 now hold
594 *
595 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
596 *
597 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
598 * of U
599 *
600 * Perform a rank-2 update of A(1:k-2,1:k-2) as
601 *
602 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
603 * = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
604 *
605 * and store L(k) and L(k+1) in columns k and k+1
606 *
607  IF( k.GT.2 ) THEN
608 * D = |A12|
609  d = dlapy2( dble( a( k-1, k ) ),
610  $ dimag( a( k-1, k ) ) )
611  d11 = dble( a( k, k ) / d )
612  d22 = dble( a( k-1, k-1 ) / d )
613  d12 = a( k-1, k ) / d
614  tt = one / ( d11*d22-one )
615 *
616  DO 30 j = k - 2, 1, -1
617 *
618 * Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
619 *
620  wkm1 = tt*( d11*a( j, k-1 )-dconjg( d12 )*
621  $ a( j, k ) )
622  wk = tt*( d22*a( j, k )-d12*a( j, k-1 ) )
623 *
624 * Perform a rank-2 update of A(1:k-2,1:k-2)
625 *
626  DO 20 i = j, 1, -1
627  a( i, j ) = a( i, j ) -
628  $ ( a( i, k ) / d )*dconjg( wk ) -
629  $ ( a( i, k-1 ) / d )*dconjg( wkm1 )
630  20 CONTINUE
631 *
632 * Store U(k) and U(k-1) in cols k and k-1 for row J
633 *
634  a( j, k ) = wk / d
635  a( j, k-1 ) = wkm1 / d
636 * (*) Make sure that diagonal element of pivot is real
637  a( j, j ) = dcmplx( dble( a( j, j ) ), zero )
638 *
639  30 CONTINUE
640 *
641  END IF
642 *
643 * Copy superdiagonal elements of D(K) to E(K) and
644 * ZERO out superdiagonal entry of A
645 *
646  e( k ) = a( k-1, k )
647  e( k-1 ) = czero
648  a( k-1, k ) = czero
649 *
650  END IF
651 *
652 * End column K is nonsingular
653 *
654  END IF
655 *
656 * Store details of the interchanges in IPIV
657 *
658  IF( kstep.EQ.1 ) THEN
659  ipiv( k ) = kp
660  ELSE
661  ipiv( k ) = -p
662  ipiv( k-1 ) = -kp
663  END IF
664 *
665 * Decrease K and return to the start of the main loop
666 *
667  k = k - kstep
668  GO TO 10
669 *
670  34 CONTINUE
671 *
672  ELSE
673 *
674 * Factorize A as L*D*L**H using the lower triangle of A
675 *
676 * Initialize the unused last entry of the subdiagonal array E.
677 *
678  e( n ) = czero
679 *
680 * K is the main loop index, increasing from 1 to N in steps of
681 * 1 or 2
682 *
683  k = 1
684  40 CONTINUE
685 *
686 * If K > N, exit from loop
687 *
688  IF( k.GT.n )
689  $ GO TO 64
690  kstep = 1
691  p = k
692 *
693 * Determine rows and columns to be interchanged and whether
694 * a 1-by-1 or 2-by-2 pivot block will be used
695 *
696  absakk = abs( dble( a( k, k ) ) )
697 *
698 * IMAX is the row-index of the largest off-diagonal element in
699 * column K, and COLMAX is its absolute value.
700 * Determine both COLMAX and IMAX.
701 *
702  IF( k.LT.n ) THEN
703  imax = k + izamax( n-k, a( k+1, k ), 1 )
704  colmax = cabs1( a( imax, k ) )
705  ELSE
706  colmax = zero
707  END IF
708 *
709  IF( max( absakk, colmax ).EQ.zero ) THEN
710 *
711 * Column K is zero or underflow: set INFO and continue
712 *
713  IF( info.EQ.0 )
714  $ info = k
715  kp = k
716  a( k, k ) = dble( a( k, k ) )
717 *
718 * Set E( K ) to zero
719 *
720  IF( k.LT.n )
721  $ e( k ) = czero
722 *
723  ELSE
724 *
725 * ============================================================
726 *
727 * BEGIN pivot search
728 *
729 * Case(1)
730 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
731 * (used to handle NaN and Inf)
732 *
733  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
734 *
735 * no interchange, use 1-by-1 pivot block
736 *
737  kp = k
738 *
739  ELSE
740 *
741  done = .false.
742 *
743 * Loop until pivot found
744 *
745  42 CONTINUE
746 *
747 * BEGIN pivot search loop body
748 *
749 *
750 * JMAX is the column-index of the largest off-diagonal
751 * element in row IMAX, and ROWMAX is its absolute value.
752 * Determine both ROWMAX and JMAX.
753 *
754  IF( imax.NE.k ) THEN
755  jmax = k - 1 + izamax( imax-k, a( imax, k ), lda )
756  rowmax = cabs1( a( imax, jmax ) )
757  ELSE
758  rowmax = zero
759  END IF
760 *
761  IF( imax.LT.n ) THEN
762  itemp = imax + izamax( n-imax, a( imax+1, imax ),
763  $ 1 )
764  dtemp = cabs1( a( itemp, imax ) )
765  IF( dtemp.GT.rowmax ) THEN
766  rowmax = dtemp
767  jmax = itemp
768  END IF
769  END IF
770 *
771 * Case(2)
772 * Equivalent to testing for
773 * ABS( DBLE( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
774 * (used to handle NaN and Inf)
775 *
776  IF( .NOT.( abs( dble( a( imax, imax ) ) )
777  $ .LT.alpha*rowmax ) ) THEN
778 *
779 * interchange rows and columns K and IMAX,
780 * use 1-by-1 pivot block
781 *
782  kp = imax
783  done = .true.
784 *
785 * Case(3)
786 * Equivalent to testing for ROWMAX.EQ.COLMAX,
787 * (used to handle NaN and Inf)
788 *
789  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
790  $ THEN
791 *
792 * interchange rows and columns K+1 and IMAX,
793 * use 2-by-2 pivot block
794 *
795  kp = imax
796  kstep = 2
797  done = .true.
798 *
799 * Case(4)
800  ELSE
801 *
802 * Pivot not found: set params and repeat
803 *
804  p = imax
805  colmax = rowmax
806  imax = jmax
807  END IF
808 *
809 *
810 * END pivot search loop body
811 *
812  IF( .NOT.done ) GOTO 42
813 *
814  END IF
815 *
816 * END pivot search
817 *
818 * ============================================================
819 *
820 * KK is the column of A where pivoting step stopped
821 *
822  kk = k + kstep - 1
823 *
824 * For only a 2x2 pivot, interchange rows and columns K and P
825 * in the trailing submatrix A(k:n,k:n)
826 *
827  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
828 * (1) Swap columnar parts
829  IF( p.LT.n )
830  $ CALL zswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
831 * (2) Swap and conjugate middle parts
832  DO 44 j = k + 1, p - 1
833  t = dconjg( a( j, k ) )
834  a( j, k ) = dconjg( a( p, j ) )
835  a( p, j ) = t
836  44 CONTINUE
837 * (3) Swap and conjugate corner elements at row-col interserction
838  a( p, k ) = dconjg( a( p, k ) )
839 * (4) Swap diagonal elements at row-col intersection
840  r1 = dble( a( k, k ) )
841  a( k, k ) = dble( a( p, p ) )
842  a( p, p ) = r1
843 *
844 * Convert lower triangle of A into L form by applying
845 * the interchanges in columns 1:k-1.
846 *
847  IF ( k.GT.1 )
848  $ CALL zswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
849 *
850  END IF
851 *
852 * For both 1x1 and 2x2 pivots, interchange rows and
853 * columns KK and KP in the trailing submatrix A(k:n,k:n)
854 *
855  IF( kp.NE.kk ) THEN
856 * (1) Swap columnar parts
857  IF( kp.LT.n )
858  $ CALL zswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
859 * (2) Swap and conjugate middle parts
860  DO 45 j = kk + 1, kp - 1
861  t = dconjg( a( j, kk ) )
862  a( j, kk ) = dconjg( a( kp, j ) )
863  a( kp, j ) = t
864  45 CONTINUE
865 * (3) Swap and conjugate corner elements at row-col interserction
866  a( kp, kk ) = dconjg( a( kp, kk ) )
867 * (4) Swap diagonal elements at row-col intersection
868  r1 = dble( a( kk, kk ) )
869  a( kk, kk ) = dble( a( kp, kp ) )
870  a( kp, kp ) = r1
871 *
872  IF( kstep.EQ.2 ) THEN
873 * (*) Make sure that diagonal element of pivot is real
874  a( k, k ) = dble( a( k, k ) )
875 * (5) Swap row elements
876  t = a( k+1, k )
877  a( k+1, k ) = a( kp, k )
878  a( kp, k ) = t
879  END IF
880 *
881 * Convert lower triangle of A into L form by applying
882 * the interchanges in columns 1:k-1.
883 *
884  IF ( k.GT.1 )
885  $ CALL zswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
886 *
887  ELSE
888 * (*) Make sure that diagonal element of pivot is real
889  a( k, k ) = dble( a( k, k ) )
890  IF( kstep.EQ.2 )
891  $ a( k+1, k+1 ) = dble( a( k+1, k+1 ) )
892  END IF
893 *
894 * Update the trailing submatrix
895 *
896  IF( kstep.EQ.1 ) THEN
897 *
898 * 1-by-1 pivot block D(k): column k of A now holds
899 *
900 * W(k) = L(k)*D(k),
901 *
902 * where L(k) is the k-th column of L
903 *
904  IF( k.LT.n ) THEN
905 *
906 * Perform a rank-1 update of A(k+1:n,k+1:n) and
907 * store L(k) in column k
908 *
909 * Handle division by a small number
910 *
911  IF( abs( dble( a( k, k ) ) ).GE.sfmin ) THEN
912 *
913 * Perform a rank-1 update of A(k+1:n,k+1:n) as
914 * A := A - L(k)*D(k)*L(k)**T
915 * = A - W(k)*(1/D(k))*W(k)**T
916 *
917  d11 = one / dble( a( k, k ) )
918  CALL zher( uplo, n-k, -d11, a( k+1, k ), 1,
919  $ a( k+1, k+1 ), lda )
920 *
921 * Store L(k) in column k
922 *
923  CALL zdscal( n-k, d11, a( k+1, k ), 1 )
924  ELSE
925 *
926 * Store L(k) in column k
927 *
928  d11 = dble( a( k, k ) )
929  DO 46 ii = k + 1, n
930  a( ii, k ) = a( ii, k ) / d11
931  46 CONTINUE
932 *
933 * Perform a rank-1 update of A(k+1:n,k+1:n) as
934 * A := A - L(k)*D(k)*L(k)**T
935 * = A - W(k)*(1/D(k))*W(k)**T
936 * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
937 *
938  CALL zher( uplo, n-k, -d11, a( k+1, k ), 1,
939  $ a( k+1, k+1 ), lda )
940  END IF
941 *
942 * Store the subdiagonal element of D in array E
943 *
944  e( k ) = czero
945 *
946  END IF
947 *
948  ELSE
949 *
950 * 2-by-2 pivot block D(k): columns k and k+1 now hold
951 *
952 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
953 *
954 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
955 * of L
956 *
957 *
958 * Perform a rank-2 update of A(k+2:n,k+2:n) as
959 *
960 * A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
961 * = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
962 *
963 * and store L(k) and L(k+1) in columns k and k+1
964 *
965  IF( k.LT.n-1 ) THEN
966 * D = |A21|
967  d = dlapy2( dble( a( k+1, k ) ),
968  $ dimag( a( k+1, k ) ) )
969  d11 = dble( a( k+1, k+1 ) ) / d
970  d22 = dble( a( k, k ) ) / d
971  d21 = a( k+1, k ) / d
972  tt = one / ( d11*d22-one )
973 *
974  DO 60 j = k + 2, n
975 *
976 * Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
977 *
978  wk = tt*( d11*a( j, k )-d21*a( j, k+1 ) )
979  wkp1 = tt*( d22*a( j, k+1 )-dconjg( d21 )*
980  $ a( j, k ) )
981 *
982 * Perform a rank-2 update of A(k+2:n,k+2:n)
983 *
984  DO 50 i = j, n
985  a( i, j ) = a( i, j ) -
986  $ ( a( i, k ) / d )*dconjg( wk ) -
987  $ ( a( i, k+1 ) / d )*dconjg( wkp1 )
988  50 CONTINUE
989 *
990 * Store L(k) and L(k+1) in cols k and k+1 for row J
991 *
992  a( j, k ) = wk / d
993  a( j, k+1 ) = wkp1 / d
994 * (*) Make sure that diagonal element of pivot is real
995  a( j, j ) = dcmplx( dble( a( j, j ) ), zero )
996 *
997  60 CONTINUE
998 *
999  END IF
1000 *
1001 * Copy subdiagonal elements of D(K) to E(K) and
1002 * ZERO out subdiagonal entry of A
1003 *
1004  e( k ) = a( k+1, k )
1005  e( k+1 ) = czero
1006  a( k+1, k ) = czero
1007 *
1008  END IF
1009 *
1010 * End column K is nonsingular
1011 *
1012  END IF
1013 *
1014 * Store details of the interchanges in IPIV
1015 *
1016  IF( kstep.EQ.1 ) THEN
1017  ipiv( k ) = kp
1018  ELSE
1019  ipiv( k ) = -p
1020  ipiv( k+1 ) = -kp
1021  END IF
1022 *
1023 * Increase K and return to the start of the main loop
1024 *
1025  k = k + kstep
1026  GO TO 40
1027 *
1028  64 CONTINUE
1029 *
1030  END IF
1031 *
1032  RETURN
1033 *
1034 * End of ZHETF2_RK
1035 *
1036  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:78
subroutine zher(UPLO, N, ALPHA, X, INCX, A, LDA)
ZHER
Definition: zher.f:135
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