LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
clahef_rook.f
Go to the documentation of this file.
1* \brief \b CLAHEF_ROOK computes a partial factorization of a complex Hermitian indefinite matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method (blocked algorithm, calling Level 3 BLAS).
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download CLAHEF_ROOK + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clahef_rook.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clahef_rook.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clahef_rook.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CLAHEF_ROOK( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
20*
21* .. Scalar Arguments ..
22* CHARACTER UPLO
23* INTEGER INFO, KB, LDA, LDW, N, NB
24* ..
25* .. Array Arguments ..
26* INTEGER IPIV( * )
27* COMPLEX A( LDA, * ), W( LDW, * )
28* ..
29*
30*
31*> \par Purpose:
32* =============
33*>
34*> \verbatim
35*>
36*> CLAHEF_ROOK computes a partial factorization of a complex Hermitian
37*> matrix A using the bounded Bunch-Kaufman ("rook") diagonal pivoting
38*> method. The partial factorization has the form:
39*>
40*> A = ( I U12 ) ( A11 0 ) ( I 0 ) if UPLO = 'U', or:
41*> ( 0 U22 ) ( 0 D ) ( U12**H U22**H )
42*>
43*> A = ( L11 0 ) ( D 0 ) ( L11**H L21**H ) if UPLO = 'L'
44*> ( L21 I ) ( 0 A22 ) ( 0 I )
45*>
46*> where the order of D is at most NB. The actual order is returned in
47*> the argument KB, and is either NB or NB-1, or N if N <= NB.
48*> Note that U**H denotes the conjugate transpose of U.
49*>
50*> CLAHEF_ROOK is an auxiliary routine called by CHETRF_ROOK. It uses
51*> blocked code (calling Level 3 BLAS) to update the submatrix
52*> A11 (if UPLO = 'U') or A22 (if UPLO = 'L').
53*> \endverbatim
54*
55* Arguments:
56* ==========
57*
58*> \param[in] UPLO
59*> \verbatim
60*> UPLO is CHARACTER*1
61*> Specifies whether the upper or lower triangular part of the
62*> Hermitian matrix A is stored:
63*> = 'U': Upper triangular
64*> = 'L': Lower triangular
65*> \endverbatim
66*>
67*> \param[in] N
68*> \verbatim
69*> N is INTEGER
70*> The order of the matrix A. N >= 0.
71*> \endverbatim
72*>
73*> \param[in] NB
74*> \verbatim
75*> NB is INTEGER
76*> The maximum number of columns of the matrix A that should be
77*> factored. NB should be at least 2 to allow for 2-by-2 pivot
78*> blocks.
79*> \endverbatim
80*>
81*> \param[out] KB
82*> \verbatim
83*> KB is INTEGER
84*> The number of columns of A that were actually factored.
85*> KB is either NB-1 or NB, or N if N <= NB.
86*> \endverbatim
87*>
88*> \param[in,out] A
89*> \verbatim
90*> A is COMPLEX array, dimension (LDA,N)
91*> On entry, the Hermitian matrix A. If UPLO = 'U', the leading
92*> n-by-n upper triangular part of A contains the upper
93*> triangular part of the matrix A, and the strictly lower
94*> triangular part of A is not referenced. If UPLO = 'L', the
95*> leading n-by-n lower triangular part of A contains the lower
96*> triangular part of the matrix A, and the strictly upper
97*> triangular part of A is not referenced.
98*> On exit, A contains details of the partial factorization.
99*> \endverbatim
100*>
101*> \param[in] LDA
102*> \verbatim
103*> LDA is INTEGER
104*> The leading dimension of the array A. LDA >= max(1,N).
105*> \endverbatim
106*>
107*> \param[out] IPIV
108*> \verbatim
109*> IPIV is INTEGER array, dimension (N)
110*> Details of the interchanges and the block structure of D.
111*>
112*> If UPLO = 'U':
113*> Only the last KB elements of IPIV are set.
114*>
115*> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
116*> interchanged and D(k,k) is a 1-by-1 diagonal block.
117*>
118*> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
119*> columns k and -IPIV(k) were interchanged and rows and
120*> columns k-1 and -IPIV(k-1) were inerchaged,
121*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
122*>
123*> If UPLO = 'L':
124*> Only the first KB elements of IPIV are set.
125*>
126*> If IPIV(k) > 0, then rows and columns k and IPIV(k)
127*> were interchanged and D(k,k) is a 1-by-1 diagonal block.
128*>
129*> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
130*> columns k and -IPIV(k) were interchanged and rows and
131*> columns k+1 and -IPIV(k+1) were inerchaged,
132*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
133*> \endverbatim
134*>
135*> \param[out] W
136*> \verbatim
137*> W is COMPLEX array, dimension (LDW,NB)
138*> \endverbatim
139*>
140*> \param[in] LDW
141*> \verbatim
142*> LDW is INTEGER
143*> The leading dimension of the array W. LDW >= max(1,N).
144*> \endverbatim
145*>
146*> \param[out] INFO
147*> \verbatim
148*> INFO is INTEGER
149*> = 0: successful exit
150*> > 0: if INFO = k, D(k,k) is exactly zero. The factorization
151*> has been completed, but the block diagonal matrix D is
152*> exactly singular.
153*> \endverbatim
154*
155* Authors:
156* ========
157*
158*> \author Univ. of Tennessee
159*> \author Univ. of California Berkeley
160*> \author Univ. of Colorado Denver
161*> \author NAG Ltd.
162*
163*> \ingroup lahef_rook
164*
165*> \par Contributors:
166* ==================
167*>
168*> \verbatim
169*>
170*> November 2013, Igor Kozachenko,
171*> Computer Science Division,
172*> University of California, Berkeley
173*>
174*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
175*> School of Mathematics,
176*> University of Manchester
177*> \endverbatim
178*
179* =====================================================================
180 SUBROUTINE clahef_rook( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
181 $ INFO )
182*
183* -- LAPACK computational routine --
184* -- LAPACK is a software package provided by Univ. of Tennessee, --
185* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
186*
187* .. Scalar Arguments ..
188 CHARACTER UPLO
189 INTEGER INFO, KB, LDA, LDW, N, NB
190* ..
191* .. Array Arguments ..
192 INTEGER IPIV( * )
193 COMPLEX A( LDA, * ), W( LDW, * )
194* ..
195*
196* =====================================================================
197*
198* .. Parameters ..
199 REAL ZERO, ONE
200 parameter( zero = 0.0e+0, one = 1.0e+0 )
201 COMPLEX CONE
202 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
203 REAL EIGHT, SEVTEN
204 parameter( eight = 8.0e+0, sevten = 17.0e+0 )
205* ..
206* .. Local Scalars ..
207 LOGICAL DONE
208 INTEGER IMAX, ITEMP, II, J, JB, JJ, JMAX, JP1, JP2, K,
209 $ kk, kkw, kp, kstep, kw, p
210 REAL ABSAKK, ALPHA, COLMAX, STEMP, R1, ROWMAX, T,
211 $ sfmin
212 COMPLEX D11, D21, D22, Z
213* ..
214* .. External Functions ..
215 LOGICAL LSAME
216 INTEGER ICAMAX
217 REAL SLAMCH
218 EXTERNAL lsame, icamax, slamch
219* ..
220* .. External Subroutines ..
221 EXTERNAL ccopy, csscal, cgemm, cgemv, clacgv,
222 $ cswap
223* ..
224* .. Intrinsic Functions ..
225 INTRINSIC abs, conjg, aimag, max, min, real, sqrt
226* ..
227* .. Statement Functions ..
228 REAL CABS1
229* ..
230* .. Statement Function definitions ..
231 cabs1( z ) = abs( real( z ) ) + abs( aimag( z ) )
232* ..
233* .. Executable Statements ..
234*
235 info = 0
236*
237* Initialize ALPHA for use in choosing pivot block size.
238*
239 alpha = ( one+sqrt( sevten ) ) / eight
240*
241* Compute machine safe minimum
242*
243 sfmin = slamch( 'S' )
244*
245 IF( lsame( uplo, 'U' ) ) THEN
246*
247* Factorize the trailing columns of A using the upper triangle
248* of A and working backwards, and compute the matrix W = U12*D
249* for use in updating A11 (note that conjg(W) is actually stored)
250*
251* K is the main loop index, decreasing from N in steps of 1 or 2
252*
253 k = n
254 10 CONTINUE
255*
256* KW is the column of W which corresponds to column K of A
257*
258 kw = nb + k - n
259*
260* Exit from loop
261*
262 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
263 $ GO TO 30
264*
265 kstep = 1
266 p = k
267*
268* Copy column K of A to column KW of W and update it
269*
270 IF( k.GT.1 )
271 $ CALL ccopy( k-1, a( 1, k ), 1, w( 1, kw ), 1 )
272 w( k, kw ) = real( a( k, k ) )
273 IF( k.LT.n ) THEN
274 CALL cgemv( 'No transpose', k, n-k, -cone, a( 1, k+1 ),
275 $ lda,
276 $ w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
277 w( k, kw ) = real( w( k, kw ) )
278 END IF
279*
280* Determine rows and columns to be interchanged and whether
281* a 1-by-1 or 2-by-2 pivot block will be used
282*
283 absakk = abs( real( w( k, kw ) ) )
284*
285* IMAX is the row-index of the largest off-diagonal element in
286* column K, and COLMAX is its absolute value.
287* Determine both COLMAX and IMAX.
288*
289 IF( k.GT.1 ) THEN
290 imax = icamax( k-1, w( 1, kw ), 1 )
291 colmax = cabs1( w( imax, kw ) )
292 ELSE
293 colmax = zero
294 END IF
295*
296 IF( max( absakk, colmax ).EQ.zero ) THEN
297*
298* Column K is zero or underflow: set INFO and continue
299*
300 IF( info.EQ.0 )
301 $ info = k
302 kp = k
303 a( k, k ) = real( w( k, kw ) )
304 IF( k.GT.1 )
305 $ CALL ccopy( k-1, w( 1, kw ), 1, a( 1, k ), 1 )
306 ELSE
307*
308* ============================================================
309*
310* BEGIN pivot search
311*
312* Case(1)
313* Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
314* (used to handle NaN and Inf)
315 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
316*
317* no interchange, use 1-by-1 pivot block
318*
319 kp = k
320*
321 ELSE
322*
323* Lop until pivot found
324*
325 done = .false.
326*
327 12 CONTINUE
328*
329* BEGIN pivot search loop body
330*
331*
332* Copy column IMAX to column KW-1 of W and update it
333*
334 IF( imax.GT.1 )
335 $ CALL ccopy( imax-1, a( 1, imax ), 1, w( 1,
336 $ kw-1 ),
337 $ 1 )
338 w( imax, kw-1 ) = real( a( imax, imax ) )
339*
340 CALL ccopy( k-imax, a( imax, imax+1 ), lda,
341 $ w( imax+1, kw-1 ), 1 )
342 CALL clacgv( k-imax, w( imax+1, kw-1 ), 1 )
343*
344 IF( k.LT.n ) THEN
345 CALL cgemv( 'No transpose', k, n-k, -cone,
346 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
347 $ cone, w( 1, kw-1 ), 1 )
348 w( imax, kw-1 ) = real( w( imax, kw-1 ) )
349 END IF
350*
351* JMAX is the column-index of the largest off-diagonal
352* element in row IMAX, and ROWMAX is its absolute value.
353* Determine both ROWMAX and JMAX.
354*
355 IF( imax.NE.k ) THEN
356 jmax = imax + icamax( k-imax, w( imax+1, kw-1 ),
357 $ 1 )
358 rowmax = cabs1( w( jmax, kw-1 ) )
359 ELSE
360 rowmax = zero
361 END IF
362*
363 IF( imax.GT.1 ) THEN
364 itemp = icamax( imax-1, w( 1, kw-1 ), 1 )
365 stemp = cabs1( w( itemp, kw-1 ) )
366 IF( stemp.GT.rowmax ) THEN
367 rowmax = stemp
368 jmax = itemp
369 END IF
370 END IF
371*
372* Case(2)
373* Equivalent to testing for
374* ABS( REAL( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
375* (used to handle NaN and Inf)
376*
377 IF( .NOT.( abs( real( w( imax,kw-1 ) ) )
378 $ .LT.alpha*rowmax ) ) THEN
379*
380* interchange rows and columns K and IMAX,
381* use 1-by-1 pivot block
382*
383 kp = imax
384*
385* copy column KW-1 of W to column KW of W
386*
387 CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
388*
389 done = .true.
390*
391* Case(3)
392* Equivalent to testing for ROWMAX.EQ.COLMAX,
393* (used to handle NaN and Inf)
394*
395 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
396 $ THEN
397*
398* interchange rows and columns K-1 and IMAX,
399* use 2-by-2 pivot block
400*
401 kp = imax
402 kstep = 2
403 done = .true.
404*
405* Case(4)
406 ELSE
407*
408* Pivot not found: set params and repeat
409*
410 p = imax
411 colmax = rowmax
412 imax = jmax
413*
414* Copy updated JMAXth (next IMAXth) column to Kth of W
415*
416 CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
417*
418 END IF
419*
420*
421* END pivot search loop body
422*
423 IF( .NOT.done ) GOTO 12
424*
425 END IF
426*
427* END pivot search
428*
429* ============================================================
430*
431* KK is the column of A where pivoting step stopped
432*
433 kk = k - kstep + 1
434*
435* KKW is the column of W which corresponds to column KK of A
436*
437 kkw = nb + kk - n
438*
439* Interchange rows and columns P and K.
440* Updated column P is already stored in column KW of W.
441*
442 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
443*
444* Copy non-updated column K to column P of submatrix A
445* at step K. No need to copy element into columns
446* K and K-1 of A for 2-by-2 pivot, since these columns
447* will be later overwritten.
448*
449 a( p, p ) = real( a( k, k ) )
450 CALL ccopy( k-1-p, a( p+1, k ), 1, a( p, p+1 ),
451 $ lda )
452 CALL clacgv( k-1-p, a( p, p+1 ), lda )
453 IF( p.GT.1 )
454 $ CALL ccopy( p-1, a( 1, k ), 1, a( 1, p ), 1 )
455*
456* Interchange rows K and P in the last K+1 to N columns of A
457* (columns K and K-1 of A for 2-by-2 pivot will be
458* later overwritten). Interchange rows K and P
459* in last KKW to NB columns of W.
460*
461 IF( k.LT.n )
462 $ CALL cswap( n-k, a( k, k+1 ), lda, a( p, k+1 ),
463 $ lda )
464 CALL cswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ),
465 $ ldw )
466 END IF
467*
468* Interchange rows and columns KP and KK.
469* Updated column KP is already stored in column KKW of W.
470*
471 IF( kp.NE.kk ) THEN
472*
473* Copy non-updated column KK to column KP of submatrix A
474* at step K. No need to copy element into column K
475* (or K and K-1 for 2-by-2 pivot) of A, since these columns
476* will be later overwritten.
477*
478 a( kp, kp ) = real( a( kk, kk ) )
479 CALL ccopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
480 $ lda )
481 CALL clacgv( kk-1-kp, a( kp, kp+1 ), lda )
482 IF( kp.GT.1 )
483 $ CALL ccopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
484*
485* Interchange rows KK and KP in last K+1 to N columns of A
486* (columns K (or K and K-1 for 2-by-2 pivot) of A will be
487* later overwritten). Interchange rows KK and KP
488* in last KKW to NB columns of W.
489*
490 IF( k.LT.n )
491 $ CALL cswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
492 $ lda )
493 CALL cswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
494 $ ldw )
495 END IF
496*
497 IF( kstep.EQ.1 ) THEN
498*
499* 1-by-1 pivot block D(k): column kw of W now holds
500*
501* W(kw) = U(k)*D(k),
502*
503* where U(k) is the k-th column of U
504*
505* (1) Store subdiag. elements of column U(k)
506* and 1-by-1 block D(k) in column k of A.
507* (NOTE: Diagonal element U(k,k) is a UNIT element
508* and not stored)
509* A(k,k) := D(k,k) = W(k,kw)
510* A(1:k-1,k) := U(1:k-1,k) = W(1:k-1,kw)/D(k,k)
511*
512* (NOTE: No need to use for Hermitian matrix
513* A( K, K ) = REAL( W( K, K) ) to separately copy diagonal
514* element D(k,k) from W (potentially saves only one load))
515 CALL ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
516 IF( k.GT.1 ) THEN
517*
518* (NOTE: No need to check if A(k,k) is NOT ZERO,
519* since that was ensured earlier in pivot search:
520* case A(k,k) = 0 falls into 2x2 pivot case(3))
521*
522* Handle division by a small number
523*
524 t = real( a( k, k ) )
525 IF( abs( t ).GE.sfmin ) THEN
526 r1 = one / t
527 CALL csscal( k-1, r1, a( 1, k ), 1 )
528 ELSE
529 DO 14 ii = 1, k-1
530 a( ii, k ) = a( ii, k ) / t
531 14 CONTINUE
532 END IF
533*
534* (2) Conjugate column W(kw)
535*
536 CALL clacgv( k-1, w( 1, kw ), 1 )
537 END IF
538*
539 ELSE
540*
541* 2-by-2 pivot block D(k): columns kw and kw-1 of W now hold
542*
543* ( W(kw-1) W(kw) ) = ( U(k-1) U(k) )*D(k)
544*
545* where U(k) and U(k-1) are the k-th and (k-1)-th columns
546* of U
547*
548* (1) Store U(1:k-2,k-1) and U(1:k-2,k) and 2-by-2
549* block D(k-1:k,k-1:k) in columns k-1 and k of A.
550* (NOTE: 2-by-2 diagonal block U(k-1:k,k-1:k) is a UNIT
551* block and not stored)
552* A(k-1:k,k-1:k) := D(k-1:k,k-1:k) = W(k-1:k,kw-1:kw)
553* A(1:k-2,k-1:k) := U(1:k-2,k:k-1:k) =
554* = W(1:k-2,kw-1:kw) * ( D(k-1:k,k-1:k)**(-1) )
555*
556 IF( k.GT.2 ) THEN
557*
558* Factor out the columns of the inverse of 2-by-2 pivot
559* block D, so that each column contains 1, to reduce the
560* number of FLOPS when we multiply panel
561* ( W(kw-1) W(kw) ) by this inverse, i.e. by D**(-1).
562*
563* D**(-1) = ( d11 cj(d21) )**(-1) =
564* ( d21 d22 )
565*
566* = 1/(d11*d22-|d21|**2) * ( ( d22) (-cj(d21) ) ) =
567* ( (-d21) ( d11 ) )
568*
569* = 1/(|d21|**2) * 1/((d11/cj(d21))*(d22/d21)-1) *
570*
571* * ( d21*( d22/d21 ) conj(d21)*( - 1 ) ) =
572* ( ( -1 ) ( d11/conj(d21) ) )
573*
574* = 1/(|d21|**2) * 1/(D22*D11-1) *
575*
576* * ( d21*( D11 ) conj(d21)*( -1 ) ) =
577* ( ( -1 ) ( D22 ) )
578*
579* = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) =
580* ( ( -1 ) ( D22 ) )
581*
582* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) =
583* ( ( -1 ) ( D22 ) )
584*
585* Handle division by a small number. (NOTE: order of
586* operations is important)
587*
588* = ( T*(( D11 )/conj(D21)) T*(( -1 )/D21 ) )
589* ( (( -1 ) ) (( D22 ) ) ),
590*
591* where D11 = d22/d21,
592* D22 = d11/conj(d21),
593* D21 = d21,
594* T = 1/(D22*D11-1).
595*
596* (NOTE: No need to check for division by ZERO,
597* since that was ensured earlier in pivot search:
598* (a) d21 != 0 in 2x2 pivot case(4),
599* since |d21| should be larger than |d11| and |d22|;
600* (b) (D22*D11 - 1) != 0, since from (a),
601* both |D11| < 1, |D22| < 1, hence |D22*D11| << 1.)
602*
603 d21 = w( k-1, kw )
604 d11 = w( k, kw ) / conjg( d21 )
605 d22 = w( k-1, kw-1 ) / d21
606 t = one / ( real( d11*d22 )-one )
607*
608* Update elements in columns A(k-1) and A(k) as
609* dot products of rows of ( W(kw-1) W(kw) ) and columns
610* of D**(-1)
611*
612 DO 20 j = 1, k - 2
613 a( j, k-1 ) = t*( ( d11*w( j, kw-1 )-w( j, kw ) ) /
614 $ d21 )
615 a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
616 $ conjg( d21 ) )
617 20 CONTINUE
618 END IF
619*
620* Copy D(k) to A
621*
622 a( k-1, k-1 ) = w( k-1, kw-1 )
623 a( k-1, k ) = w( k-1, kw )
624 a( k, k ) = w( k, kw )
625*
626* (2) Conjugate columns W(kw) and W(kw-1)
627*
628 CALL clacgv( k-1, w( 1, kw ), 1 )
629 CALL clacgv( k-2, w( 1, kw-1 ), 1 )
630*
631 END IF
632*
633 END IF
634*
635* Store details of the interchanges in IPIV
636*
637 IF( kstep.EQ.1 ) THEN
638 ipiv( k ) = kp
639 ELSE
640 ipiv( k ) = -p
641 ipiv( k-1 ) = -kp
642 END IF
643*
644* Decrease K and return to the start of the main loop
645*
646 k = k - kstep
647 GO TO 10
648*
649 30 CONTINUE
650*
651* Update the upper triangle of A11 (= A(1:k,1:k)) as
652*
653* A11 := A11 - U12*D*U12**H = A11 - U12*W**H
654*
655* computing blocks of NB columns at a time (note that conjg(W) is
656* actually stored)
657*
658 DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
659 jb = min( nb, k-j+1 )
660*
661* Update the upper triangle of the diagonal block
662*
663 DO 40 jj = j, j + jb - 1
664 a( jj, jj ) = real( a( jj, jj ) )
665 CALL cgemv( 'No transpose', jj-j+1, n-k, -cone,
666 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
667 $ a( j, jj ), 1 )
668 a( jj, jj ) = real( a( jj, jj ) )
669 40 CONTINUE
670*
671* Update the rectangular superdiagonal block
672*
673 IF( j.GE.2 )
674 $ CALL cgemm( 'No transpose', 'Transpose', j-1, jb, n-k,
675 $ -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
676 $ cone, a( 1, j ), lda )
677 50 CONTINUE
678*
679* Put U12 in standard form by partially undoing the interchanges
680* in of rows in columns k+1:n looping backwards from k+1 to n
681*
682 j = k + 1
683 60 CONTINUE
684*
685* Undo the interchanges (if any) of rows J and JP2
686* (or J and JP2, and J+1 and JP1) at each step J
687*
688 kstep = 1
689 jp1 = 1
690* (Here, J is a diagonal index)
691 jj = j
692 jp2 = ipiv( j )
693 IF( jp2.LT.0 ) THEN
694 jp2 = -jp2
695* (Here, J is a diagonal index)
696 j = j + 1
697 jp1 = -ipiv( j )
698 kstep = 2
699 END IF
700* (NOTE: Here, J is used to determine row length. Length N-J+1
701* of the rows to swap back doesn't include diagonal element)
702 j = j + 1
703 IF( jp2.NE.jj .AND. j.LE.n )
704 $ CALL cswap( n-j+1, a( jp2, j ), lda, a( jj, j ), lda )
705 jj = jj + 1
706 IF( kstep.EQ.2 .AND. jp1.NE.jj .AND. j.LE.n )
707 $ CALL cswap( n-j+1, a( jp1, j ), lda, a( jj, j ), lda )
708 IF( j.LT.n )
709 $ GO TO 60
710*
711* Set KB to the number of columns factorized
712*
713 kb = n - k
714*
715 ELSE
716*
717* Factorize the leading columns of A using the lower triangle
718* of A and working forwards, and compute the matrix W = L21*D
719* for use in updating A22 (note that conjg(W) is actually stored)
720*
721* K is the main loop index, increasing from 1 in steps of 1 or 2
722*
723 k = 1
724 70 CONTINUE
725*
726* Exit from loop
727*
728 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
729 $ GO TO 90
730*
731 kstep = 1
732 p = k
733*
734* Copy column K of A to column K of W and update column K of W
735*
736 w( k, k ) = real( a( k, k ) )
737 IF( k.LT.n )
738 $ CALL ccopy( n-k, a( k+1, k ), 1, w( k+1, k ), 1 )
739 IF( k.GT.1 ) THEN
740 CALL cgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
741 $ lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
742 w( k, k ) = real( w( k, k ) )
743 END IF
744*
745* Determine rows and columns to be interchanged and whether
746* a 1-by-1 or 2-by-2 pivot block will be used
747*
748 absakk = abs( real( w( k, k ) ) )
749*
750* IMAX is the row-index of the largest off-diagonal element in
751* column K, and COLMAX is its absolute value.
752* Determine both COLMAX and IMAX.
753*
754 IF( k.LT.n ) THEN
755 imax = k + icamax( n-k, w( k+1, k ), 1 )
756 colmax = cabs1( w( imax, k ) )
757 ELSE
758 colmax = zero
759 END IF
760*
761 IF( max( absakk, colmax ).EQ.zero ) THEN
762*
763* Column K is zero or underflow: set INFO and continue
764*
765 IF( info.EQ.0 )
766 $ info = k
767 kp = k
768 a( k, k ) = real( w( k, k ) )
769 IF( k.LT.n )
770 $ CALL ccopy( n-k, w( k+1, k ), 1, a( k+1, k ), 1 )
771 ELSE
772*
773* ============================================================
774*
775* BEGIN pivot search
776*
777* Case(1)
778* Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
779* (used to handle NaN and Inf)
780*
781 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
782*
783* no interchange, use 1-by-1 pivot block
784*
785 kp = k
786*
787 ELSE
788*
789 done = .false.
790*
791* Loop until pivot found
792*
793 72 CONTINUE
794*
795* BEGIN pivot search loop body
796*
797*
798* Copy column IMAX to column k+1 of W and update it
799*
800 CALL ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ),
801 $ 1)
802 CALL clacgv( imax-k, w( k, k+1 ), 1 )
803 w( imax, k+1 ) = real( a( imax, imax ) )
804*
805 IF( imax.LT.n )
806 $ CALL ccopy( n-imax, a( imax+1, imax ), 1,
807 $ w( imax+1, k+1 ), 1 )
808*
809 IF( k.GT.1 ) THEN
810 CALL cgemv( 'No transpose', n-k+1, k-1, -cone,
811 $ a( k, 1 ), lda, w( imax, 1 ), ldw,
812 $ cone, w( k, k+1 ), 1 )
813 w( imax, k+1 ) = real( w( imax, k+1 ) )
814 END IF
815*
816* JMAX is the column-index of the largest off-diagonal
817* element in row IMAX, and ROWMAX is its absolute value.
818* Determine both ROWMAX and JMAX.
819*
820 IF( imax.NE.k ) THEN
821 jmax = k - 1 + icamax( imax-k, w( k, k+1 ), 1 )
822 rowmax = cabs1( w( jmax, k+1 ) )
823 ELSE
824 rowmax = zero
825 END IF
826*
827 IF( imax.LT.n ) THEN
828 itemp = imax + icamax( n-imax, w( imax+1, k+1 ),
829 $ 1)
830 stemp = cabs1( w( itemp, k+1 ) )
831 IF( stemp.GT.rowmax ) THEN
832 rowmax = stemp
833 jmax = itemp
834 END IF
835 END IF
836*
837* Case(2)
838* Equivalent to testing for
839* ABS( REAL( W( IMAX,K+1 ) ) ).GE.ALPHA*ROWMAX
840* (used to handle NaN and Inf)
841*
842 IF( .NOT.( abs( real( w( imax,k+1 ) ) )
843 $ .LT.alpha*rowmax ) ) THEN
844*
845* interchange rows and columns K and IMAX,
846* use 1-by-1 pivot block
847*
848 kp = imax
849*
850* copy column K+1 of W to column K of W
851*
852 CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ),
853 $ 1 )
854*
855 done = .true.
856*
857* Case(3)
858* Equivalent to testing for ROWMAX.EQ.COLMAX,
859* (used to handle NaN and Inf)
860*
861 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
862 $ THEN
863*
864* interchange rows and columns K+1 and IMAX,
865* use 2-by-2 pivot block
866*
867 kp = imax
868 kstep = 2
869 done = .true.
870*
871* Case(4)
872 ELSE
873*
874* Pivot not found: set params and repeat
875*
876 p = imax
877 colmax = rowmax
878 imax = jmax
879*
880* Copy updated JMAXth (next IMAXth) column to Kth of W
881*
882 CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ),
883 $ 1 )
884*
885 END IF
886*
887*
888* End pivot search loop body
889*
890 IF( .NOT.done ) GOTO 72
891*
892 END IF
893*
894* END pivot search
895*
896* ============================================================
897*
898* KK is the column of A where pivoting step stopped
899*
900 kk = k + kstep - 1
901*
902* Interchange rows and columns P and K (only for 2-by-2 pivot).
903* Updated column P is already stored in column K of W.
904*
905 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
906*
907* Copy non-updated column KK-1 to column P of submatrix A
908* at step K. No need to copy element into columns
909* K and K+1 of A for 2-by-2 pivot, since these columns
910* will be later overwritten.
911*
912 a( p, p ) = real( a( k, k ) )
913 CALL ccopy( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
914 CALL clacgv( p-k-1, a( p, k+1 ), lda )
915 IF( p.LT.n )
916 $ CALL ccopy( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
917*
918* Interchange rows K and P in first K-1 columns of A
919* (columns K and K+1 of A for 2-by-2 pivot will be
920* later overwritten). Interchange rows K and P
921* in first KK columns of W.
922*
923 IF( k.GT.1 )
924 $ CALL cswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
925 CALL cswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
926 END IF
927*
928* Interchange rows and columns KP and KK.
929* Updated column KP is already stored in column KK of W.
930*
931 IF( kp.NE.kk ) THEN
932*
933* Copy non-updated column KK to column KP of submatrix A
934* at step K. No need to copy element into column K
935* (or K and K+1 for 2-by-2 pivot) of A, since these columns
936* will be later overwritten.
937*
938 a( kp, kp ) = real( a( kk, kk ) )
939 CALL ccopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
940 $ lda )
941 CALL clacgv( kp-kk-1, a( kp, kk+1 ), lda )
942 IF( kp.LT.n )
943 $ CALL ccopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ),
944 $ 1 )
945*
946* Interchange rows KK and KP in first K-1 columns of A
947* (column K (or K and K+1 for 2-by-2 pivot) of A will be
948* later overwritten). Interchange rows KK and KP
949* in first KK columns of W.
950*
951 IF( k.GT.1 )
952 $ CALL cswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
953 CALL cswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
954 END IF
955*
956 IF( kstep.EQ.1 ) THEN
957*
958* 1-by-1 pivot block D(k): column k of W now holds
959*
960* W(k) = L(k)*D(k),
961*
962* where L(k) is the k-th column of L
963*
964* (1) Store subdiag. elements of column L(k)
965* and 1-by-1 block D(k) in column k of A.
966* (NOTE: Diagonal element L(k,k) is a UNIT element
967* and not stored)
968* A(k,k) := D(k,k) = W(k,k)
969* A(k+1:N,k) := L(k+1:N,k) = W(k+1:N,k)/D(k,k)
970*
971* (NOTE: No need to use for Hermitian matrix
972* A( K, K ) = REAL( W( K, K) ) to separately copy diagonal
973* element D(k,k) from W (potentially saves only one load))
974 CALL ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
975 IF( k.LT.n ) THEN
976*
977* (NOTE: No need to check if A(k,k) is NOT ZERO,
978* since that was ensured earlier in pivot search:
979* case A(k,k) = 0 falls into 2x2 pivot case(3))
980*
981* Handle division by a small number
982*
983 t = real( a( k, k ) )
984 IF( abs( t ).GE.sfmin ) THEN
985 r1 = one / t
986 CALL csscal( n-k, r1, a( k+1, k ), 1 )
987 ELSE
988 DO 74 ii = k + 1, n
989 a( ii, k ) = a( ii, k ) / t
990 74 CONTINUE
991 END IF
992*
993* (2) Conjugate column W(k)
994*
995 CALL clacgv( n-k, w( k+1, k ), 1 )
996 END IF
997*
998 ELSE
999*
1000* 2-by-2 pivot block D(k): columns k and k+1 of W now hold
1001*
1002* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
1003*
1004* where L(k) and L(k+1) are the k-th and (k+1)-th columns
1005* of L
1006*
1007* (1) Store L(k+2:N,k) and L(k+2:N,k+1) and 2-by-2
1008* block D(k:k+1,k:k+1) in columns k and k+1 of A.
1009* NOTE: 2-by-2 diagonal block L(k:k+1,k:k+1) is a UNIT
1010* block and not stored.
1011* A(k:k+1,k:k+1) := D(k:k+1,k:k+1) = W(k:k+1,k:k+1)
1012* A(k+2:N,k:k+1) := L(k+2:N,k:k+1) =
1013* = W(k+2:N,k:k+1) * ( D(k:k+1,k:k+1)**(-1) )
1014*
1015 IF( k.LT.n-1 ) THEN
1016*
1017* Factor out the columns of the inverse of 2-by-2 pivot
1018* block D, so that each column contains 1, to reduce the
1019* number of FLOPS when we multiply panel
1020* ( W(kw-1) W(kw) ) by this inverse, i.e. by D**(-1).
1021*
1022* D**(-1) = ( d11 cj(d21) )**(-1) =
1023* ( d21 d22 )
1024*
1025* = 1/(d11*d22-|d21|**2) * ( ( d22) (-cj(d21) ) ) =
1026* ( (-d21) ( d11 ) )
1027*
1028* = 1/(|d21|**2) * 1/((d11/cj(d21))*(d22/d21)-1) *
1029*
1030* * ( d21*( d22/d21 ) conj(d21)*( - 1 ) ) =
1031* ( ( -1 ) ( d11/conj(d21) ) )
1032*
1033* = 1/(|d21|**2) * 1/(D22*D11-1) *
1034*
1035* * ( d21*( D11 ) conj(d21)*( -1 ) ) =
1036* ( ( -1 ) ( D22 ) )
1037*
1038* = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) =
1039* ( ( -1 ) ( D22 ) )
1040*
1041* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) =
1042* ( ( -1 ) ( D22 ) )
1043*
1044* Handle division by a small number. (NOTE: order of
1045* operations is important)
1046*
1047* = ( T*(( D11 )/conj(D21)) T*(( -1 )/D21 ) )
1048* ( (( -1 ) ) (( D22 ) ) ),
1049*
1050* where D11 = d22/d21,
1051* D22 = d11/conj(d21),
1052* D21 = d21,
1053* T = 1/(D22*D11-1).
1054*
1055* (NOTE: No need to check for division by ZERO,
1056* since that was ensured earlier in pivot search:
1057* (a) d21 != 0 in 2x2 pivot case(4),
1058* since |d21| should be larger than |d11| and |d22|;
1059* (b) (D22*D11 - 1) != 0, since from (a),
1060* both |D11| < 1, |D22| < 1, hence |D22*D11| << 1.)
1061*
1062 d21 = w( k+1, k )
1063 d11 = w( k+1, k+1 ) / d21
1064 d22 = w( k, k ) / conjg( d21 )
1065 t = one / ( real( d11*d22 )-one )
1066*
1067* Update elements in columns A(k) and A(k+1) as
1068* dot products of rows of ( W(k) W(k+1) ) and columns
1069* of D**(-1)
1070*
1071 DO 80 j = k + 2, n
1072 a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
1073 $ conjg( d21 ) )
1074 a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
1075 $ d21 )
1076 80 CONTINUE
1077 END IF
1078*
1079* Copy D(k) to A
1080*
1081 a( k, k ) = w( k, k )
1082 a( k+1, k ) = w( k+1, k )
1083 a( k+1, k+1 ) = w( k+1, k+1 )
1084*
1085* (2) Conjugate columns W(k) and W(k+1)
1086*
1087 CALL clacgv( n-k, w( k+1, k ), 1 )
1088 CALL clacgv( n-k-1, w( k+2, k+1 ), 1 )
1089*
1090 END IF
1091*
1092 END IF
1093*
1094* Store details of the interchanges in IPIV
1095*
1096 IF( kstep.EQ.1 ) THEN
1097 ipiv( k ) = kp
1098 ELSE
1099 ipiv( k ) = -p
1100 ipiv( k+1 ) = -kp
1101 END IF
1102*
1103* Increase K and return to the start of the main loop
1104*
1105 k = k + kstep
1106 GO TO 70
1107*
1108 90 CONTINUE
1109*
1110* Update the lower triangle of A22 (= A(k:n,k:n)) as
1111*
1112* A22 := A22 - L21*D*L21**H = A22 - L21*W**H
1113*
1114* computing blocks of NB columns at a time (note that conjg(W) is
1115* actually stored)
1116*
1117 DO 110 j = k, n, nb
1118 jb = min( nb, n-j+1 )
1119*
1120* Update the lower triangle of the diagonal block
1121*
1122 DO 100 jj = j, j + jb - 1
1123 a( jj, jj ) = real( a( jj, jj ) )
1124 CALL cgemv( 'No transpose', j+jb-jj, k-1, -cone,
1125 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
1126 $ a( jj, jj ), 1 )
1127 a( jj, jj ) = real( a( jj, jj ) )
1128 100 CONTINUE
1129*
1130* Update the rectangular subdiagonal block
1131*
1132 IF( j+jb.LE.n )
1133 $ CALL cgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
1134 $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
1135 $ ldw, cone, a( j+jb, j ), lda )
1136 110 CONTINUE
1137*
1138* Put L21 in standard form by partially undoing the interchanges
1139* of rows in columns 1:k-1 looping backwards from k-1 to 1
1140*
1141 j = k - 1
1142 120 CONTINUE
1143*
1144* Undo the interchanges (if any) of rows J and JP2
1145* (or J and JP2, and J-1 and JP1) at each step J
1146*
1147 kstep = 1
1148 jp1 = 1
1149* (Here, J is a diagonal index)
1150 jj = j
1151 jp2 = ipiv( j )
1152 IF( jp2.LT.0 ) THEN
1153 jp2 = -jp2
1154* (Here, J is a diagonal index)
1155 j = j - 1
1156 jp1 = -ipiv( j )
1157 kstep = 2
1158 END IF
1159* (NOTE: Here, J is used to determine row length. Length J
1160* of the rows to swap back doesn't include diagonal element)
1161 j = j - 1
1162 IF( jp2.NE.jj .AND. j.GE.1 )
1163 $ CALL cswap( j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
1164 jj = jj -1
1165 IF( kstep.EQ.2 .AND. jp1.NE.jj .AND. j.GE.1 )
1166 $ CALL cswap( j, a( jp1, 1 ), lda, a( jj, 1 ), lda )
1167 IF( j.GT.1 )
1168 $ GO TO 120
1169*
1170* Set KB to the number of columns factorized
1171*
1172 kb = k - 1
1173*
1174 END IF
1175 RETURN
1176*
1177* End of CLAHEF_ROOK
1178*
1179 END
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
clacgv
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:72
clahef_rook
subroutine clahef_rook(uplo, n, nb, kb, a, lda, ipiv, w, ldw, info)
Download CLAHEF_ROOK + dependencies [TGZ] [ZIP] [TXT]
Definition clahef_rook.f:182
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81