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