LAPACK 3.3.0

csytf2.f

Go to the documentation of this file.
00001       SUBROUTINE CSYTF2( 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 *  CSYTF2 computes the factorization of a complex symmetric 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 transpose of U, and D is symmetric and
00027 *  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 *          symmetric 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 symmetric 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.209 and l.377
00083 *         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
00084 *    by
00085 *         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN
00086 *
00087 *  1-96 - Based on modifications by J. Lewis, Boeing Computer Services
00088 *         Company
00089 *
00090 *  If UPLO = 'U', then A = U*D*U', where
00091 *     U = P(n)*U(n)* ... *P(k)U(k)* ...,
00092 *  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
00093 *  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
00094 *  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
00095 *  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
00096 *  that if the diagonal block D(k) is of order s (s = 1 or 2), then
00097 *
00098 *             (   I    v    0   )   k-s
00099 *     U(k) =  (   0    I    0   )   s
00100 *             (   0    0    I   )   n-k
00101 *                k-s   s   n-k
00102 *
00103 *  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
00104 *  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
00105 *  and A(k,k), and v overwrites A(1:k-2,k-1:k).
00106 *
00107 *  If UPLO = 'L', then A = L*D*L', where
00108 *     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
00109 *  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
00110 *  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
00111 *  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
00112 *  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
00113 *  that if the diagonal block D(k) is of order s (s = 1 or 2), then
00114 *
00115 *             (   I    0     0   )  k-1
00116 *     L(k) =  (   0    I     0   )  s
00117 *             (   0    v     I   )  n-k-s+1
00118 *                k-1   s  n-k-s+1
00119 *
00120 *  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
00121 *  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
00122 *  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
00123 *
00124 *  =====================================================================
00125 *
00126 *     .. Parameters ..
00127       REAL               ZERO, ONE
00128       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00129       REAL               EIGHT, SEVTEN
00130       PARAMETER          ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
00131       COMPLEX            CONE
00132       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
00133 *     ..
00134 *     .. Local Scalars ..
00135       LOGICAL            UPPER
00136       INTEGER            I, IMAX, J, JMAX, K, KK, KP, KSTEP
00137       REAL               ABSAKK, ALPHA, COLMAX, ROWMAX
00138       COMPLEX            D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, Z
00139 *     ..
00140 *     .. External Functions ..
00141       LOGICAL            LSAME, SISNAN
00142       INTEGER            ICAMAX
00143       EXTERNAL           LSAME, ICAMAX, SISNAN
00144 *     ..
00145 *     .. External Subroutines ..
00146       EXTERNAL           CSCAL, CSWAP, CSYR, XERBLA
00147 *     ..
00148 *     .. Intrinsic Functions ..
00149       INTRINSIC          ABS, AIMAG, MAX, REAL, SQRT
00150 *     ..
00151 *     .. Statement Functions ..
00152       REAL               CABS1
00153 *     ..
00154 *     .. Statement Function definitions ..
00155       CABS1( Z ) = ABS( REAL( Z ) ) + ABS( AIMAG( Z ) )
00156 *     ..
00157 *     .. Executable Statements ..
00158 *
00159 *     Test the input parameters.
00160 *
00161       INFO = 0
00162       UPPER = LSAME( UPLO, 'U' )
00163       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00164          INFO = -1
00165       ELSE IF( N.LT.0 ) THEN
00166          INFO = -2
00167       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00168          INFO = -4
00169       END IF
00170       IF( INFO.NE.0 ) THEN
00171          CALL XERBLA( 'CSYTF2', -INFO )
00172          RETURN
00173       END IF
00174 *
00175 *     Initialize ALPHA for use in choosing pivot block size.
00176 *
00177       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
00178 *
00179       IF( UPPER ) THEN
00180 *
00181 *        Factorize A as U*D*U' using the upper triangle of A
00182 *
00183 *        K is the main loop index, decreasing from N to 1 in steps of
00184 *        1 or 2
00185 *
00186          K = N
00187    10    CONTINUE
00188 *
00189 *        If K < 1, exit from loop
00190 *
00191          IF( K.LT.1 )
00192      $      GO TO 70
00193          KSTEP = 1
00194 *
00195 *        Determine rows and columns to be interchanged and whether
00196 *        a 1-by-1 or 2-by-2 pivot block will be used
00197 *
00198          ABSAKK = CABS1( A( K, K ) )
00199 *
00200 *        IMAX is the row-index of the largest off-diagonal element in
00201 *        column K, and COLMAX is its absolute value
00202 *
00203          IF( K.GT.1 ) THEN
00204             IMAX = ICAMAX( K-1, A( 1, K ), 1 )
00205             COLMAX = CABS1( A( IMAX, K ) )
00206          ELSE
00207             COLMAX = ZERO
00208          END IF
00209 *
00210          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. SISNAN(ABSAKK) ) THEN
00211 *
00212 *           Column K is zero or NaN: set INFO and continue
00213 *
00214             IF( INFO.EQ.0 )
00215      $         INFO = K
00216             KP = K
00217          ELSE
00218             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
00219 *
00220 *              no interchange, use 1-by-1 pivot block
00221 *
00222                KP = K
00223             ELSE
00224 *
00225 *              JMAX is the column-index of the largest off-diagonal
00226 *              element in row IMAX, and ROWMAX is its absolute value
00227 *
00228                JMAX = IMAX + ICAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
00229                ROWMAX = CABS1( A( IMAX, JMAX ) )
00230                IF( IMAX.GT.1 ) THEN
00231                   JMAX = ICAMAX( IMAX-1, A( 1, IMAX ), 1 )
00232                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
00233                END IF
00234 *
00235                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
00236 *
00237 *                 no interchange, use 1-by-1 pivot block
00238 *
00239                   KP = K
00240                ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
00241 *
00242 *                 interchange rows and columns K and IMAX, use 1-by-1
00243 *                 pivot block
00244 *
00245                   KP = IMAX
00246                ELSE
00247 *
00248 *                 interchange rows and columns K-1 and IMAX, use 2-by-2
00249 *                 pivot block
00250 *
00251                   KP = IMAX
00252                   KSTEP = 2
00253                END IF
00254             END IF
00255 *
00256             KK = K - KSTEP + 1
00257             IF( KP.NE.KK ) THEN
00258 *
00259 *              Interchange rows and columns KK and KP in the leading
00260 *              submatrix A(1:k,1:k)
00261 *
00262                CALL CSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
00263                CALL CSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ),
00264      $                     LDA )
00265                T = A( KK, KK )
00266                A( KK, KK ) = A( KP, KP )
00267                A( KP, KP ) = T
00268                IF( KSTEP.EQ.2 ) THEN
00269                   T = A( K-1, K )
00270                   A( K-1, K ) = A( KP, K )
00271                   A( KP, K ) = T
00272                END IF
00273             END IF
00274 *
00275 *           Update the leading submatrix
00276 *
00277             IF( KSTEP.EQ.1 ) THEN
00278 *
00279 *              1-by-1 pivot block D(k): column k now holds
00280 *
00281 *              W(k) = U(k)*D(k)
00282 *
00283 *              where U(k) is the k-th column of U
00284 *
00285 *              Perform a rank-1 update of A(1:k-1,1:k-1) as
00286 *
00287 *              A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)'
00288 *
00289                R1 = CONE / A( K, K )
00290                CALL CSYR( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
00291 *
00292 *              Store U(k) in column k
00293 *
00294                CALL CSCAL( K-1, R1, A( 1, K ), 1 )
00295             ELSE
00296 *
00297 *              2-by-2 pivot block D(k): columns k and k-1 now hold
00298 *
00299 *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
00300 *
00301 *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
00302 *              of U
00303 *
00304 *              Perform a rank-2 update of A(1:k-2,1:k-2) as
00305 *
00306 *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )'
00307 *                 = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )'
00308 *
00309                IF( K.GT.2 ) THEN
00310 *
00311                   D12 = A( K-1, K )
00312                   D22 = A( K-1, K-1 ) / D12
00313                   D11 = A( K, K ) / D12
00314                   T = CONE / ( D11*D22-CONE )
00315                   D12 = T / D12
00316 *
00317                   DO 30 J = K - 2, 1, -1
00318                      WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) )
00319                      WK = D12*( D22*A( J, K )-A( J, K-1 ) )
00320                      DO 20 I = J, 1, -1
00321                         A( I, J ) = A( I, J ) - A( I, K )*WK -
00322      $                              A( I, K-1 )*WKM1
00323    20                CONTINUE
00324                      A( J, K ) = WK
00325                      A( J, K-1 ) = WKM1
00326    30             CONTINUE
00327 *
00328                END IF
00329 *
00330             END IF
00331          END IF
00332 *
00333 *        Store details of the interchanges in IPIV
00334 *
00335          IF( KSTEP.EQ.1 ) THEN
00336             IPIV( K ) = KP
00337          ELSE
00338             IPIV( K ) = -KP
00339             IPIV( K-1 ) = -KP
00340          END IF
00341 *
00342 *        Decrease K and return to the start of the main loop
00343 *
00344          K = K - KSTEP
00345          GO TO 10
00346 *
00347       ELSE
00348 *
00349 *        Factorize A as L*D*L' using the lower triangle of A
00350 *
00351 *        K is the main loop index, increasing from 1 to N in steps of
00352 *        1 or 2
00353 *
00354          K = 1
00355    40    CONTINUE
00356 *
00357 *        If K > N, exit from loop
00358 *
00359          IF( K.GT.N )
00360      $      GO TO 70
00361          KSTEP = 1
00362 *
00363 *        Determine rows and columns to be interchanged and whether
00364 *        a 1-by-1 or 2-by-2 pivot block will be used
00365 *
00366          ABSAKK = CABS1( A( K, K ) )
00367 *
00368 *        IMAX is the row-index of the largest off-diagonal element in
00369 *        column K, and COLMAX is its absolute value
00370 *
00371          IF( K.LT.N ) THEN
00372             IMAX = K + ICAMAX( N-K, A( K+1, K ), 1 )
00373             COLMAX = CABS1( A( IMAX, K ) )
00374          ELSE
00375             COLMAX = ZERO
00376          END IF
00377 *
00378          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. SISNAN(ABSAKK) ) THEN
00379 *
00380 *           Column K is zero or NaN: set INFO and continue
00381 *
00382             IF( INFO.EQ.0 )
00383      $         INFO = K
00384             KP = K
00385          ELSE
00386             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
00387 *
00388 *              no interchange, use 1-by-1 pivot block
00389 *
00390                KP = K
00391             ELSE
00392 *
00393 *              JMAX is the column-index of the largest off-diagonal
00394 *              element in row IMAX, and ROWMAX is its absolute value
00395 *
00396                JMAX = K - 1 + ICAMAX( IMAX-K, A( IMAX, K ), LDA )
00397                ROWMAX = CABS1( A( IMAX, JMAX ) )
00398                IF( IMAX.LT.N ) THEN
00399                   JMAX = IMAX + ICAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
00400                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
00401                END IF
00402 *
00403                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
00404 *
00405 *                 no interchange, use 1-by-1 pivot block
00406 *
00407                   KP = K
00408                ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
00409 *
00410 *                 interchange rows and columns K and IMAX, use 1-by-1
00411 *                 pivot block
00412 *
00413                   KP = IMAX
00414                ELSE
00415 *
00416 *                 interchange rows and columns K+1 and IMAX, use 2-by-2
00417 *                 pivot block
00418 *
00419                   KP = IMAX
00420                   KSTEP = 2
00421                END IF
00422             END IF
00423 *
00424             KK = K + KSTEP - 1
00425             IF( KP.NE.KK ) THEN
00426 *
00427 *              Interchange rows and columns KK and KP in the trailing
00428 *              submatrix A(k:n,k:n)
00429 *
00430                IF( KP.LT.N )
00431      $            CALL CSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
00432                CALL CSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
00433      $                     LDA )
00434                T = A( KK, KK )
00435                A( KK, KK ) = A( KP, KP )
00436                A( KP, KP ) = T
00437                IF( KSTEP.EQ.2 ) THEN
00438                   T = A( K+1, K )
00439                   A( K+1, K ) = A( KP, K )
00440                   A( KP, K ) = T
00441                END IF
00442             END IF
00443 *
00444 *           Update the trailing submatrix
00445 *
00446             IF( KSTEP.EQ.1 ) THEN
00447 *
00448 *              1-by-1 pivot block D(k): column k now holds
00449 *
00450 *              W(k) = L(k)*D(k)
00451 *
00452 *              where L(k) is the k-th column of L
00453 *
00454                IF( K.LT.N ) THEN
00455 *
00456 *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
00457 *
00458 *                 A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)'
00459 *
00460                   R1 = CONE / A( K, K )
00461                   CALL CSYR( UPLO, N-K, -R1, A( K+1, K ), 1,
00462      $                       A( K+1, K+1 ), LDA )
00463 *
00464 *                 Store L(k) in column K
00465 *
00466                   CALL CSCAL( N-K, R1, A( K+1, K ), 1 )
00467                END IF
00468             ELSE
00469 *
00470 *              2-by-2 pivot block D(k)
00471 *
00472                IF( K.LT.N-1 ) THEN
00473 *
00474 *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
00475 *
00476 *                 A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )'
00477 *                    = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )'
00478 *
00479 *                 where L(k) and L(k+1) are the k-th and (k+1)-th
00480 *                 columns of L
00481 *
00482                   D21 = A( K+1, K )
00483                   D11 = A( K+1, K+1 ) / D21
00484                   D22 = A( K, K ) / D21
00485                   T = CONE / ( D11*D22-CONE )
00486                   D21 = T / D21
00487 *
00488                   DO 60 J = K + 2, N
00489                      WK = D21*( D11*A( J, K )-A( J, K+1 ) )
00490                      WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) )
00491                      DO 50 I = J, N
00492                         A( I, J ) = A( I, J ) - A( I, K )*WK -
00493      $                              A( I, K+1 )*WKP1
00494    50                CONTINUE
00495                      A( J, K ) = WK
00496                      A( J, K+1 ) = WKP1
00497    60             CONTINUE
00498                END IF
00499             END IF
00500          END IF
00501 *
00502 *        Store details of the interchanges in IPIV
00503 *
00504          IF( KSTEP.EQ.1 ) THEN
00505             IPIV( K ) = KP
00506          ELSE
00507             IPIV( K ) = -KP
00508             IPIV( K+1 ) = -KP
00509          END IF
00510 *
00511 *        Increase K and return to the start of the main loop
00512 *
00513          K = K + KSTEP
00514          GO TO 40
00515 *
00516       END IF
00517 *
00518    70 CONTINUE
00519       RETURN
00520 *
00521 *     End of CSYTF2
00522 *
00523       END
 All Files Functions