LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zsytf2_rk.f
Go to the documentation of this file.
1*> \brief \b ZSYTF2_RK computes the factorization of a complex symmetric 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 ZSYTF2_RK + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytf2_rk.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytf2_rk.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytf2_rk.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZSYTF2_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*> ZSYTF2_RK computes the factorization of a complex symmetric matrix A
38*> using the bounded Bunch-Kaufman (rook) diagonal pivoting method:
39*>
40*> A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T),
41*>
42*> where U (or L) is unit upper (or lower) triangular matrix,
43*> U**T (or L**T) is the transpose of U (or L), P is a permutation
44*> matrix, P**T is the transpose of P, and D is symmetric 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*> symmetric 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 symmetric 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 symmetric 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 symmetric 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 symmetric 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 complex16SYcomputational
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 zsytf2_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 CONE, CZERO
263 parameter( cone = ( 1.0d+0, 0.0d+0 ),
264 $ czero = ( 0.0d+0, 0.0d+0 ) )
265* ..
266* .. Local Scalars ..
267 LOGICAL UPPER, DONE
268 INTEGER I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP,
269 $ P, II
270 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, ROWMAX, DTEMP, SFMIN
271 COMPLEX*16 D11, D12, D21, D22, T, WK, WKM1, WKP1, Z
272* ..
273* .. External Functions ..
274 LOGICAL LSAME
275 INTEGER IZAMAX
276 DOUBLE PRECISION DLAMCH
277 EXTERNAL lsame, izamax, dlamch
278* ..
279* .. External Subroutines ..
280 EXTERNAL zscal, zswap, zsyr, xerbla
281* ..
282* .. Intrinsic Functions ..
283 INTRINSIC abs, max, sqrt, dimag, dble
284* ..
285* .. Statement Functions ..
286 DOUBLE PRECISION CABS1
287* ..
288* .. Statement Function definitions ..
289 cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
290* ..
291* .. Executable Statements ..
292*
293* Test the input parameters.
294*
295 info = 0
296 upper = lsame( uplo, 'U' )
297 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
298 info = -1
299 ELSE IF( n.LT.0 ) THEN
300 info = -2
301 ELSE IF( lda.LT.max( 1, n ) ) THEN
302 info = -4
303 END IF
304 IF( info.NE.0 ) THEN
305 CALL xerbla( 'ZSYTF2_RK', -info )
306 RETURN
307 END IF
308*
309* Initialize ALPHA for use in choosing pivot block size.
310*
311 alpha = ( one+sqrt( sevten ) ) / eight
312*
313* Compute machine safe minimum
314*
315 sfmin = dlamch( 'S' )
316*
317 IF( upper ) THEN
318*
319* Factorize A as U*D*U**T using the upper triangle of A
320*
321* Initialize the first entry of array E, where superdiagonal
322* elements of D are stored
323*
324 e( 1 ) = czero
325*
326* K is the main loop index, decreasing from N to 1 in steps of
327* 1 or 2
328*
329 k = n
330 10 CONTINUE
331*
332* If K < 1, exit from loop
333*
334 IF( k.LT.1 )
335 $ GO TO 34
336 kstep = 1
337 p = k
338*
339* Determine rows and columns to be interchanged and whether
340* a 1-by-1 or 2-by-2 pivot block will be used
341*
342 absakk = cabs1( a( k, k ) )
343*
344* IMAX is the row-index of the largest off-diagonal element in
345* column K, and COLMAX is its absolute value.
346* Determine both COLMAX and IMAX.
347*
348 IF( k.GT.1 ) THEN
349 imax = izamax( k-1, a( 1, k ), 1 )
350 colmax = cabs1( a( imax, k ) )
351 ELSE
352 colmax = zero
353 END IF
354*
355 IF( (max( absakk, colmax ).EQ.zero) ) THEN
356*
357* Column K is zero or underflow: set INFO and continue
358*
359 IF( info.EQ.0 )
360 $ info = k
361 kp = k
362*
363* Set E( K ) to zero
364*
365 IF( k.GT.1 )
366 $ e( k ) = czero
367*
368 ELSE
369*
370* Test for interchange
371*
372* Equivalent to testing for (used to handle NaN and Inf)
373* ABSAKK.GE.ALPHA*COLMAX
374*
375 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
376*
377* no interchange,
378* use 1-by-1 pivot block
379*
380 kp = k
381 ELSE
382*
383 done = .false.
384*
385* Loop until pivot found
386*
387 12 CONTINUE
388*
389* Begin pivot search loop body
390*
391* JMAX is the column-index of the largest off-diagonal
392* element in row IMAX, and ROWMAX is its absolute value.
393* Determine both ROWMAX and JMAX.
394*
395 IF( imax.NE.k ) THEN
396 jmax = imax + izamax( k-imax, a( imax, imax+1 ),
397 $ lda )
398 rowmax = cabs1( a( imax, jmax ) )
399 ELSE
400 rowmax = zero
401 END IF
402*
403 IF( imax.GT.1 ) THEN
404 itemp = izamax( imax-1, a( 1, imax ), 1 )
405 dtemp = cabs1( a( itemp, imax ) )
406 IF( dtemp.GT.rowmax ) THEN
407 rowmax = dtemp
408 jmax = itemp
409 END IF
410 END IF
411*
412* Equivalent to testing for (used to handle NaN and Inf)
413* ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
414*
415 IF( .NOT.( cabs1( a( imax, imax ) ).LT.alpha*rowmax ))
416 $ THEN
417*
418* interchange rows and columns K and IMAX,
419* use 1-by-1 pivot block
420*
421 kp = imax
422 done = .true.
423*
424* Equivalent to testing for ROWMAX .EQ. COLMAX,
425* used to handle NaN and Inf
426*
427 ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
428*
429* interchange rows and columns K+1 and IMAX,
430* use 2-by-2 pivot block
431*
432 kp = imax
433 kstep = 2
434 done = .true.
435 ELSE
436*
437* Pivot NOT found, set variables and repeat
438*
439 p = imax
440 colmax = rowmax
441 imax = jmax
442 END IF
443*
444* End pivot search loop body
445*
446 IF( .NOT. done ) GOTO 12
447*
448 END IF
449*
450* Swap TWO rows and TWO columns
451*
452* First swap
453*
454 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
455*
456* Interchange rows and column K and P in the leading
457* submatrix A(1:k,1:k) if we have a 2-by-2 pivot
458*
459 IF( p.GT.1 )
460 $ CALL zswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
461 IF( p.LT.(k-1) )
462 $ CALL zswap( k-p-1, a( p+1, k ), 1, a( p, p+1 ),
463 $ lda )
464 t = a( k, k )
465 a( k, k ) = a( p, p )
466 a( p, p ) = t
467*
468* Convert upper triangle of A into U form by applying
469* the interchanges in columns k+1:N.
470*
471 IF( k.LT.n )
472 $ CALL zswap( n-k, a( k, k+1 ), lda, a( p, k+1 ), lda )
473*
474 END IF
475*
476* Second swap
477*
478 kk = k - kstep + 1
479 IF( kp.NE.kk ) THEN
480*
481* Interchange rows and columns KK and KP in the leading
482* submatrix A(1:k,1:k)
483*
484 IF( kp.GT.1 )
485 $ CALL zswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
486 IF( ( kk.GT.1 ) .AND. ( kp.LT.(kk-1) ) )
487 $ CALL zswap( kk-kp-1, a( kp+1, kk ), 1, a( kp, kp+1 ),
488 $ lda )
489 t = a( kk, kk )
490 a( kk, kk ) = a( kp, kp )
491 a( kp, kp ) = t
492 IF( kstep.EQ.2 ) THEN
493 t = a( k-1, k )
494 a( k-1, k ) = a( kp, k )
495 a( kp, k ) = t
496 END IF
497*
498* Convert upper triangle of A into U form by applying
499* the interchanges in columns k+1:N.
500*
501 IF( k.LT.n )
502 $ CALL zswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
503 $ lda )
504*
505 END IF
506*
507* Update the leading submatrix
508*
509 IF( kstep.EQ.1 ) THEN
510*
511* 1-by-1 pivot block D(k): column k now holds
512*
513* W(k) = U(k)*D(k)
514*
515* where U(k) is the k-th column of U
516*
517 IF( k.GT.1 ) THEN
518*
519* Perform a rank-1 update of A(1:k-1,1:k-1) and
520* store U(k) in column k
521*
522 IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
523*
524* Perform a rank-1 update of A(1:k-1,1:k-1) as
525* A := A - U(k)*D(k)*U(k)**T
526* = A - W(k)*1/D(k)*W(k)**T
527*
528 d11 = cone / a( k, k )
529 CALL zsyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
530*
531* Store U(k) in column k
532*
533 CALL zscal( k-1, d11, a( 1, k ), 1 )
534 ELSE
535*
536* Store L(k) in column K
537*
538 d11 = a( k, k )
539 DO 16 ii = 1, k - 1
540 a( ii, k ) = a( ii, k ) / d11
541 16 CONTINUE
542*
543* Perform a rank-1 update of A(k+1:n,k+1:n) as
544* A := A - U(k)*D(k)*U(k)**T
545* = A - W(k)*(1/D(k))*W(k)**T
546* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
547*
548 CALL zsyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
549 END IF
550*
551* Store the superdiagonal element of D in array E
552*
553 e( k ) = czero
554*
555 END IF
556*
557 ELSE
558*
559* 2-by-2 pivot block D(k): columns k and k-1 now hold
560*
561* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
562*
563* where U(k) and U(k-1) are the k-th and (k-1)-th columns
564* of U
565*
566* Perform a rank-2 update of A(1:k-2,1:k-2) as
567*
568* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
569* = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
570*
571* and store L(k) and L(k+1) in columns k and k+1
572*
573 IF( k.GT.2 ) THEN
574*
575 d12 = a( k-1, k )
576 d22 = a( k-1, k-1 ) / d12
577 d11 = a( k, k ) / d12
578 t = cone / ( d11*d22-cone )
579*
580 DO 30 j = k - 2, 1, -1
581*
582 wkm1 = t*( d11*a( j, k-1 )-a( j, k ) )
583 wk = t*( d22*a( j, k )-a( j, k-1 ) )
584*
585 DO 20 i = j, 1, -1
586 a( i, j ) = a( i, j ) - (a( i, k ) / d12 )*wk -
587 $ ( a( i, k-1 ) / d12 )*wkm1
588 20 CONTINUE
589*
590* Store U(k) and U(k-1) in cols k and k-1 for row J
591*
592 a( j, k ) = wk / d12
593 a( j, k-1 ) = wkm1 / d12
594*
595 30 CONTINUE
596*
597 END IF
598*
599* Copy superdiagonal elements of D(K) to E(K) and
600* ZERO out superdiagonal entry of A
601*
602 e( k ) = a( k-1, k )
603 e( k-1 ) = czero
604 a( k-1, k ) = czero
605*
606 END IF
607*
608* End column K is nonsingular
609*
610 END IF
611*
612* Store details of the interchanges in IPIV
613*
614 IF( kstep.EQ.1 ) THEN
615 ipiv( k ) = kp
616 ELSE
617 ipiv( k ) = -p
618 ipiv( k-1 ) = -kp
619 END IF
620*
621* Decrease K and return to the start of the main loop
622*
623 k = k - kstep
624 GO TO 10
625*
626 34 CONTINUE
627*
628 ELSE
629*
630* Factorize A as L*D*L**T using the lower triangle of A
631*
632* Initialize the unused last entry of the subdiagonal array E.
633*
634 e( n ) = czero
635*
636* K is the main loop index, increasing from 1 to N in steps of
637* 1 or 2
638*
639 k = 1
640 40 CONTINUE
641*
642* If K > N, exit from loop
643*
644 IF( k.GT.n )
645 $ GO TO 64
646 kstep = 1
647 p = k
648*
649* Determine rows and columns to be interchanged and whether
650* a 1-by-1 or 2-by-2 pivot block will be used
651*
652 absakk = cabs1( a( k, k ) )
653*
654* IMAX is the row-index of the largest off-diagonal element in
655* column K, and COLMAX is its absolute value.
656* Determine both COLMAX and IMAX.
657*
658 IF( k.LT.n ) THEN
659 imax = k + izamax( n-k, a( k+1, k ), 1 )
660 colmax = cabs1( a( imax, k ) )
661 ELSE
662 colmax = zero
663 END IF
664*
665 IF( ( max( absakk, colmax ).EQ.zero ) ) THEN
666*
667* Column K is zero or underflow: set INFO and continue
668*
669 IF( info.EQ.0 )
670 $ info = k
671 kp = k
672*
673* Set E( K ) to zero
674*
675 IF( k.LT.n )
676 $ e( k ) = czero
677*
678 ELSE
679*
680* Test for interchange
681*
682* Equivalent to testing for (used to handle NaN and Inf)
683* ABSAKK.GE.ALPHA*COLMAX
684*
685 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
686*
687* no interchange, use 1-by-1 pivot block
688*
689 kp = k
690*
691 ELSE
692*
693 done = .false.
694*
695* Loop until pivot found
696*
697 42 CONTINUE
698*
699* Begin pivot search loop body
700*
701* JMAX is the column-index of the largest off-diagonal
702* element in row IMAX, and ROWMAX is its absolute value.
703* Determine both ROWMAX and JMAX.
704*
705 IF( imax.NE.k ) THEN
706 jmax = k - 1 + izamax( imax-k, a( imax, k ), lda )
707 rowmax = cabs1( a( imax, jmax ) )
708 ELSE
709 rowmax = zero
710 END IF
711*
712 IF( imax.LT.n ) THEN
713 itemp = imax + izamax( n-imax, a( imax+1, imax ),
714 $ 1 )
715 dtemp = cabs1( a( itemp, imax ) )
716 IF( dtemp.GT.rowmax ) THEN
717 rowmax = dtemp
718 jmax = itemp
719 END IF
720 END IF
721*
722* Equivalent to testing for (used to handle NaN and Inf)
723* ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
724*
725 IF( .NOT.( cabs1( a( imax, imax ) ).LT.alpha*rowmax ))
726 $ THEN
727*
728* interchange rows and columns K and IMAX,
729* use 1-by-1 pivot block
730*
731 kp = imax
732 done = .true.
733*
734* Equivalent to testing for ROWMAX .EQ. COLMAX,
735* used to handle NaN and Inf
736*
737 ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
738*
739* interchange rows and columns K+1 and IMAX,
740* use 2-by-2 pivot block
741*
742 kp = imax
743 kstep = 2
744 done = .true.
745 ELSE
746*
747* Pivot NOT found, set variables and repeat
748*
749 p = imax
750 colmax = rowmax
751 imax = jmax
752 END IF
753*
754* End pivot search loop body
755*
756 IF( .NOT. done ) GOTO 42
757*
758 END IF
759*
760* Swap TWO rows and TWO columns
761*
762* First swap
763*
764 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
765*
766* Interchange rows and column K and P in the trailing
767* submatrix A(k:n,k:n) if we have a 2-by-2 pivot
768*
769 IF( p.LT.n )
770 $ CALL zswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
771 IF( p.GT.(k+1) )
772 $ CALL zswap( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
773 t = a( k, k )
774 a( k, k ) = a( p, p )
775 a( p, p ) = t
776*
777* Convert lower triangle of A into L form by applying
778* the interchanges in columns 1:k-1.
779*
780 IF ( k.GT.1 )
781 $ CALL zswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
782*
783 END IF
784*
785* Second swap
786*
787 kk = k + kstep - 1
788 IF( kp.NE.kk ) THEN
789*
790* Interchange rows and columns KK and KP in the trailing
791* submatrix A(k:n,k:n)
792*
793 IF( kp.LT.n )
794 $ CALL zswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
795 IF( ( kk.LT.n ) .AND. ( kp.GT.(kk+1) ) )
796 $ CALL zswap( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
797 $ lda )
798 t = a( kk, kk )
799 a( kk, kk ) = a( kp, kp )
800 a( kp, kp ) = t
801 IF( kstep.EQ.2 ) THEN
802 t = a( k+1, k )
803 a( k+1, k ) = a( kp, k )
804 a( kp, k ) = t
805 END IF
806*
807* Convert lower triangle of A into L form by applying
808* the interchanges in columns 1:k-1.
809*
810 IF ( k.GT.1 )
811 $ CALL zswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
812*
813 END IF
814*
815* Update the trailing submatrix
816*
817 IF( kstep.EQ.1 ) THEN
818*
819* 1-by-1 pivot block D(k): column k now holds
820*
821* W(k) = L(k)*D(k)
822*
823* where L(k) is the k-th column of L
824*
825 IF( k.LT.n ) THEN
826*
827* Perform a rank-1 update of A(k+1:n,k+1:n) and
828* store L(k) in column k
829*
830 IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
831*
832* Perform a rank-1 update of A(k+1:n,k+1:n) as
833* A := A - L(k)*D(k)*L(k)**T
834* = A - W(k)*(1/D(k))*W(k)**T
835*
836 d11 = cone / a( k, k )
837 CALL zsyr( uplo, n-k, -d11, a( k+1, k ), 1,
838 $ a( k+1, k+1 ), lda )
839*
840* Store L(k) in column k
841*
842 CALL zscal( n-k, d11, a( k+1, k ), 1 )
843 ELSE
844*
845* Store L(k) in column k
846*
847 d11 = a( k, k )
848 DO 46 ii = k + 1, n
849 a( ii, k ) = a( ii, k ) / d11
850 46 CONTINUE
851*
852* Perform a rank-1 update of A(k+1:n,k+1:n) as
853* A := A - L(k)*D(k)*L(k)**T
854* = A - W(k)*(1/D(k))*W(k)**T
855* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
856*
857 CALL zsyr( uplo, n-k, -d11, a( k+1, k ), 1,
858 $ a( k+1, k+1 ), lda )
859 END IF
860*
861* Store the subdiagonal element of D in array E
862*
863 e( k ) = czero
864*
865 END IF
866*
867 ELSE
868*
869* 2-by-2 pivot block D(k): columns k and k+1 now hold
870*
871* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
872*
873* where L(k) and L(k+1) are the k-th and (k+1)-th columns
874* of L
875*
876*
877* Perform a rank-2 update of A(k+2:n,k+2:n) as
878*
879* A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
880* = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
881*
882* and store L(k) and L(k+1) in columns k and k+1
883*
884 IF( k.LT.n-1 ) THEN
885*
886 d21 = a( k+1, k )
887 d11 = a( k+1, k+1 ) / d21
888 d22 = a( k, k ) / d21
889 t = cone / ( d11*d22-cone )
890*
891 DO 60 j = k + 2, n
892*
893* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
894*
895 wk = t*( d11*a( j, k )-a( j, k+1 ) )
896 wkp1 = t*( d22*a( j, k+1 )-a( j, k ) )
897*
898* Perform a rank-2 update of A(k+2:n,k+2:n)
899*
900 DO 50 i = j, n
901 a( i, j ) = a( i, j ) - ( a( i, k ) / d21 )*wk -
902 $ ( a( i, k+1 ) / d21 )*wkp1
903 50 CONTINUE
904*
905* Store L(k) and L(k+1) in cols k and k+1 for row J
906*
907 a( j, k ) = wk / d21
908 a( j, k+1 ) = wkp1 / d21
909*
910 60 CONTINUE
911*
912 END IF
913*
914* Copy subdiagonal elements of D(K) to E(K) and
915* ZERO out subdiagonal entry of A
916*
917 e( k ) = a( k+1, k )
918 e( k+1 ) = czero
919 a( k+1, k ) = czero
920*
921 END IF
922*
923* End column K is nonsingular
924*
925 END IF
926*
927* Store details of the interchanges in IPIV
928*
929 IF( kstep.EQ.1 ) THEN
930 ipiv( k ) = kp
931 ELSE
932 ipiv( k ) = -p
933 ipiv( k+1 ) = -kp
934 END IF
935*
936* Increase K and return to the start of the main loop
937*
938 k = k + kstep
939 GO TO 40
940*
941 64 CONTINUE
942*
943 END IF
944*
945 RETURN
946*
947* End of ZSYTF2_RK
948*
949 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine zsyr(UPLO, N, ALPHA, X, INCX, A, LDA)
ZSYR performs the symmetric rank-1 update of a complex symmetric matrix.
Definition: zsyr.f:135
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