LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
csytf2_rook.f
Go to the documentation of this file.
1 *> \brief \b CSYTF2_ROOK computes the factorization of a complex symmetric 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
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytf2_rook.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytf2_rook.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytf2_rook.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CSYTF2_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 *> CSYTF2_ROOK computes the factorization of a complex symmetric matrix A
39 *> using the bounded Bunch-Kaufman ("rook") diagonal pivoting method:
40 *>
41 *> A = U*D*U**T or A = L*D*L**T
42 *>
43 *> where U (or L) is a product of permutation and unit upper (lower)
44 *> triangular matrices, U**T is the transpose of U, and D is symmetric and
45 *> 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 *> symmetric 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 symmetric 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)
96 *> were 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 complexSYcomputational
135 *
136 *> \par Further Details:
137 * =====================
138 *>
139 *> \verbatim
140 *>
141 *> If UPLO = 'U', then A = U*D*U**T, 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**T, 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 abd , USA
192 *> \endverbatim
193 *
194 * =====================================================================
195  SUBROUTINE csytf2_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  COMPLEX CONE
219  parameter ( cone = ( 1.0e+0, 0.0e+0 ) )
220 * ..
221 * .. Local Scalars ..
222  LOGICAL UPPER, DONE
223  INTEGER I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP,
224  \$ p, ii
225  REAL ABSAKK, ALPHA, COLMAX, ROWMAX, STEMP, SFMIN
226  COMPLEX D11, D12, D21, D22, T, WK, WKM1, WKP1, Z
227 * ..
228 * .. External Functions ..
229  LOGICAL LSAME
230  INTEGER ICAMAX
231  REAL SLAMCH
232  EXTERNAL lsame, icamax, slamch
233 * ..
234 * .. External Subroutines ..
235  EXTERNAL cscal, cswap, csyr, xerbla
236 * ..
237 * .. Intrinsic Functions ..
238  INTRINSIC abs, max, sqrt, aimag, real
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( 'CSYTF2_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**T 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 = cabs1( 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  ELSE
313 *
314 * Test for interchange
315 *
316 * Equivalent to testing for (used to handle NaN and Inf)
317 * ABSAKK.GE.ALPHA*COLMAX
318 *
319  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
320 *
321 * no interchange,
322 * use 1-by-1 pivot block
323 *
324  kp = k
325  ELSE
326 *
327  done = .false.
328 *
329 * Loop until pivot found
330 *
331  12 CONTINUE
332 *
333 * Begin pivot search loop body
334 *
335 * JMAX is the column-index of the largest off-diagonal
336 * element in row IMAX, and ROWMAX is its absolute value.
337 * Determine both ROWMAX and JMAX.
338 *
339  IF( imax.NE.k ) THEN
340  jmax = imax + icamax( k-imax, a( imax, imax+1 ),
341  \$ lda )
342  rowmax = cabs1( a( imax, jmax ) )
343  ELSE
344  rowmax = zero
345  END IF
346 *
347  IF( imax.GT.1 ) THEN
348  itemp = icamax( imax-1, a( 1, imax ), 1 )
349  stemp = cabs1( a( itemp, imax ) )
350  IF( stemp.GT.rowmax ) THEN
351  rowmax = stemp
352  jmax = itemp
353  END IF
354  END IF
355 *
356 * Equivalent to testing for (used to handle NaN and Inf)
357 * CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
358 *
359  IF( .NOT.( cabs1(a( imax, imax )).LT.alpha*rowmax ) )
360  \$ THEN
361 *
362 * interchange rows and columns K and IMAX,
363 * use 1-by-1 pivot block
364 *
365  kp = imax
366  done = .true.
367 *
368 * Equivalent to testing for ROWMAX .EQ. COLMAX,
369 * used to handle NaN and Inf
370 *
371  ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
372 *
373 * interchange rows and columns K+1 and IMAX,
374 * use 2-by-2 pivot block
375 *
376  kp = imax
377  kstep = 2
378  done = .true.
379  ELSE
380 *
382 *
383  p = imax
384  colmax = rowmax
385  imax = jmax
386  END IF
387 *
388 * End pivot search loop body
389 *
390  IF( .NOT. done ) GOTO 12
391 *
392  END IF
393 *
394 * Swap TWO rows and TWO columns
395 *
396 * First swap
397 *
398  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
399 *
400 * Interchange rows and column K and P in the leading
401 * submatrix A(1:k,1:k) if we have a 2-by-2 pivot
402 *
403  IF( p.GT.1 )
404  \$ CALL cswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
405  IF( p.LT.(k-1) )
406  \$ CALL cswap( k-p-1, a( p+1, k ), 1, a( p, p+1 ),
407  \$ lda )
408  t = a( k, k )
409  a( k, k ) = a( p, p )
410  a( p, p ) = t
411  END IF
412 *
413 * Second swap
414 *
415  kk = k - kstep + 1
416  IF( kp.NE.kk ) THEN
417 *
418 * Interchange rows and columns KK and KP in the leading
419 * submatrix A(1:k,1:k)
420 *
421  IF( kp.GT.1 )
422  \$ CALL cswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
423  IF( ( kk.GT.1 ) .AND. ( kp.LT.(kk-1) ) )
424  \$ CALL cswap( kk-kp-1, a( kp+1, kk ), 1, a( kp, kp+1 ),
425  \$ lda )
426  t = a( kk, kk )
427  a( kk, kk ) = a( kp, kp )
428  a( kp, kp ) = t
429  IF( kstep.EQ.2 ) THEN
430  t = a( k-1, k )
431  a( k-1, k ) = a( kp, k )
432  a( kp, k ) = t
433  END IF
434  END IF
435 *
436 * Update the leading submatrix
437 *
438  IF( kstep.EQ.1 ) THEN
439 *
440 * 1-by-1 pivot block D(k): column k now holds
441 *
442 * W(k) = U(k)*D(k)
443 *
444 * where U(k) is the k-th column of U
445 *
446  IF( k.GT.1 ) THEN
447 *
448 * Perform a rank-1 update of A(1:k-1,1:k-1) and
449 * store U(k) in column k
450 *
451  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
452 *
453 * Perform a rank-1 update of A(1:k-1,1:k-1) as
454 * A := A - U(k)*D(k)*U(k)**T
455 * = A - W(k)*1/D(k)*W(k)**T
456 *
457  d11 = cone / a( k, k )
458  CALL csyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
459 *
460 * Store U(k) in column k
461 *
462  CALL cscal( k-1, d11, a( 1, k ), 1 )
463  ELSE
464 *
465 * Store L(k) in column K
466 *
467  d11 = a( k, k )
468  DO 16 ii = 1, k - 1
469  a( ii, k ) = a( ii, k ) / d11
470  16 CONTINUE
471 *
472 * Perform a rank-1 update of A(k+1:n,k+1:n) as
473 * A := A - U(k)*D(k)*U(k)**T
474 * = A - W(k)*(1/D(k))*W(k)**T
475 * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
476 *
477  CALL csyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
478  END IF
479  END IF
480 *
481  ELSE
482 *
483 * 2-by-2 pivot block D(k): columns k and k-1 now hold
484 *
485 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
486 *
487 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
488 * of U
489 *
490 * Perform a rank-2 update of A(1:k-2,1:k-2) as
491 *
492 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
493 * = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
494 *
495 * and store L(k) and L(k+1) in columns k and k+1
496 *
497  IF( k.GT.2 ) THEN
498 *
499  d12 = a( k-1, k )
500  d22 = a( k-1, k-1 ) / d12
501  d11 = a( k, k ) / d12
502  t = cone / ( d11*d22-cone )
503 *
504  DO 30 j = k - 2, 1, -1
505 *
506  wkm1 = t*( d11*a( j, k-1 )-a( j, k ) )
507  wk = t*( d22*a( j, k )-a( j, k-1 ) )
508 *
509  DO 20 i = j, 1, -1
510  a( i, j ) = a( i, j ) - (a( i, k ) / d12 )*wk -
511  \$ ( a( i, k-1 ) / d12 )*wkm1
512  20 CONTINUE
513 *
514 * Store U(k) and U(k-1) in cols k and k-1 for row J
515 *
516  a( j, k ) = wk / d12
517  a( j, k-1 ) = wkm1 / d12
518 *
519  30 CONTINUE
520 *
521  END IF
522 *
523  END IF
524  END IF
525 *
526 * Store details of the interchanges in IPIV
527 *
528  IF( kstep.EQ.1 ) THEN
529  ipiv( k ) = kp
530  ELSE
531  ipiv( k ) = -p
532  ipiv( k-1 ) = -kp
533  END IF
534 *
535 * Decrease K and return to the start of the main loop
536 *
537  k = k - kstep
538  GO TO 10
539 *
540  ELSE
541 *
542 * Factorize A as L*D*L**T using the lower triangle of A
543 *
544 * K is the main loop index, increasing from 1 to N in steps of
545 * 1 or 2
546 *
547  k = 1
548  40 CONTINUE
549 *
550 * If K > N, exit from loop
551 *
552  IF( k.GT.n )
553  \$ GO TO 70
554  kstep = 1
555  p = k
556 *
557 * Determine rows and columns to be interchanged and whether
558 * a 1-by-1 or 2-by-2 pivot block will be used
559 *
560  absakk = cabs1( a( k, k ) )
561 *
562 * IMAX is the row-index of the largest off-diagonal element in
563 * column K, and COLMAX is its absolute value.
564 * Determine both COLMAX and IMAX.
565 *
566  IF( k.LT.n ) THEN
567  imax = k + icamax( n-k, a( k+1, k ), 1 )
568  colmax = cabs1( a( imax, k ) )
569  ELSE
570  colmax = zero
571  END IF
572 *
573  IF( ( max( absakk, colmax ).EQ.zero ) ) THEN
574 *
575 * Column K is zero or underflow: set INFO and continue
576 *
577  IF( info.EQ.0 )
578  \$ info = k
579  kp = k
580  ELSE
581 *
582 * Test for interchange
583 *
584 * Equivalent to testing for (used to handle NaN and Inf)
585 * ABSAKK.GE.ALPHA*COLMAX
586 *
587  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
588 *
589 * no interchange, use 1-by-1 pivot block
590 *
591  kp = k
592  ELSE
593 *
594  done = .false.
595 *
596 * Loop until pivot found
597 *
598  42 CONTINUE
599 *
600 * Begin pivot search loop body
601 *
602 * JMAX is the column-index of the largest off-diagonal
603 * element in row IMAX, and ROWMAX is its absolute value.
604 * Determine both ROWMAX and JMAX.
605 *
606  IF( imax.NE.k ) THEN
607  jmax = k - 1 + icamax( imax-k, a( imax, k ), lda )
608  rowmax = cabs1( a( imax, jmax ) )
609  ELSE
610  rowmax = zero
611  END IF
612 *
613  IF( imax.LT.n ) THEN
614  itemp = imax + icamax( n-imax, a( imax+1, imax ),
615  \$ 1 )
616  stemp = cabs1( a( itemp, imax ) )
617  IF( stemp.GT.rowmax ) THEN
618  rowmax = stemp
619  jmax = itemp
620  END IF
621  END IF
622 *
623 * Equivalent to testing for (used to handle NaN and Inf)
624 * CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
625 *
626  IF( .NOT.( cabs1(a( imax, imax )).LT.alpha*rowmax ) )
627  \$ THEN
628 *
629 * interchange rows and columns K and IMAX,
630 * use 1-by-1 pivot block
631 *
632  kp = imax
633  done = .true.
634 *
635 * Equivalent to testing for ROWMAX .EQ. COLMAX,
636 * used to handle NaN and Inf
637 *
638  ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
639 *
640 * interchange rows and columns K+1 and IMAX,
641 * use 2-by-2 pivot block
642 *
643  kp = imax
644  kstep = 2
645  done = .true.
646  ELSE
647 *
649 *
650  p = imax
651  colmax = rowmax
652  imax = jmax
653  END IF
654 *
655 * End pivot search loop body
656 *
657  IF( .NOT. done ) GOTO 42
658 *
659  END IF
660 *
661 * Swap TWO rows and TWO columns
662 *
663 * First swap
664 *
665  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
666 *
667 * Interchange rows and column K and P in the trailing
668 * submatrix A(k:n,k:n) if we have a 2-by-2 pivot
669 *
670  IF( p.LT.n )
671  \$ CALL cswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
672  IF( p.GT.(k+1) )
673  \$ CALL cswap( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
674  t = a( k, k )
675  a( k, k ) = a( p, p )
676  a( p, p ) = t
677  END IF
678 *
679 * Second swap
680 *
681  kk = k + kstep - 1
682  IF( kp.NE.kk ) THEN
683 *
684 * Interchange rows and columns KK and KP in the trailing
685 * submatrix A(k:n,k:n)
686 *
687  IF( kp.LT.n )
688  \$ CALL cswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
689  IF( ( kk.LT.n ) .AND. ( kp.GT.(kk+1) ) )
690  \$ CALL cswap( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
691  \$ lda )
692  t = a( kk, kk )
693  a( kk, kk ) = a( kp, kp )
694  a( kp, kp ) = t
695  IF( kstep.EQ.2 ) THEN
696  t = a( k+1, k )
697  a( k+1, k ) = a( kp, k )
698  a( kp, k ) = t
699  END IF
700  END IF
701 *
702 * Update the trailing submatrix
703 *
704  IF( kstep.EQ.1 ) THEN
705 *
706 * 1-by-1 pivot block D(k): column k now holds
707 *
708 * W(k) = L(k)*D(k)
709 *
710 * where L(k) is the k-th column of L
711 *
712  IF( k.LT.n ) THEN
713 *
714 * Perform a rank-1 update of A(k+1:n,k+1:n) and
715 * store L(k) in column k
716 *
717  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
718 *
719 * Perform a rank-1 update of A(k+1:n,k+1:n) as
720 * A := A - L(k)*D(k)*L(k)**T
721 * = A - W(k)*(1/D(k))*W(k)**T
722 *
723  d11 = cone / a( k, k )
724  CALL csyr( uplo, n-k, -d11, a( k+1, k ), 1,
725  \$ a( k+1, k+1 ), lda )
726 *
727 * Store L(k) in column k
728 *
729  CALL cscal( n-k, d11, a( k+1, k ), 1 )
730  ELSE
731 *
732 * Store L(k) in column k
733 *
734  d11 = a( k, k )
735  DO 46 ii = k + 1, n
736  a( ii, k ) = a( ii, k ) / d11
737  46 CONTINUE
738 *
739 * Perform a rank-1 update of A(k+1:n,k+1:n) as
740 * A := A - L(k)*D(k)*L(k)**T
741 * = A - W(k)*(1/D(k))*W(k)**T
742 * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
743 *
744  CALL csyr( uplo, n-k, -d11, a( k+1, k ), 1,
745  \$ a( k+1, k+1 ), lda )
746  END IF
747  END IF
748 *
749  ELSE
750 *
751 * 2-by-2 pivot block D(k): columns k and k+1 now hold
752 *
753 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
754 *
755 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
756 * of L
757 *
758 *
759 * Perform a rank-2 update of A(k+2:n,k+2:n) as
760 *
761 * A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
762 * = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
763 *
764 * and store L(k) and L(k+1) in columns k and k+1
765 *
766  IF( k.LT.n-1 ) THEN
767 *
768  d21 = a( k+1, k )
769  d11 = a( k+1, k+1 ) / d21
770  d22 = a( k, k ) / d21
771  t = cone / ( d11*d22-cone )
772 *
773  DO 60 j = k + 2, n
774 *
775 * Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
776 *
777  wk = t*( d11*a( j, k )-a( j, k+1 ) )
778  wkp1 = t*( d22*a( j, k+1 )-a( j, k ) )
779 *
780 * Perform a rank-2 update of A(k+2:n,k+2:n)
781 *
782  DO 50 i = j, n
783  a( i, j ) = a( i, j ) - ( a( i, k ) / d21 )*wk -
784  \$ ( a( i, k+1 ) / d21 )*wkp1
785  50 CONTINUE
786 *
787 * Store L(k) and L(k+1) in cols k and k+1 for row J
788 *
789  a( j, k ) = wk / d21
790  a( j, k+1 ) = wkp1 / d21
791 *
792  60 CONTINUE
793 *
794  END IF
795 *
796  END IF
797  END IF
798 *
799 * Store details of the interchanges in IPIV
800 *
801  IF( kstep.EQ.1 ) THEN
802  ipiv( k ) = kp
803  ELSE
804  ipiv( k ) = -p
805  ipiv( k+1 ) = -kp
806  END IF
807 *
808 * Increase K and return to the start of the main loop
809 *
810  k = k + kstep
811  GO TO 40
812 *
813  END IF
814 *
815  70 CONTINUE
816 *
817  RETURN
818 *
819 * End of CSYTF2_ROOK
820 *
821  END
subroutine csyr(UPLO, N, ALPHA, X, INCX, A, LDA)
CSYR performs the symmetric rank-1 update of a complex symmetric matrix.
Definition: csyr.f:137
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:54
subroutine csytf2_rook(UPLO, N, A, LDA, IPIV, INFO)
CSYTF2_ROOK computes the factorization of a complex symmetric indefinite matrix using the bounded Bun...
Definition: csytf2_rook.f:196
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52