LAPACK: SRC/chetf2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CHETF2( 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            A( LDA, * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  CHETF2 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 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.392
00083 *         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
00084 *    by
00085 *         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(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       REAL               ZERO, ONE
00129       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00130       REAL               EIGHT, SEVTEN
00131       PARAMETER          ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
00132 *     ..
00133 *     .. Local Scalars ..
00134       LOGICAL            UPPER
00135       INTEGER            I, IMAX, J, JMAX, K, KK, KP, KSTEP
00136       REAL               ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
00137      $                   TT
00138       COMPLEX            D12, D21, T, WK, WKM1, WKP1, ZDUM
00139 *     ..
00140 *     .. External Functions ..
00141       LOGICAL            LSAME, SISNAN
00142       INTEGER            ICAMAX
00143       REAL               SLAPY2
00144       EXTERNAL           LSAME, ICAMAX, SLAPY2, SISNAN
00145 *     ..
00146 *     .. External Subroutines ..
00147       EXTERNAL           CHER, CSSCAL, CSWAP, XERBLA
00148 *     ..
00149 *     .. Intrinsic Functions ..
00150       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, REAL, SQRT
00151 *     ..
00152 *     .. Statement Functions ..
00153       REAL               CABS1
00154 *     ..
00155 *     .. Statement Function definitions ..
00156       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( 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( 'CHETF2', -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( REAL( 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 = ICAMAX( 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. SISNAN(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 ) = REAL( 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 + ICAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
00231                ROWMAX = CABS1( A( IMAX, JMAX ) )
00232                IF( IMAX.GT.1 ) THEN
00233                   JMAX = ICAMAX( 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( REAL( 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 CSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
00266                DO 20 J = KP + 1, KK - 1
00267                   T = CONJG( A( J, KK ) )
00268                   A( J, KK ) = CONJG( A( KP, J ) )
00269                   A( KP, J ) = T
00270    20          CONTINUE
00271                A( KP, KK ) = CONJG( A( KP, KK ) )
00272                R1 = REAL( A( KK, KK ) )
00273                A( KK, KK ) = REAL( A( KP, KP ) )
00274                A( KP, KP ) = R1
00275                IF( KSTEP.EQ.2 ) THEN
00276                   A( K, K ) = REAL( 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 ) = REAL( A( K, K ) )
00283                IF( KSTEP.EQ.2 )
00284      $            A( K-1, K-1 ) = REAL( 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 / REAL( A( K, K ) )
00302                CALL CHER( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
00303 *
00304 *              Store U(k) in column k
00305 *
00306                CALL CSSCAL( 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 = SLAPY2( REAL( A( K-1, K ) ),
00324      $                AIMAG( A( K-1, K ) ) )
00325                   D22 = REAL( A( K-1, K-1 ) ) / D
00326                   D11 = REAL( 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 )-CONJG( D12 )*A( J, K ) )
00333                      WK = D*( D22*A( J, K )-D12*A( J, K-1 ) )
00334                      DO 30 I = J, 1, -1
00335                         A( I, J ) = A( I, J ) - A( I, K )*CONJG( WK ) -
00336      $                              A( I, K-1 )*CONJG( WKM1 )
00337    30                CONTINUE
00338                      A( J, K ) = WK
00339                      A( J, K-1 ) = WKM1
00340                      A( J, J ) = CMPLX( REAL( A( J, J ) ), 0.0E+0 )
00341    40             CONTINUE
00342 *
00343                END IF
00344 *
00345             END IF
00346          END IF
00347 *
00348 *        Store details of the interchanges in IPIV
00349 *
00350          IF( KSTEP.EQ.1 ) THEN
00351             IPIV( K ) = KP
00352          ELSE
00353             IPIV( K ) = -KP
00354             IPIV( K-1 ) = -KP
00355          END IF
00356 *
00357 *        Decrease K and return to the start of the main loop
00358 *
00359          K = K - KSTEP
00360          GO TO 10
00361 *
00362       ELSE
00363 *
00364 *        Factorize A as L*D*L' using the lower triangle of A
00365 *
00366 *        K is the main loop index, increasing from 1 to N in steps of
00367 *        1 or 2
00368 *
00369          K = 1
00370    50    CONTINUE
00371 *
00372 *        If K > N, exit from loop
00373 *
00374          IF( K.GT.N )
00375      $      GO TO 90
00376          KSTEP = 1
00377 *
00378 *        Determine rows and columns to be interchanged and whether
00379 *        a 1-by-1 or 2-by-2 pivot block will be used
00380 *
00381          ABSAKK = ABS( REAL( A( K, K ) ) )
00382 *
00383 *        IMAX is the row-index of the largest off-diagonal element in
00384 *        column K, and COLMAX is its absolute value
00385 *
00386          IF( K.LT.N ) THEN
00387             IMAX = K + ICAMAX( N-K, A( K+1, K ), 1 )
00388             COLMAX = CABS1( A( IMAX, K ) )
00389          ELSE
00390             COLMAX = ZERO
00391          END IF
00392 *
00393          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN
00394 *
00395 *           Column K is zero or contains a NaN: set INFO and continue
00396 *
00397             IF( INFO.EQ.0 )
00398      $         INFO = K
00399             KP = K
00400             A( K, K ) = REAL( A( K, K ) )
00401          ELSE
00402             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
00403 *
00404 *              no interchange, use 1-by-1 pivot block
00405 *
00406                KP = K
00407             ELSE
00408 *
00409 *              JMAX is the column-index of the largest off-diagonal
00410 *              element in row IMAX, and ROWMAX is its absolute value
00411 *
00412                JMAX = K - 1 + ICAMAX( IMAX-K, A( IMAX, K ), LDA )
00413                ROWMAX = CABS1( A( IMAX, JMAX ) )
00414                IF( IMAX.LT.N ) THEN
00415                   JMAX = IMAX + ICAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
00416                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
00417                END IF
00418 *
00419                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
00420 *
00421 *                 no interchange, use 1-by-1 pivot block
00422 *
00423                   KP = K
00424                ELSE IF( ABS( REAL( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
00425      $                   THEN
00426 *
00427 *                 interchange rows and columns K and IMAX, use 1-by-1
00428 *                 pivot block
00429 *
00430                   KP = IMAX
00431                ELSE
00432 *
00433 *                 interchange rows and columns K+1 and IMAX, use 2-by-2
00434 *                 pivot block
00435 *
00436                   KP = IMAX
00437                   KSTEP = 2
00438                END IF
00439             END IF
00440 *
00441             KK = K + KSTEP - 1
00442             IF( KP.NE.KK ) THEN
00443 *
00444 *              Interchange rows and columns KK and KP in the trailing
00445 *              submatrix A(k:n,k:n)
00446 *
00447                IF( KP.LT.N )
00448      $            CALL CSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
00449                DO 60 J = KK + 1, KP - 1
00450                   T = CONJG( A( J, KK ) )
00451                   A( J, KK ) = CONJG( A( KP, J ) )
00452                   A( KP, J ) = T
00453    60          CONTINUE
00454                A( KP, KK ) = CONJG( A( KP, KK ) )
00455                R1 = REAL( A( KK, KK ) )
00456                A( KK, KK ) = REAL( A( KP, KP ) )
00457                A( KP, KP ) = R1
00458                IF( KSTEP.EQ.2 ) THEN
00459                   A( K, K ) = REAL( A( K, K ) )
00460                   T = A( K+1, K )
00461                   A( K+1, K ) = A( KP, K )
00462                   A( KP, K ) = T
00463                END IF
00464             ELSE
00465                A( K, K ) = REAL( A( K, K ) )
00466                IF( KSTEP.EQ.2 )
00467      $            A( K+1, K+1 ) = REAL( A( K+1, K+1 ) )
00468             END IF
00469 *
00470 *           Update the trailing submatrix
00471 *
00472             IF( KSTEP.EQ.1 ) THEN
00473 *
00474 *              1-by-1 pivot block D(k): column k now holds
00475 *
00476 *              W(k) = L(k)*D(k)
00477 *
00478 *              where L(k) is the k-th column of L
00479 *
00480                IF( K.LT.N ) THEN
00481 *
00482 *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
00483 *
00484 *                 A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)'
00485 *
00486                   R1 = ONE / REAL( A( K, K ) )
00487                   CALL CHER( UPLO, N-K, -R1, A( K+1, K ), 1,
00488      $                       A( K+1, K+1 ), LDA )
00489 *
00490 *                 Store L(k) in column K
00491 *
00492                   CALL CSSCAL( N-K, R1, A( K+1, K ), 1 )
00493                END IF
00494             ELSE
00495 *
00496 *              2-by-2 pivot block D(k)
00497 *
00498                IF( K.LT.N-1 ) THEN
00499 *
00500 *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
00501 *
00502 *                 A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )'
00503 *                    = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )'
00504 *
00505 *                 where L(k) and L(k+1) are the k-th and (k+1)-th
00506 *                 columns of L
00507 *
00508                   D = SLAPY2( REAL( A( K+1, K ) ),
00509      $                        AIMAG( A( K+1, K ) ) )
00510                   D11 = REAL( A( K+1, K+1 ) ) / D
00511                   D22 = REAL( A( K, K ) ) / D
00512                   TT = ONE / ( D11*D22-ONE )
00513                   D21 = A( K+1, K ) / D
00514                   D =  TT / D
00515 *
00516                   DO 80 J = K + 2, N
00517                      WK = D*( D11*A( J, K )-D21*A( J, K+1 ) )
00518                      WKP1 = D*( D22*A( J, K+1 )-CONJG( D21 )*A( J, K ) )
00519                      DO 70 I = J, N
00520                         A( I, J ) = A( I, J ) - A( I, K )*CONJG( WK ) -
00521      $                              A( I, K+1 )*CONJG( WKP1 )
00522    70                CONTINUE
00523                      A( J, K ) = WK
00524                      A( J, K+1 ) = WKP1
00525                      A( J, J ) = CMPLX( REAL( A( J, J ) ), 0.0E+0 )
00526    80             CONTINUE
00527                END IF
00528             END IF
00529          END IF
00530 *
00531 *        Store details of the interchanges in IPIV
00532 *
00533          IF( KSTEP.EQ.1 ) THEN
00534             IPIV( K ) = KP
00535          ELSE
00536             IPIV( K ) = -KP
00537             IPIV( K+1 ) = -KP
00538          END IF
00539 *
00540 *        Increase K and return to the start of the main loop
00541 *
00542          K = K + KSTEP
00543          GO TO 50
00544 *
00545       END IF
00546 *
00547    90 CONTINUE
00548       RETURN
00549 *
00550 *     End of CHETF2
00551 *
00552       END