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