LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ chetf2()

subroutine chetf2 ( character uplo,
integer n,
complex, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
integer info )

CHETF2 computes the factorization of a complex Hermitian matrix, using the diagonal pivoting method (unblocked algorithm calling Level 2 BLAS).

Download CHETF2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> CHETF2 computes the factorization of a complex Hermitian matrix A
!> using the Bunch-Kaufman diagonal pivoting method:
!>
!>    A = U*D*U**H  or  A = L*D*L**H
!>
!> where U (or L) is a product of permutation and unit upper (lower)
!> triangular matrices, U**H is the conjugate transpose of U, and D is
!> Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
!>
!> This is the unblocked version of the algorithm, calling Level 2 BLAS.
!> 
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          Specifies whether the upper or lower triangular part of the
!>          Hermitian matrix A is stored:
!>          = 'U':  Upper triangular
!>          = 'L':  Lower triangular
!> 
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!> 
[in,out]A
!>          A is COMPLEX array, dimension (LDA,N)
!>          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
!>          n-by-n upper triangular part of A contains the upper
!>          triangular part of the matrix A, and the strictly lower
!>          triangular part of A is not referenced.  If UPLO = 'L', the
!>          leading n-by-n lower triangular part of A contains the lower
!>          triangular part of the matrix A, and the strictly upper
!>          triangular part of A is not referenced.
!>
!>          On exit, the block diagonal matrix D and the multipliers used
!>          to obtain the factor U or L (see below for further details).
!> 
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!> 
[out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          Details of the interchanges and the block structure of D.
!>
!>          If UPLO = 'U':
!>             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
!>             interchanged and D(k,k) is a 1-by-1 diagonal block.
!>
!>             If IPIV(k) = IPIV(k-1) < 0, then rows and columns
!>             k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
!>             is a 2-by-2 diagonal block.
!>
!>          If UPLO = 'L':
!>             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
!>             interchanged and D(k,k) is a 1-by-1 diagonal block.
!>
!>             If IPIV(k) = IPIV(k+1) < 0, then rows and columns
!>             k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
!>             is a 2-by-2 diagonal block.
!> 
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -k, the k-th argument had an illegal value
!>          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
!>               has been completed, but the block diagonal matrix D is
!>               exactly singular, and division by zero will occur if it
!>               is used to solve a system of equations.
!> 
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  09-29-06 - patch from
!>    Bobby Cheng, MathWorks
!>
!>    Replace l.210 and l.392
!>         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
!>    by
!>         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN
!>
!>  01-01-96 - Based on modifications by
!>    J. Lewis, Boeing Computer Services Company
!>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
!>
!>  If UPLO = 'U', then A = U*D*U**H, where
!>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
!>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
!>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
!>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
!>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
!>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
!>
!>             (   I    v    0   )   k-s
!>     U(k) =  (   0    I    0   )   s
!>             (   0    0    I   )   n-k
!>                k-s   s   n-k
!>
!>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
!>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
!>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
!>
!>  If UPLO = 'L', then A = L*D*L**H, where
!>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
!>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
!>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
!>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
!>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
!>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
!>
!>             (   I    0     0   )  k-1
!>     L(k) =  (   0    I     0   )  s
!>             (   0    v     I   )  n-k-s+1
!>                k-1   s  n-k-s+1
!>
!>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
!>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
!>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
!> 

Definition at line 183 of file chetf2.f.

184*
185* -- LAPACK computational routine --
186* -- LAPACK is a software package provided by Univ. of Tennessee, --
187* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188*
189* .. Scalar Arguments ..
190 CHARACTER UPLO
191 INTEGER INFO, LDA, N
192* ..
193* .. Array Arguments ..
194 INTEGER IPIV( * )
195 COMPLEX A( LDA, * )
196* ..
197*
198* =====================================================================
199*
200* .. Parameters ..
201 REAL ZERO, ONE
202 parameter( zero = 0.0e+0, one = 1.0e+0 )
203 REAL EIGHT, SEVTEN
204 parameter( eight = 8.0e+0, sevten = 17.0e+0 )
205* ..
206* .. Local Scalars ..
207 LOGICAL UPPER
208 INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP
209 REAL ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
210 $ TT
211 COMPLEX D12, D21, T, WK, WKM1, WKP1, ZDUM
212* ..
213* .. External Functions ..
214 LOGICAL LSAME, SISNAN
215 INTEGER ICAMAX
216 REAL SLAPY2
217 EXTERNAL lsame, icamax, slapy2, sisnan
218* ..
219* .. External Subroutines ..
220 EXTERNAL cher, csscal, cswap, xerbla
221* ..
222* .. Intrinsic Functions ..
223 INTRINSIC abs, aimag, cmplx, conjg, max, real, sqrt
224* ..
225* .. Statement Functions ..
226 REAL CABS1
227* ..
228* .. Statement Function definitions ..
229 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
230* ..
231* .. Executable Statements ..
232*
233* Test the input parameters.
234*
235 info = 0
236 upper = lsame( uplo, 'U' )
237 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
238 info = -1
239 ELSE IF( n.LT.0 ) THEN
240 info = -2
241 ELSE IF( lda.LT.max( 1, n ) ) THEN
242 info = -4
243 END IF
244 IF( info.NE.0 ) THEN
245 CALL xerbla( 'CHETF2', -info )
246 RETURN
247 END IF
248*
249* Initialize ALPHA for use in choosing pivot block size.
250*
251 alpha = ( one+sqrt( sevten ) ) / eight
252*
253 IF( upper ) THEN
254*
255* Factorize A as U*D*U**H using the upper triangle of A
256*
257* K is the main loop index, decreasing from N to 1 in steps of
258* 1 or 2
259*
260 k = n
261 10 CONTINUE
262*
263* If K < 1, exit from loop
264*
265 IF( k.LT.1 )
266 $ GO TO 90
267 kstep = 1
268*
269* Determine rows and columns to be interchanged and whether
270* a 1-by-1 or 2-by-2 pivot block will be used
271*
272 absakk = abs( real( a( k, k ) ) )
273*
274* IMAX is the row-index of the largest off-diagonal element in
275* column K, and COLMAX is its absolute value.
276* Determine both COLMAX and IMAX.
277*
278 IF( k.GT.1 ) THEN
279 imax = icamax( k-1, a( 1, k ), 1 )
280 colmax = cabs1( a( imax, k ) )
281 ELSE
282 colmax = zero
283 END IF
284*
285 IF( (max( absakk, colmax ).EQ.zero) .OR.
286 $ sisnan(absakk) ) THEN
287*
288* Column K is or underflow, or contains a NaN:
289* set INFO and continue
290*
291 IF( info.EQ.0 )
292 $ info = k
293 kp = k
294 a( k, k ) = real( a( k, k ) )
295 ELSE
296 IF( absakk.GE.alpha*colmax ) THEN
297*
298* no interchange, use 1-by-1 pivot block
299*
300 kp = k
301 ELSE
302*
303* JMAX is the column-index of the largest off-diagonal
304* element in row IMAX, and ROWMAX is its absolute value
305*
306 jmax = imax + icamax( k-imax, a( imax, imax+1 ), lda )
307 rowmax = cabs1( a( imax, jmax ) )
308 IF( imax.GT.1 ) THEN
309 jmax = icamax( imax-1, a( 1, imax ), 1 )
310 rowmax = max( rowmax, cabs1( a( jmax, imax ) ) )
311 END IF
312*
313 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
314*
315* no interchange, use 1-by-1 pivot block
316*
317 kp = k
318 ELSE IF( abs( real( a( imax, imax ) ) ).GE.alpha*rowmax )
319 $ THEN
320*
321* interchange rows and columns K and IMAX, use 1-by-1
322* pivot block
323*
324 kp = imax
325 ELSE
326*
327* interchange rows and columns K-1 and IMAX, use 2-by-2
328* pivot block
329*
330 kp = imax
331 kstep = 2
332 END IF
333 END IF
334*
335 kk = k - kstep + 1
336 IF( kp.NE.kk ) THEN
337*
338* Interchange rows and columns KK and KP in the leading
339* submatrix A(1:k,1:k)
340*
341 CALL cswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
342 DO 20 j = kp + 1, kk - 1
343 t = conjg( a( j, kk ) )
344 a( j, kk ) = conjg( a( kp, j ) )
345 a( kp, j ) = t
346 20 CONTINUE
347 a( kp, kk ) = conjg( a( kp, kk ) )
348 r1 = real( a( kk, kk ) )
349 a( kk, kk ) = real( a( kp, kp ) )
350 a( kp, kp ) = r1
351 IF( kstep.EQ.2 ) THEN
352 a( k, k ) = real( a( k, k ) )
353 t = a( k-1, k )
354 a( k-1, k ) = a( kp, k )
355 a( kp, k ) = t
356 END IF
357 ELSE
358 a( k, k ) = real( a( k, k ) )
359 IF( kstep.EQ.2 )
360 $ a( k-1, k-1 ) = real( a( k-1, k-1 ) )
361 END IF
362*
363* Update the leading submatrix
364*
365 IF( kstep.EQ.1 ) THEN
366*
367* 1-by-1 pivot block D(k): column k now holds
368*
369* W(k) = U(k)*D(k)
370*
371* where U(k) is the k-th column of U
372*
373* Perform a rank-1 update of A(1:k-1,1:k-1) as
374*
375* A := A - U(k)*D(k)*U(k)**H = A - W(k)*1/D(k)*W(k)**H
376*
377 r1 = one / real( a( k, k ) )
378 CALL cher( uplo, k-1, -r1, a( 1, k ), 1, a, lda )
379*
380* Store U(k) in column k
381*
382 CALL csscal( k-1, r1, a( 1, k ), 1 )
383 ELSE
384*
385* 2-by-2 pivot block D(k): columns k and k-1 now hold
386*
387* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
388*
389* where U(k) and U(k-1) are the k-th and (k-1)-th columns
390* of U
391*
392* Perform a rank-2 update of A(1:k-2,1:k-2) as
393*
394* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**H
395* = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**H
396*
397 IF( k.GT.2 ) THEN
398*
399 d = slapy2( real( a( k-1, k ) ),
400 $ aimag( a( k-1, k ) ) )
401 d22 = real( a( k-1, k-1 ) ) / d
402 d11 = real( a( k, k ) ) / d
403 tt = one / ( d11*d22-one )
404 d12 = a( k-1, k ) / d
405 d = tt / d
406*
407 DO 40 j = k - 2, 1, -1
408 wkm1 = d*( d11*a( j, k-1 )-conjg( d12 )*a( j, k ) )
409 wk = d*( d22*a( j, k )-d12*a( j, k-1 ) )
410 DO 30 i = j, 1, -1
411 a( i, j ) = a( i, j ) - a( i, k )*conjg( wk ) -
412 $ a( i, k-1 )*conjg( wkm1 )
413 30 CONTINUE
414 a( j, k ) = wk
415 a( j, k-1 ) = wkm1
416 a( j, j ) = cmplx( real( a( j, j ) ), 0.0e+0 )
417 40 CONTINUE
418*
419 END IF
420*
421 END IF
422 END IF
423*
424* Store details of the interchanges in IPIV
425*
426 IF( kstep.EQ.1 ) THEN
427 ipiv( k ) = kp
428 ELSE
429 ipiv( k ) = -kp
430 ipiv( k-1 ) = -kp
431 END IF
432*
433* Decrease K and return to the start of the main loop
434*
435 k = k - kstep
436 GO TO 10
437*
438 ELSE
439*
440* Factorize A as L*D*L**H using the lower triangle of A
441*
442* K is the main loop index, increasing from 1 to N in steps of
443* 1 or 2
444*
445 k = 1
446 50 CONTINUE
447*
448* If K > N, exit from loop
449*
450 IF( k.GT.n )
451 $ GO TO 90
452 kstep = 1
453*
454* Determine rows and columns to be interchanged and whether
455* a 1-by-1 or 2-by-2 pivot block will be used
456*
457 absakk = abs( real( a( k, k ) ) )
458*
459* IMAX is the row-index of the largest off-diagonal element in
460* column K, and COLMAX is its absolute value.
461* Determine both COLMAX and IMAX.
462*
463 IF( k.LT.n ) THEN
464 imax = k + icamax( n-k, a( k+1, k ), 1 )
465 colmax = cabs1( a( imax, k ) )
466 ELSE
467 colmax = zero
468 END IF
469*
470 IF( (max( absakk, colmax ).EQ.zero) .OR.
471 $ sisnan(absakk) ) THEN
472*
473* Column K is zero or underflow, contains a NaN:
474* set INFO and continue
475*
476 IF( info.EQ.0 )
477 $ info = k
478 kp = k
479 a( k, k ) = real( a( k, k ) )
480 ELSE
481 IF( absakk.GE.alpha*colmax ) THEN
482*
483* no interchange, use 1-by-1 pivot block
484*
485 kp = k
486 ELSE
487*
488* JMAX is the column-index of the largest off-diagonal
489* element in row IMAX, and ROWMAX is its absolute value
490*
491 jmax = k - 1 + icamax( imax-k, a( imax, k ), lda )
492 rowmax = cabs1( a( imax, jmax ) )
493 IF( imax.LT.n ) THEN
494 jmax = imax + icamax( n-imax, a( imax+1, imax ),
495 $ 1 )
496 rowmax = max( rowmax, cabs1( a( jmax, imax ) ) )
497 END IF
498*
499 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) ) THEN
500*
501* no interchange, use 1-by-1 pivot block
502*
503 kp = k
504 ELSE IF( abs( real( a( imax, imax ) ) ).GE.alpha*rowmax )
505 $ THEN
506*
507* interchange rows and columns K and IMAX, use 1-by-1
508* pivot block
509*
510 kp = imax
511 ELSE
512*
513* interchange rows and columns K+1 and IMAX, use 2-by-2
514* pivot block
515*
516 kp = imax
517 kstep = 2
518 END IF
519 END IF
520*
521 kk = k + kstep - 1
522 IF( kp.NE.kk ) THEN
523*
524* Interchange rows and columns KK and KP in the trailing
525* submatrix A(k:n,k:n)
526*
527 IF( kp.LT.n )
528 $ CALL cswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ),
529 $ 1 )
530 DO 60 j = kk + 1, kp - 1
531 t = conjg( a( j, kk ) )
532 a( j, kk ) = conjg( a( kp, j ) )
533 a( kp, j ) = t
534 60 CONTINUE
535 a( kp, kk ) = conjg( a( kp, kk ) )
536 r1 = real( a( kk, kk ) )
537 a( kk, kk ) = real( a( kp, kp ) )
538 a( kp, kp ) = r1
539 IF( kstep.EQ.2 ) THEN
540 a( k, k ) = real( a( k, k ) )
541 t = a( k+1, k )
542 a( k+1, k ) = a( kp, k )
543 a( kp, k ) = t
544 END IF
545 ELSE
546 a( k, k ) = real( a( k, k ) )
547 IF( kstep.EQ.2 )
548 $ a( k+1, k+1 ) = real( a( k+1, k+1 ) )
549 END IF
550*
551* Update the trailing submatrix
552*
553 IF( kstep.EQ.1 ) THEN
554*
555* 1-by-1 pivot block D(k): column k now holds
556*
557* W(k) = L(k)*D(k)
558*
559* where L(k) is the k-th column of L
560*
561 IF( k.LT.n ) THEN
562*
563* Perform a rank-1 update of A(k+1:n,k+1:n) as
564*
565* A := A - L(k)*D(k)*L(k)**H = A - W(k)*(1/D(k))*W(k)**H
566*
567 r1 = one / real( a( k, k ) )
568 CALL cher( uplo, n-k, -r1, a( k+1, k ), 1,
569 $ a( k+1, k+1 ), lda )
570*
571* Store L(k) in column K
572*
573 CALL csscal( n-k, r1, a( k+1, k ), 1 )
574 END IF
575 ELSE
576*
577* 2-by-2 pivot block D(k)
578*
579 IF( k.LT.n-1 ) THEN
580*
581* Perform a rank-2 update of A(k+2:n,k+2:n) as
582*
583* A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**H
584* = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**H
585*
586* where L(k) and L(k+1) are the k-th and (k+1)-th
587* columns of L
588*
589 d = slapy2( real( a( k+1, k ) ),
590 $ aimag( a( k+1, k ) ) )
591 d11 = real( a( k+1, k+1 ) ) / d
592 d22 = real( a( k, k ) ) / d
593 tt = one / ( d11*d22-one )
594 d21 = a( k+1, k ) / d
595 d = tt / d
596*
597 DO 80 j = k + 2, n
598 wk = d*( d11*a( j, k )-d21*a( j, k+1 ) )
599 wkp1 = d*( d22*a( j, k+1 )-conjg( d21 )*a( j, k ) )
600 DO 70 i = j, n
601 a( i, j ) = a( i, j ) - a( i, k )*conjg( wk ) -
602 $ a( i, k+1 )*conjg( wkp1 )
603 70 CONTINUE
604 a( j, k ) = wk
605 a( j, k+1 ) = wkp1
606 a( j, j ) = cmplx( real( a( j, j ) ), 0.0e+0 )
607 80 CONTINUE
608 END IF
609 END IF
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 ) = -kp
618 ipiv( k+1 ) = -kp
619 END IF
620*
621* Increase K and return to the start of the main loop
622*
623 k = k + kstep
624 GO TO 50
625*
626 END IF
627*
628 90 CONTINUE
629 RETURN
630*
631* End of CHETF2
632*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cher(uplo, n, alpha, x, incx, a, lda)
CHER
Definition cher.f:135
integer function icamax(n, cx, incx)
ICAMAX
Definition icamax.f:71
logical function sisnan(sin)
SISNAN tests input for NaN.
Definition sisnan.f:57
real function slapy2(x, y)
SLAPY2 returns sqrt(x2+y2).
Definition slapy2.f:61
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: