Annotation of rpl/lapack/lapack/zgtcon.f, revision 1.8

1.1       bertrand    1:       SUBROUTINE ZGTCON( NORM, N, DL, D, DU, DU2, IPIV, ANORM, RCOND,
                      2:      $                   WORK, INFO )
                      3: *
1.8     ! bertrand    4: *  -- LAPACK routine (version 3.3.1) --
1.1       bertrand    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                      6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8     ! bertrand    7: *  -- April 2011                                                      --
1.1       bertrand    8: *
                      9: *     Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
                     10: *
                     11: *     .. Scalar Arguments ..
                     12:       CHARACTER          NORM
                     13:       INTEGER            INFO, N
                     14:       DOUBLE PRECISION   ANORM, RCOND
                     15: *     ..
                     16: *     .. Array Arguments ..
                     17:       INTEGER            IPIV( * )
                     18:       COMPLEX*16         D( * ), DL( * ), DU( * ), DU2( * ), WORK( * )
                     19: *     ..
                     20: *
                     21: *  Purpose
                     22: *  =======
                     23: *
                     24: *  ZGTCON estimates the reciprocal of the condition number of a complex
                     25: *  tridiagonal matrix A using the LU factorization as computed by
                     26: *  ZGTTRF.
                     27: *
                     28: *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
                     29: *  condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
                     30: *
                     31: *  Arguments
                     32: *  =========
                     33: *
                     34: *  NORM    (input) CHARACTER*1
                     35: *          Specifies whether the 1-norm condition number or the
                     36: *          infinity-norm condition number is required:
                     37: *          = '1' or 'O':  1-norm;
                     38: *          = 'I':         Infinity-norm.
                     39: *
                     40: *  N       (input) INTEGER
                     41: *          The order of the matrix A.  N >= 0.
                     42: *
                     43: *  DL      (input) COMPLEX*16 array, dimension (N-1)
                     44: *          The (n-1) multipliers that define the matrix L from the
                     45: *          LU factorization of A as computed by ZGTTRF.
                     46: *
                     47: *  D       (input) COMPLEX*16 array, dimension (N)
                     48: *          The n diagonal elements of the upper triangular matrix U from
                     49: *          the LU factorization of A.
                     50: *
                     51: *  DU      (input) COMPLEX*16 array, dimension (N-1)
                     52: *          The (n-1) elements of the first superdiagonal of U.
                     53: *
                     54: *  DU2     (input) COMPLEX*16 array, dimension (N-2)
                     55: *          The (n-2) elements of the second superdiagonal of U.
                     56: *
                     57: *  IPIV    (input) INTEGER array, dimension (N)
                     58: *          The pivot indices; for 1 <= i <= n, row i of the matrix was
                     59: *          interchanged with row IPIV(i).  IPIV(i) will always be either
                     60: *          i or i+1; IPIV(i) = i indicates a row interchange was not
                     61: *          required.
                     62: *
                     63: *  ANORM   (input) DOUBLE PRECISION
                     64: *          If NORM = '1' or 'O', the 1-norm of the original matrix A.
                     65: *          If NORM = 'I', the infinity-norm of the original matrix A.
                     66: *
                     67: *  RCOND   (output) DOUBLE PRECISION
                     68: *          The reciprocal of the condition number of the matrix A,
                     69: *          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
                     70: *          estimate of the 1-norm of inv(A) computed in this routine.
                     71: *
                     72: *  WORK    (workspace) COMPLEX*16 array, dimension (2*N)
                     73: *
                     74: *  INFO    (output) INTEGER
                     75: *          = 0:  successful exit
                     76: *          < 0:  if INFO = -i, the i-th argument had an illegal value
                     77: *
                     78: *  =====================================================================
                     79: *
                     80: *     .. Parameters ..
                     81:       DOUBLE PRECISION   ONE, ZERO
                     82:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
                     83: *     ..
                     84: *     .. Local Scalars ..
                     85:       LOGICAL            ONENRM
                     86:       INTEGER            I, KASE, KASE1
                     87:       DOUBLE PRECISION   AINVNM
                     88: *     ..
                     89: *     .. Local Arrays ..
                     90:       INTEGER            ISAVE( 3 )
                     91: *     ..
                     92: *     .. External Functions ..
                     93:       LOGICAL            LSAME
                     94:       EXTERNAL           LSAME
                     95: *     ..
                     96: *     .. External Subroutines ..
                     97:       EXTERNAL           XERBLA, ZGTTRS, ZLACN2
                     98: *     ..
                     99: *     .. Intrinsic Functions ..
                    100:       INTRINSIC          DCMPLX
                    101: *     ..
                    102: *     .. Executable Statements ..
                    103: *
                    104: *     Test the input arguments.
                    105: *
                    106:       INFO = 0
                    107:       ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
                    108:       IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
                    109:          INFO = -1
                    110:       ELSE IF( N.LT.0 ) THEN
                    111:          INFO = -2
                    112:       ELSE IF( ANORM.LT.ZERO ) THEN
                    113:          INFO = -8
                    114:       END IF
                    115:       IF( INFO.NE.0 ) THEN
                    116:          CALL XERBLA( 'ZGTCON', -INFO )
                    117:          RETURN
                    118:       END IF
                    119: *
                    120: *     Quick return if possible
                    121: *
                    122:       RCOND = ZERO
                    123:       IF( N.EQ.0 ) THEN
                    124:          RCOND = ONE
                    125:          RETURN
                    126:       ELSE IF( ANORM.EQ.ZERO ) THEN
                    127:          RETURN
                    128:       END IF
                    129: *
                    130: *     Check that D(1:N) is non-zero.
                    131: *
                    132:       DO 10 I = 1, N
                    133:          IF( D( I ).EQ.DCMPLX( ZERO ) )
                    134:      $      RETURN
                    135:    10 CONTINUE
                    136: *
                    137:       AINVNM = ZERO
                    138:       IF( ONENRM ) THEN
                    139:          KASE1 = 1
                    140:       ELSE
                    141:          KASE1 = 2
                    142:       END IF
                    143:       KASE = 0
                    144:    20 CONTINUE
                    145:       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
                    146:       IF( KASE.NE.0 ) THEN
                    147:          IF( KASE.EQ.KASE1 ) THEN
                    148: *
                    149: *           Multiply by inv(U)*inv(L).
                    150: *
                    151:             CALL ZGTTRS( 'No transpose', N, 1, DL, D, DU, DU2, IPIV,
                    152:      $                   WORK, N, INFO )
                    153:          ELSE
                    154: *
1.8     ! bertrand  155: *           Multiply by inv(L**H)*inv(U**H).
1.1       bertrand  156: *
                    157:             CALL ZGTTRS( 'Conjugate transpose', N, 1, DL, D, DU, DU2,
                    158:      $                   IPIV, WORK, N, INFO )
                    159:          END IF
                    160:          GO TO 20
                    161:       END IF
                    162: *
                    163: *     Compute the estimate of the reciprocal condition number.
                    164: *
                    165:       IF( AINVNM.NE.ZERO )
                    166:      $   RCOND = ( ONE / AINVNM ) / ANORM
                    167: *
                    168:       RETURN
                    169: *
                    170: *     End of ZGTCON
                    171: *
                    172:       END

CVSweb interface <joel.bertrand@systella.fr>