LAPACK: SRC/zhetf2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZHETF2( UPLO, N, A, LDA, IPIV, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            INFO, LDA, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            IPIV( * )
00014       COMPLEX*16         A( LDA, * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  ZHETF2 computes the factorization of a complex Hermitian matrix A
00021 *  using the Bunch-Kaufman diagonal pivoting method:
00022 *
00023 *     A = U*D*U'  or  A = L*D*L'
00024 *
00025 *  where U (or L) is a product of permutation and unit upper (lower)
00026 *  triangular matrices, U' is the conjugate transpose of U, and D is
00027 *  Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
00028 *
00029 *  This is the unblocked version of the algorithm, calling Level 2 BLAS.
00030 *
00031 *  Arguments
00032 *  =========
00033 *
00034 *  UPLO    (input) CHARACTER*1
00035 *          Specifies whether the upper or lower triangular part of the
00036 *          Hermitian matrix A is stored:
00037 *          = 'U':  Upper triangular
00038 *          = 'L':  Lower triangular
00039 *
00040 *  N       (input) INTEGER
00041 *          The order of the matrix A.  N >= 0.
00042 *
00043 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
00044 *          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
00045 *          n-by-n upper triangular part of A contains the upper
00046 *          triangular part of the matrix A, and the strictly lower
00047 *          triangular part of A is not referenced.  If UPLO = 'L', the
00048 *          leading n-by-n lower triangular part of A contains the lower
00049 *          triangular part of the matrix A, and the strictly upper
00050 *          triangular part of A is not referenced.
00051 *
00052 *          On exit, the block diagonal matrix D and the multipliers used
00053 *          to obtain the factor U or L (see below for further details).
00054 *
00055 *  LDA     (input) INTEGER
00056 *          The leading dimension of the array A.  LDA >= max(1,N).
00057 *
00058 *  IPIV    (output) INTEGER array, dimension (N)
00059 *          Details of the interchanges and the block structure of D.
00060 *          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
00061 *          interchanged and D(k,k) is a 1-by-1 diagonal block.
00062 *          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
00063 *          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
00064 *          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
00065 *          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
00066 *          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
00067 *
00068 *  INFO    (output) INTEGER
00069 *          = 0: successful exit
00070 *          < 0: if INFO = -k, the k-th argument had an illegal value
00071 *          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
00072 *               has been completed, but the block diagonal matrix D is
00073 *               exactly singular, and division by zero will occur if it
00074 *               is used to solve a system of equations.
00075 *
00076 *  Further Details
00077 *  ===============
00078 *
00079 *  09-29-06 - patch from
00080 *    Bobby Cheng, MathWorks
00081 *
00082 *    Replace l.210 and l.393
00083 *         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
00084 *    by
00085 *         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
00086 *
00087 *  01-01-96 - Based on modifications by
00088 *    J. Lewis, Boeing Computer Services Company
00089 *    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
00090 *
00091 *  If UPLO = 'U', then A = U*D*U', where
00092 *     U = P(n)*U(n)* ... *P(k)U(k)* ...,
00093 *  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
00094 *  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
00095 *  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
00096 *  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
00097 *  that if the diagonal block D(k) is of order s (s = 1 or 2), then
00098 *
00099 *             (   I    v    0   )   k-s
00100 *     U(k) =  (   0    I    0   )   s
00101 *             (   0    0    I   )   n-k
00102 *                k-s   s   n-k
00103 *
00104 *  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
00105 *  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
00106 *  and A(k,k), and v overwrites A(1:k-2,k-1:k).
00107 *
00108 *  If UPLO = 'L', then A = L*D*L', where
00109 *     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
00110 *  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
00111 *  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
00112 *  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
00113 *  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
00114 *  that if the diagonal block D(k) is of order s (s = 1 or 2), then
00115 *
00116 *             (   I    0     0   )  k-1
00117 *     L(k) =  (   0    I     0   )  s
00118 *             (   0    v     I   )  n-k-s+1
00119 *                k-1   s  n-k-s+1
00120 *
00121 *  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
00122 *  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
00123 *  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
00124 *
00125 *  =====================================================================
00126 *
00127 *     .. Parameters ..
00128       DOUBLE PRECISION   ZERO, ONE
00129       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00130       DOUBLE PRECISION   EIGHT, SEVTEN
00131       PARAMETER          ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
00132 *     ..
00133 *     .. Local Scalars ..
00134       LOGICAL            UPPER
00135       INTEGER            I, IMAX, J, JMAX, K, KK, KP, KSTEP
00136       DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
00137      $                   TT
00138       COMPLEX*16         D12, D21, T, WK, WKM1, WKP1, ZDUM
00139 *     ..
00140 *     .. External Functions ..
00141       LOGICAL            LSAME, DISNAN
00142       INTEGER            IZAMAX
00143       DOUBLE PRECISION   DLAPY2
00144       EXTERNAL           LSAME, IZAMAX, DLAPY2, DISNAN
00145 *     ..
00146 *     .. External Subroutines ..
00147       EXTERNAL           XERBLA, ZDSCAL, ZHER, ZSWAP
00148 *     ..
00149 *     .. Intrinsic Functions ..
00150       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
00151 *     ..
00152 *     .. Statement Functions ..
00153       DOUBLE PRECISION   CABS1
00154 *     ..
00155 *     .. Statement Function definitions ..
00156       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
00157 *     ..
00158 *     .. Executable Statements ..
00159 *
00160 *     Test the input parameters.
00161 *
00162       INFO = 0
00163       UPPER = LSAME( UPLO, 'U' )
00164       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00165          INFO = -1
00166       ELSE IF( N.LT.0 ) THEN
00167          INFO = -2
00168       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00169          INFO = -4
00170       END IF
00171       IF( INFO.NE.0 ) THEN
00172          CALL XERBLA( 'ZHETF2', -INFO )
00173          RETURN
00174       END IF
00175 *
00176 *     Initialize ALPHA for use in choosing pivot block size.
00177 *
00178       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
00179 *
00180       IF( UPPER ) THEN
00181 *
00182 *        Factorize A as U*D*U' using the upper triangle of A
00183 *
00184 *        K is the main loop index, decreasing from N to 1 in steps of
00185 *        1 or 2
00186 *
00187          K = N
00188    10    CONTINUE
00189 *
00190 *        If K < 1, exit from loop
00191 *
00192          IF( K.LT.1 )
00193      $      GO TO 90
00194          KSTEP = 1
00195 *
00196 *        Determine rows and columns to be interchanged and whether
00197 *        a 1-by-1 or 2-by-2 pivot block will be used
00198 *
00199          ABSAKK = ABS( DBLE( A( K, K ) ) )
00200 *
00201 *        IMAX is the row-index of the largest off-diagonal element in
00202 *        column K, and COLMAX is its absolute value
00203 *
00204          IF( K.GT.1 ) THEN
00205             IMAX = IZAMAX( K-1, A( 1, K ), 1 )
00206             COLMAX = CABS1( A( IMAX, K ) )
00207          ELSE
00208             COLMAX = ZERO
00209          END IF
00210 *
00211          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
00212 *
00213 *           Column K is zero or contains a NaN: set INFO and continue
00214 *
00215             IF( INFO.EQ.0 )
00216      $         INFO = K
00217             KP = K
00218             A( K, K ) = DBLE( A( K, K ) )
00219          ELSE
00220             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
00221 *
00222 *              no interchange, use 1-by-1 pivot block
00223 *
00224                KP = K
00225             ELSE
00226 *
00227 *              JMAX is the column-index of the largest off-diagonal
00228 *              element in row IMAX, and ROWMAX is its absolute value
00229 *
00230                JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
00231                ROWMAX = CABS1( A( IMAX, JMAX ) )
00232                IF( IMAX.GT.1 ) THEN
00233                   JMAX = IZAMAX( IMAX-1, A( 1, IMAX ), 1 )
00234                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
00235                END IF
00236 *
00237                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
00238 *
00239 *                 no interchange, use 1-by-1 pivot block
00240 *
00241                   KP = K
00242                ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
00243      $                   THEN
00244 *
00245 *                 interchange rows and columns K and IMAX, use 1-by-1
00246 *                 pivot block
00247 *
00248                   KP = IMAX
00249                ELSE
00250 *
00251 *                 interchange rows and columns K-1 and IMAX, use 2-by-2
00252 *                 pivot block
00253 *
00254                   KP = IMAX
00255                   KSTEP = 2
00256                END IF
00257             END IF
00258 *
00259             KK = K - KSTEP + 1
00260             IF( KP.NE.KK ) THEN
00261 *
00262 *              Interchange rows and columns KK and KP in the leading
00263 *              submatrix A(1:k,1:k)
00264 *
00265                CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
00266                DO 20 J = KP + 1, KK - 1
00267                   T = DCONJG( A( J, KK ) )
00268                   A( J, KK ) = DCONJG( A( KP, J ) )
00269                   A( KP, J ) = T
00270    20          CONTINUE
00271                A( KP, KK ) = DCONJG( A( KP, KK ) )
00272                R1 = DBLE( A( KK, KK ) )
00273                A( KK, KK ) = DBLE( A( KP, KP ) )
00274                A( KP, KP ) = R1
00275                IF( KSTEP.EQ.2 ) THEN
00276                   A( K, K ) = DBLE( A( K, K ) )
00277                   T = A( K-1, K )
00278                   A( K-1, K ) = A( KP, K )
00279                   A( KP, K ) = T
00280                END IF
00281             ELSE
00282                A( K, K ) = DBLE( A( K, K ) )
00283                IF( KSTEP.EQ.2 )
00284      $            A( K-1, K-1 ) = DBLE( A( K-1, K-1 ) )
00285             END IF
00286 *
00287 *           Update the leading submatrix
00288 *
00289             IF( KSTEP.EQ.1 ) THEN
00290 *
00291 *              1-by-1 pivot block D(k): column k now holds
00292 *
00293 *              W(k) = U(k)*D(k)
00294 *
00295 *              where U(k) is the k-th column of U
00296 *
00297 *              Perform a rank-1 update of A(1:k-1,1:k-1) as
00298 *
00299 *              A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)'
00300 *
00301                R1 = ONE / DBLE( A( K, K ) )
00302                CALL ZHER( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
00303 *
00304 *              Store U(k) in column k
00305 *
00306                CALL ZDSCAL( K-1, R1, A( 1, K ), 1 )
00307             ELSE
00308 *
00309 *              2-by-2 pivot block D(k): columns k and k-1 now hold
00310 *
00311 *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
00312 *
00313 *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
00314 *              of U
00315 *
00316 *              Perform a rank-2 update of A(1:k-2,1:k-2) as
00317 *
00318 *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )'
00319 *                 = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )'
00320 *
00321                IF( K.GT.2 ) THEN
00322 *
00323                   D = DLAPY2( DBLE( A( K-1, K ) ),
00324      $                DIMAG( A( K-1, K ) ) )
00325                   D22 = DBLE( A( K-1, K-1 ) ) / D
00326                   D11 = DBLE( A( K, K ) ) / D
00327                   TT = ONE / ( D11*D22-ONE )
00328                   D12 = A( K-1, K ) / D
00329                   D = TT / D
00330 *
00331                   DO 40 J = K - 2, 1, -1
00332                      WKM1 = D*( D11*A( J, K-1 )-DCONJG( D12 )*
00333      $                      A( J, K ) )
00334                      WK = D*( D22*A( J, K )-D12*A( J, K-1 ) )
00335                      DO 30 I = J, 1, -1
00336                         A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
00337      $                              A( I, K-1 )*DCONJG( WKM1 )
00338    30                CONTINUE
00339                      A( J, K ) = WK
00340                      A( J, K-1 ) = WKM1
00341                      A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
00342    40             CONTINUE
00343 *
00344                END IF
00345 *
00346             END IF
00347          END IF
00348 *
00349 *        Store details of the interchanges in IPIV
00350 *
00351          IF( KSTEP.EQ.1 ) THEN
00352             IPIV( K ) = KP
00353          ELSE
00354             IPIV( K ) = -KP
00355             IPIV( K-1 ) = -KP
00356          END IF
00357 *
00358 *        Decrease K and return to the start of the main loop
00359 *
00360          K = K - KSTEP
00361          GO TO 10
00362 *
00363       ELSE
00364 *
00365 *        Factorize A as L*D*L' using the lower triangle of A
00366 *
00367 *        K is the main loop index, increasing from 1 to N in steps of
00368 *        1 or 2
00369 *
00370          K = 1
00371    50    CONTINUE
00372 *
00373 *        If K > N, exit from loop
00374 *
00375          IF( K.GT.N )
00376      $      GO TO 90
00377          KSTEP = 1
00378 *
00379 *        Determine rows and columns to be interchanged and whether
00380 *        a 1-by-1 or 2-by-2 pivot block will be used
00381 *
00382          ABSAKK = ABS( DBLE( A( K, K ) ) )
00383 *
00384 *        IMAX is the row-index of the largest off-diagonal element in
00385 *        column K, and COLMAX is its absolute value
00386 *
00387          IF( K.LT.N ) THEN
00388             IMAX = K + IZAMAX( N-K, A( K+1, K ), 1 )
00389             COLMAX = CABS1( A( IMAX, K ) )
00390          ELSE
00391             COLMAX = ZERO
00392          END IF
00393 *
00394          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
00395 *
00396 *           Column K is zero or contains a NaN: set INFO and continue
00397 *
00398             IF( INFO.EQ.0 )
00399      $         INFO = K
00400             KP = K
00401             A( K, K ) = DBLE( A( K, K ) )
00402          ELSE
00403             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
00404 *
00405 *              no interchange, use 1-by-1 pivot block
00406 *
00407                KP = K
00408             ELSE
00409 *
00410 *              JMAX is the column-index of the largest off-diagonal
00411 *              element in row IMAX, and ROWMAX is its absolute value
00412 *
00413                JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
00414                ROWMAX = CABS1( A( IMAX, JMAX ) )
00415                IF( IMAX.LT.N ) THEN
00416                   JMAX = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
00417                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
00418                END IF
00419 *
00420                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
00421 *
00422 *                 no interchange, use 1-by-1 pivot block
00423 *
00424                   KP = K
00425                ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
00426      $                   THEN
00427 *
00428 *                 interchange rows and columns K and IMAX, use 1-by-1
00429 *                 pivot block
00430 *
00431                   KP = IMAX
00432                ELSE
00433 *
00434 *                 interchange rows and columns K+1 and IMAX, use 2-by-2
00435 *                 pivot block
00436 *
00437                   KP = IMAX
00438                   KSTEP = 2
00439                END IF
00440             END IF
00441 *
00442             KK = K + KSTEP - 1
00443             IF( KP.NE.KK ) THEN
00444 *
00445 *              Interchange rows and columns KK and KP in the trailing
00446 *              submatrix A(k:n,k:n)
00447 *
00448                IF( KP.LT.N )
00449      $            CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
00450                DO 60 J = KK + 1, KP - 1
00451                   T = DCONJG( A( J, KK ) )
00452                   A( J, KK ) = DCONJG( A( KP, J ) )
00453                   A( KP, J ) = T
00454    60          CONTINUE
00455                A( KP, KK ) = DCONJG( A( KP, KK ) )
00456                R1 = DBLE( A( KK, KK ) )
00457                A( KK, KK ) = DBLE( A( KP, KP ) )
00458                A( KP, KP ) = R1
00459                IF( KSTEP.EQ.2 ) THEN
00460                   A( K, K ) = DBLE( A( K, K ) )
00461                   T = A( K+1, K )
00462                   A( K+1, K ) = A( KP, K )
00463                   A( KP, K ) = T
00464                END IF
00465             ELSE
00466                A( K, K ) = DBLE( A( K, K ) )
00467                IF( KSTEP.EQ.2 )
00468      $            A( K+1, K+1 ) = DBLE( A( K+1, K+1 ) )
00469             END IF
00470 *
00471 *           Update the trailing submatrix
00472 *
00473             IF( KSTEP.EQ.1 ) THEN
00474 *
00475 *              1-by-1 pivot block D(k): column k now holds
00476 *
00477 *              W(k) = L(k)*D(k)
00478 *
00479 *              where L(k) is the k-th column of L
00480 *
00481                IF( K.LT.N ) THEN
00482 *
00483 *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
00484 *
00485 *                 A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)'
00486 *
00487                   R1 = ONE / DBLE( A( K, K ) )
00488                   CALL ZHER( UPLO, N-K, -R1, A( K+1, K ), 1,
00489      $                       A( K+1, K+1 ), LDA )
00490 *
00491 *                 Store L(k) in column K
00492 *
00493                   CALL ZDSCAL( N-K, R1, A( K+1, K ), 1 )
00494                END IF
00495             ELSE
00496 *
00497 *              2-by-2 pivot block D(k)
00498 *
00499                IF( K.LT.N-1 ) THEN
00500 *
00501 *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
00502 *
00503 *                 A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )'
00504 *                    = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )'
00505 *
00506 *                 where L(k) and L(k+1) are the k-th and (k+1)-th
00507 *                 columns of L
00508 *
00509                   D = DLAPY2( DBLE( A( K+1, K ) ),
00510      $                DIMAG( A( K+1, K ) ) )
00511                   D11 = DBLE( A( K+1, K+1 ) ) / D
00512                   D22 = DBLE( A( K, K ) ) / D
00513                   TT = ONE / ( D11*D22-ONE )
00514                   D21 = A( K+1, K ) / D
00515                   D = TT / D
00516 *
00517                   DO 80 J = K + 2, N
00518                      WK = D*( D11*A( J, K )-D21*A( J, K+1 ) )
00519                      WKP1 = D*( D22*A( J, K+1 )-DCONJG( D21 )*
00520      $                      A( J, K ) )
00521                      DO 70 I = J, N
00522                         A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
00523      $                              A( I, K+1 )*DCONJG( WKP1 )
00524    70                CONTINUE
00525                      A( J, K ) = WK
00526                      A( J, K+1 ) = WKP1
00527                      A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
00528    80             CONTINUE
00529                END IF
00530             END IF
00531          END IF
00532 *
00533 *        Store details of the interchanges in IPIV
00534 *
00535          IF( KSTEP.EQ.1 ) THEN
00536             IPIV( K ) = KP
00537          ELSE
00538             IPIV( K ) = -KP
00539             IPIV( K+1 ) = -KP
00540          END IF
00541 *
00542 *        Increase K and return to the start of the main loop
00543 *
00544          K = K + KSTEP
00545          GO TO 50
00546 *
00547       END IF
00548 *
00549    90 CONTINUE
00550       RETURN
00551 *
00552 *     End of ZHETF2
00553 *
00554       END